微分

　数値積分の反対は数値微分です。数値積分よりも数値微分の計算は簡単です。

　導関数の標準的な定義は次のとおりです。

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

　[0,3] の間の \(f(x) = -x^4 + 4x^3 + 2\) について考えてみましょう。

import matplotlib.pyplot as plt
import numpy as np

plt.rcParams['font.family'] = "Hiragino sans" #for Mac
#plt.rcParams['font.family'] = "MS Gothic" #for Windows

#define the function
f = lambda x: -x**4 + 4*x**3 + 2

#define the parameters
x_min, x_max = 0, 3
npoints = 100
x0 = (x_min+x_max)/2
deltax = 0.5

#draw the plot
x = np.linspace(x_min, x_max, npoints)
plt.plot(x,f(x))

#draw the derivative
x_1 = np.array([x0-deltax, x0, x0+deltax])
plt.plot(x_1,f(x_1),'-o')

for point in x_1:
    linex, liney = [point, point], [0, f(point)]
    plt.plot(linex, liney, color='black', linewidth=2.0)

plt.text(x_1[0]+deltax/2, f(x_1[0]+deltax/2)+1, r'後方差分', {'ha': 'left', 'va': 'bottom'}, rotation=90)
plt.text(x_1[1]+deltax/2, f(x_1[1]+deltax/2)+1, r'前方差分', {'ha': 'left', 'va': 'bottom'}, rotation=90)

plt.xlabel('x')
plt.ylabel('$f(x)$')
plt.xlim([0,3])
plt.show()

前方差分, 後方差分, 中央差分

　点 \(x\) で数値微分することを考えます。プログラミングで数値微分では、微分を差分で近似します。

前方差分

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

後方差分

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

　通常、前方差分と後方差分はほぼ同じ答えを示し、多くの場合どちらを使用してもかまいませんが、ほとんどの場合で前方差分を使用します。

中央差分

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

　点 \(x\) で関数の導関数に不連続がある場合は注意が必要です。

演習：前方差分と後方差分のプログラムを作成せよ

import numpy as np
import matplotlib.pyplot as plt

#define the function
f = lambda x: -x**4 + 4*x**3 + 2

#-----------------------------------------------------
#define the derivative(f, x, h) here
#-----------------------------------------------------

h = 0.001
xrange = np.arange(0, 3, h)
data = []

for x in xrange:
  data.append(derivative(f, x, h))

plt.plot(xrange, data)
#plt.plot(xrange, xrange**3 - 2*xrange**2 +1) # 解析的に求めた例も追加した
plt.show()

２階微分

　関数 \(f(x)\) の 2次導関数の数値近似を導出することもできます。2次導関数は1次導関数の導関数であるため、1次導関数の式を2回適用することで計算できます。

　たとえば、中心差分の式から始めて、\(x + h/2\) と \(x − h/2\) で一次導関数の式を次のように書くことができます。

\[f'(x+h/2) \approx \frac{f(x+h)-f(x)}{h}\]

\[f'(x-h/2) \approx \frac{f(x)-f(x-h)}{h}\]

\[f"(x) \approx \frac{f(x+h)-2f(x)+f(x-h)}{h^2}\]

課題：二次導関数を計算するプログラムを作成せよ

ばらつきが大きいデータの微分

　グラフにプロットすると下の図のように見える量の測定値があるとします。

import numpy as np
import matplotlib.pyplot as plt

x,y = np.loadtxt('noise.dat', unpack=True)
#data = np.column_stack([x,y])
#np.savetxt('noise.dat', data)

plt.plot(x,y)
plt.show()

　実験で得られるデータは、このようにノイズが含まれます。曲線の全体的な形状は分かりますが滑らかではありません。

　ここで、この曲線の一次導関数を計算したい場合、各点での前方差分を計算し、取得した値をプロットするプログラムを作成します。

def derivative(x, h):
    res = x[1:] - x[:-1]
    return res/(2*h)

#for x in xrange:
#  data.append(derivative(func, x, h))

#誤差が多い方法
yd = derivative(y, x[1] - x[0])

plt.plot(yd)
plt.show()

　みてわかるように、微分をとることでノイズの問題がさらに悪化しました。もうこの曲線の形状を確認することは不可能です。これは数値微分によくある問題です。微分をとる曲線にノイズがある場合、微分をとることによってノイズが大幅に誇張され、結果が役に立たなくなることが多々あります。

　残念ながらこの種の問題は実験で取得したデータでは一般的であり、これが数値微分が数値積分よりも使用されない理由の1つです。この問題を軽減するためにできることはいくつかありますが、それらはすべて結果の精度を低下させることになります。

• h の値を増やす

• 導関数が必要な点の近くのデータの一部に曲線を当てはめる

• 微分する前に他の方法でデータを平滑化する