備忘録_数値微分

備忘録_数値微分

微分

微分とはある瞬間の変化の量を表したものである

数式で表すと・・・

     \[ \frac{df(x)}{dx} = \lim_{h \to 0} \frac{f(x+h)-f(x)}{h}\]

  • \frac{df(x)}{dx}
    f(x)xについてのxに対するf(x)(微分)を表す記号
  • \displaystyle\lim_{h \to 0}
    「小さな変化」であるhを限りなく0に近づけるという意味

導関数と微分|数学|苦手解決Q&A|進研ゼミ高校講座|ベネッセコーポレーションにわかりやすく書いてあります。

そのままPythonで実装すると・・・

def numerical_diff(f, x):
    h = 10e-50
    return (f(x+h) - f(x)) / h
  • x
    関数fへの引数
  • h
    できる限り小さな値
  • 10e-50
    「0.00・・・1」の0が50個続く数

丸め誤差が問題になる

丸め誤差とは、少数の小さな範囲において数値が省略されることで、最終的な計算結果に誤差が生じることを言う。

Pythonでの丸め誤差の例

np.float32(1e-50)
0.0
真の微分と数値微分の値は異なる為問題になる
x+hとxの間での関数fの差分を計算しているが、そもそも、この計算には誤差が生じる。「真の微分」は、xの位置での関数の傾きに対応するが、今回の実装で行っている微分は(x+h)xの間の傾きに対応する。そのため、真の微分と今回の実装の値は厳密には一致しない。この差はhを無限に0へと近づけることができない為である。
  1. 丸め誤差の問題の対応
    微小な値hとして10^{-4}を用いる。
  2. 真の微分と数値微分の値は異なる為問題の対応
    xを中心してその前後の差分を計算することで、誤差を減らす。(中心差分)

2つの問題を対応した形で実装すると・・・

def numerical_diff(f, x):
    h = 1e-4 #0.0001
    return (f(x+h) - f(x-h)) / (2*h)

数値微分の例

今回の例で使用する2次関数

     \[ x = 0.01x^{2}+0.1x \]

上記の2次関数をPythonで実装

def function_1(x):
    return 0.01*x**2 + 0.1*x
import numpy as np

グラフの生成

import numpy as np
import matplotlib.pylab as plt
x = np.arange(0.0, 20.0, 0.1)
y=function_1(x)
plt.xlabel("x")
plt.ylabel("f(x)")
plt.plot(x, y)
plt.show()

x=5、x=10のときの微分を計算

今回の微分→xに対するf(x)の変化の量→関数の傾き

numerical_diff(function_1, 5)
0.1999999999990898
numerical_diff(function_1, 10)
0.2999999999986347

尚、普通に微分で解くと、0.2、0.3になると思うので、だいたい合っていることがわかる。

偏微分(へんびぶん)

偏微分とは複数の変数からなる関数の微分である

今回の例で使用する関数

     \[ f(x^{}_{0},x^{}_{1})=x^{2}_{0}+x^{2}_{1} \]

偏微分を数式で表すと・・・

     \[ \frac{\partial f}{x^{}_{0}} \qquad \frac{\partial f}{x^{}_{1}} \]

Pythonで実装した場合

def function_2(x):
    return x[0]**2 + x[1]**2

グラフの生成

import numpy as np
import matplotlib.pylab as plt
from mpl_toolkits.mplot3d import Axes3D # 3次元グラフを書く際、使用。

x = np.meshgrid(np.arange(-3, 3, 0.1), np.arange(-3, 3, 0.1))
z = function_2(x)
fig = plt.figure()
ax = Axes3D(fig)
ax.plot_wireframe(x[0], x[1], z)

# 表示範囲の設定
ax.set_xlim(-3.5, 3.5)
ax.set_ylim(-3.5, 3.5)

# 軸ラベルの設定
ax.set_xlabel("x0")
ax.set_ylabel("x1")
ax.set_zlabel("fx")
plt.show()
  • np.meshgrid
    各座標の要素列から格子座標を作成する
  • fig = plt.figure()、ax = Axes3D(fig)
    二次元のグラフを作成し、その後、三次元のグラフを作成
  • ax.plot_wireframe(x[0], x[1], z)
    プロットする

 

偏微分を求める

変数がひとつだけの関数を定義して、その関数について微分を求めるようこと偏微分を求めることできる。(今回は変数以外を特定の値に固定するために、新しい関数を定義)

x^{}{0} = 3x^{}{1} = 4のときのx^{}{0} に対する偏微分

def function_tmp1(x0):
    return x0*x0 + 4.0**2.0

numerical_diff(function_tmp1,3.0)
6.00000000000378

x^{}{0} = 3x^{}{1} = 4のときのx^{}{1} に対する偏微分

def function_tmp2(x1):
    return 3.0**2.0 + x1*x1
numerical_diff(function_tmp2,4.0)
7.999999999999119