Abstract
みなさん微分していますか?私はしています。 数値微分は自動微分のテストに使えるし実装が楽で嬉しいです。 なので精度よくやりたくなってきます。 ところで前進差分の打ち切り誤差が有限の幅
調べたのですが出てこなかったので導出してみます。
前進差分・中心差分の精度の確認
まずは前進差分と中心差分の精度を確認してみます。
前進差分は以下のような計算で得られます。
これが一次精度、つまり誤差
上の式の右辺をテイラー展開すると、
中心差分についても、全く同様に右辺をテイラー展開することで確かめられます。
Pythonでそれぞれ実装してみます。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
matplotlib.use("module://matplotlib_pyodide.html5_canvas_backend")
def numerical_diff(f, x, h):
return (f(x+h) - f(x)) / h
def numerical_diff_center(f, x, h):
return (f(x+h) - f(x-h)) / (2*h)
def f(x):
return np.sin(x)
h = np.logspace(0, -20, num=100)
grad_forward = numerical_diff(f, np.pi/3, h)
grad_center = numerical_diff_center(f, np.pi/3, h)
grad_true = np.cos(np.pi/3)
error_forward = np.abs(grad_forward - grad_true)
error_center = np.abs(grad_center - grad_true)
plt.plot(h, error_forward, label='forward')
plt.plot(h, error_center, label='center')
plt.xscale('log')
plt.yscale('log')
plt.xlabel('h')
plt.ylabel('grad')
plt.legend()
plt.show()
plt.close()
よく言われる、打ち切り誤差と桁落ちのトレードオフが確認できます
任意の k 階 n 次精度数値微分
本題を導出しようと思います。 まずは議論のしやすさのために特に
任意の
は
つまり、 各
証明してみます。
これを
したがって、
と
ここで、
したがって、
正則性について
ヴァンデルモンド行列は、という各列が等比数列のようになっている行列のことを言います。 (各行の派閥もいるらしいです。)
ヴァンデルモンド行列は、その行列式が
となることが知られていて、
したがって今回の場合は
というわけで確かに存在したうえに一意で、また具体的に計算することができました。
実際に計算してみようと思います。
右から
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
import numpy as np
import matplotlib.pyplot as plt
def central_weight(n, k):
p = n // 2
A = np.zeros((n+1, n+1), dtype=float)
r = np.arange(-p, p+1)
for i in range(n+1):
A[i] = r**i
w = np.linalg.inv(A)[:, k]
return w
def numerical_diff(f, x, n, h=1e-5, k=1):
w = central_weight(n, k)
p = n // 2
w = central_weight(n, k)
x = np.arange(-p, p+1) * h + x
y = f(x)
return np.dot(w, y) * np.prod(np.arange(1, k+1)) / h**k
def f(x):
return np.sin(x)
def df(x):
return np.cos(x)
x = np.pi / 3
print('w_2,1 =', central_weight(2, 1))
print('w_4,1 =', central_weight(4, 1))
print('w_4,2 =', central_weight(4, 2))
print(numerical_diff(f, x, 2))
print(numerical_diff(f, x, 4))
print(numerical_diff(f, x, 4, k=2))
と、いい感じに計算できています。
桁落ちとのトレードオフも確認してみます。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
n = range(2, 12, 2)
x = np.pi / 3
h = np.logspace(0, -10, 100)
label = [f'O(h^{i})' for i in n]
err = np.zeros((len(n), len(h)))
for i, _n in enumerate(n):
for j, _h in enumerate(h):
err[i, j] = abs(numerical_diff(f, x, _n, h=_h) - df(x))
plt.figure(figsize=(8, 6))
for i in range(len(n)):
plt.loglog(h, err[i], label=label[i])
plt.legend()
plt.grid()
plt.xlabel('h')
plt.ylabel('Error')
plt.show()
打ち切り誤差の変化によって最適な点が変化していることがわかります。
感想
線形代数が役に立ってすごいまた、上のアルゴリズムは逆行列を求めるところがボトルネックで 時間計算量は