テーラー展開
テーラー展開可能な関数$f$を考えます。
\begin{align}
f(x) &= \sum_{n=0}^{\infty} \frac{1}{n!}f^{(n)}(a)(x-a)^n \\
f^{(n)}(x) &= \frac{\partial^n f}{\partial x^n}
\end{align}
差分公式
$f(x-mh)$を点$x$近傍でテーラー展開します。（$m$は整数）
f(x-mh) = \sum_{n=0}^{\infty} \frac{(-m)^n}{n!}f^{(n)}(x) h^n
これを$m=1,\cdots M$について$n \leq M$の範囲で展開し、これを行列形式で並べます。
\begin{pmatrix}
f(x) \\ f(x-h) \\ f(x-2h) \\ \vdots \\ f(x-Mh)
\end{pmatrix}
= A_M
\begin{pmatrix}
f(x) \\ f^{(1)}(x)h \\ f^{(2)}(x)h^2 \\ \vdots \\ f^{(M)}(x)h^{M}
\end{pmatrix} + O(h^{M+1})
ここで係数行列$A_M$の各成分は
(A_M)_{ab} = \begin{cases}
\frac{(-a)^b}{b!} \quad (a>0) \\
1 \quad(a=b=0) \\
0 \quad (b>a=0)
\end{cases}
です。
このとき、逆行列を使って
\begin{pmatrix}
f(x) \\ f^{(1)}(x)h \\ f^{(2)}(x)h^2 \\ \vdots \\ f^{(M)}(x)h^{M}
\end{pmatrix}
= A_M^{-1}
\begin{pmatrix}
f(x) \\ f(x-h) \\ f(x-2h) \\ \vdots \\ f(x-Mh)
\end{pmatrix} + O(h^{M+1})
となるので、これを解くことで$f^{(n)}(x)$を$O(h^{M+1-n})$の精度で求めることができます。よって、mathematicaなどを使って$A_M^{-1}$さえ計算できれば、任意のオーダーで微分を差分近似できます。
係数行列
係数行列$A_M$は次のように表せます。
A_M = \begin{pmatrix}
1 & 0 & 0 & 0 & \cdots \\
1 & -1 & (-1)^2 & (-1)^3 & \cdots \\
1 & -2 & (-2)^2 & (-2)^3 & \cdots \\
1 & -3 & (-3)^2 & (-3)^3 & \cdots \\
& & \vdots & & \\
\end{pmatrix}
\begin{pmatrix}
1 \\
& 1 \\
& & 1/2 \\
& & & \ddots \\
& & & & 1/n! \\
& & & & & \ddots
\end{pmatrix}
次のmathematicaコードで$A_n$の逆行列を計算します。
diff_approx.nb
f[a_, b_] = Which[
a == 0 && b == 0, 1,
True, (-a)^b/(b!)
];
n = 3;
Mat = Table[f[i, j], {i, 0, n}, {j, 0, n}];
MatInv = Inverse[Mat];
Mat // MatrixForm
MatInv // MatrixForm
例
M=2
2次のオーダーまでテーラー展開してみます。
\begin{align}
\begin{pmatrix}
f(x) \\ f^{(1)}(x)h \\ f^{(2)}(x)h^2
\end{pmatrix}
&= \begin{pmatrix}
1 & 0 & 0 \\
1 & -1 & 1/2 \\
1 & -2 & 2
\end{pmatrix}^{-1}
\begin{pmatrix}
f(x) \\ f(x-h) \\ f(x-2h)
\end{pmatrix} + O(h^{M+1}) \\
&= \frac{1}{2} \begin{pmatrix}
2 & 0 & 0 \\
3 & -4 & 1 \\
2 & -4 & 2
\end{pmatrix}
\begin{pmatrix}
f(x) \\ f(x-h) \\ f(x-2h)
\end{pmatrix} + O(h^{M+1})
\end{align}
よって
\begin{align}
f^{(1)}(x) &= \frac{3f(x)-4f(x-h)+f(x-2h)}{2h} + O(h^2) \\
f^{(2)}(x) &= \frac{f(x)-2f(x-h)+f(x-2h)}{h^2} + O(h)
\end{align}
と計算できます。
M=3
\begin{align}
A_3 &= \begin{pmatrix}
1 & 0 & 0 & 0 \\
1 & -1 & 1/2 & -1/6 \\
1 & -2 & 2 & -4/3 \\
1 & -3 & 9/2 & -9/2
\end{pmatrix} \\ \\
A_3^{-1} &= \begin{pmatrix}
1 & 0 & 0 & 0 \\
11/6 & -3 & 3/2 & -1/3 \\
2 & -5 & 4 & -1 \\
1 & -3 & 3 & -1
\end{pmatrix}
\end{align}
なので、
\begin{align}
f^{(1)}(x) &= \frac{ 11f(x) -18f(x-h) + 9f(x-2h) - 2f(x-3h) }{6h} + O(h^3)\\
f^{(2)}(x) &= \frac{ 2f(x) -5f(x-h) +4f(x-2h) -f(x-3h) }{h^2} + O(h^2)\\
f^{(3)}(x) &= \frac{ f(x) -3f(x-h) +3f(x-2h) -f(x-3h) }{h^3} + O(h)
\end{align}
です。
M=4
\begin{align}
A_4 &= \begin{pmatrix}
1 & 0 & 0 & 0 & 0 \\
1 & -1 & 1/2 & -1/6 & 1/24 \\
1 & -2 & 2 & -4/3 & 2/3 \\
1 & -3 & 9/2 & -9/2 & 27/8 \\
1 & -4 & 8 & -32/3 & 32/3
\end{pmatrix} \\ \\
A_4^{-1} &= \begin{pmatrix}
1 & 0 & 0 & 0 \\
25/12 & -4 & 3 & -4/3 & 1/4 \\
35/12 & -26/3 & 19/2 & -14/3 & 11/12 \\
5/2 & -9 & 12 & -7 & 3/2 \\
1 & -4 & 6 & -4 & 1
\end{pmatrix}
\end{align}
だから、
f^{(3)}(x) = \frac{ 5f(x) -18f(x-h) +24f(x-2h) -14f(x-3h) +3f(x-4h)}{2h^3} + O(h^2)
です。
大抵$h \sim 10^{-3}$程度で計算することが多いのでこのくらいまでやっておけば３次の導関数までが$\sim O(10^{-6})$以下の誤差に収まって精度よく計算できると思います。
数値実験
$M=4$の場合で$f(x) = \cos(x)$の微分を３次微分まで数値的に計算してみます。$x \in[0,2\pi]$の範囲を1000等分した点列を計算しました。($h=2\pi/1000$)
また、解析解との誤差の最大値はそれぞれ以下の通りです。
\begin{align}
&\frac{df}{dx}: 3.12 \times 10^{-10} \\
&\frac{d^2f}{dx^2}: 2.07 \times 10^{-7} \\
&\frac{d^3f}{dx^3}: 6.911 \times 10^{-5}
\end{align}
$h^2 \sim 4 \times 10^{-5}$なので、おおよそその範囲の誤差が収まっています。(pythonの精度との兼ね合いもありますが)
使用したコード
diff_approx.py
import numpy as np
from matplotlib import pyplot as plt
def diff(func, x, h):
f0 = (25/12)*func(x)
f1 = -4*func(x-h)
f2 = 3*func(x-2*h)
f3 = -(4/3)*func(x-3*h)
f4 = (1/4)*func(x-4*h)
return (f0+f1+f2+f3+f4)/h
def diff2(func, x, h):
f0 = (35/12)*func(x)
f1 = -(26/3)*func(x-h)
f2 = (19/2)*func(x-2*h)
f3 = -(14/3)*func(x-3*h)
f4 = (11/12)*func(x-4*h)
return (f0+f1+f2+f3+f4)/h/h
def diff3(func, x, h):
f0 = (5/2)*func(x)
f1 = -9*func(x-h)
f2 = 12*func(x-2*h)
f3 = -7*func(x-3*h)
f4 = (3/2)*func(x-4*h)
return (f0+f1+f2+f3+f4)/h/h/h
L = 1000
h = 2*np.pi/L
cos_arr = np.cos(np.linspace(0,2*np.pi,L))
sin_arr = np.sin(np.linspace(0,2*np.pi,L))
diff_arr = np.array([ diff(np.cos,x,h) for x in np.linspace(0,2*np.pi,L) ])
diff2_arr = np.array([ diff2(np.cos,x,h) for x in np.linspace(0,2*np.pi,L) ])
diff3_arr = np.array([ diff3(np.cos,x,h) for x in np.linspace(0,2*np.pi,L) ])
plt.plot(cos_arr)
plt.plot(sin_arr)
plt.plot(-cos_arr)
plt.plot(-sin_arr)
plt.plot(diff_arr,lw=0,marker="v",ms=0.5)
plt.plot(diff2_arr,lw=0,marker="x",ms=0.5)
plt.plot(diff3_arr,lw=0,marker="+",ms=0.5)
print(max(abs(diff_arr+sin_arr)),max(abs(diff2_arr+cos_arr)),max(abs(diff3_arr-sin_arr)))