n階の数値微分

# テーラー展開

テーラー展開可能な関数$f$を考えます。

```math
\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$は整数）

```math
f(x-mh) = \sum_{n=0}^{\infty} \frac{(-m)^n}{n!}f^{(n)}(x) h^n
```

これを$m=1,\cdots M$について$n \leq M$の範囲で展開し、これを行列形式で並べます。

```math
\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$の各成分は

```math
(A_M)_{ab} = \begin{cases}
    \frac{(-a)^b}{b!} \quad (a>0) \\
    1 \quad(a=b=0) \\
    0 \quad (b>a=0)
\end{cases}

```

です。

このとき、逆行列を使って

```math
\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$は次のように表せます。

```math
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次のオーダーまでテーラー展開してみます。

```math
\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}
```

よって

```math
\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

```math
\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}
```

なので、

```math
\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

```math
\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}
```

だから、

```math
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$)

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2445529/190949bf-fa48-f313-a973-d0ca13cfedea.png)

また、解析解との誤差の最大値はそれぞれ以下の通りです。

```math
\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)))
```