数値微分¶

In [1]:
import numpy as np
from scipy.misc import derivative # 数値微分
import matplotlib.pyplot as plt
import pandas as pd
In [2]:
import scienceplots
plt.style.use(['science', 'notebook'])
plt.rcParams['font.family'] = 'Times New Roman' # font familyの設定
plt.rcParams['font.family'] = 'serif' # font familyの設定
plt.rcParams['mathtext.fontset'] = 'cm' # math fontの設定

中心差分法¶

関数$f(x)$の1階微分を次のように近似する． \begin{gather} f'(x) \simeq \frac{f(x+h)-f(x-h)}{2h} \end{gather} 精度は $O(h^2)$．前進差分や後退差分に比べて誤差が小さい．

In [3]:
def my_diff1(func, x, h): # 通常の中心差分公式（Central Difference Method）
    scalar_input = np.isscalar(x)   # ← スカラーかどうか判定しておく！
    x = np.asarray(x) # 引数xをNumPy配列（ndarray）に変換．スカラーでも配列でも対応
    fp = func(x+h)
    fm = func(x-h)
    d1 = (fp-fm)/2.0/h
    return d1.item() if scalar_input else d1

リチャードソン外挿¶

4次精度の中心差分式で，関数$f(x)$の1階微分の近式式は次のようになる． \begin{gather} f'(x) \simeq \frac{-f(x+2h)+8f(x+h)-8f(x-h)+f(x-2h)}{12h} \end{gather} 精度は $O(h^4)$．

In [4]:
def my_diff2(func, x, h): # リチャードソンの外挿の式（Richardson extrapolation）
    scalar_input = np.isscalar(x)   # ← スカラーかどうか判定しておく！
    x = np.asarray(x) # 引数xをNumPy配列（ndarray）に変換．スカラーでも配列でも対応
    fp = func(x+h)
    fm = func(x-h)
    fpp = func(x+h+h)
    fmm = func(x-h-h)
    d2 = (-fpp+8.0*(fp-fm)+fmm)/12.0/h
    return d2.item() if scalar_input else d2

6次精度

In [5]:
def seven_point_central_difference(f, x, h=1e-5): # 6次精度
    scalar_input = np.isscalar(x)
    x = np.asarray(x)
    d = (-f(x - 3*h) + 9*f(x - 2*h) - 45*f(x - h)
         + 45*f(x + h) - 9*f(x + 2*h) + f(x + 3*h)) / (60*h)
    return d.item() if scalar_input else d
In [6]:
def func1(x):
    y = 1.0 + x + np.sin(x)
    return y
In [7]:
nx = 6
x = np.linspace(0.0, 1.0, nx)
h = 0.01
df1 = np.array([my_diff1(func1, xi, h) for xi in x]) # 中心差分公式
df2 = np.array([my_diff2(func1, xi, h) for xi in x]) # リチャードソンの外挿の式
df3 = np.array([seven_point_central_difference(func1, xi, h) for xi in x]) # 高次中央差分法（7点）による1階微分
In [8]:
print("f1(x) = 1+x+sin(x)")
print(f"{'x':>6s} {'diff1':>10s} {'error1':>10s}"\
                f" {'diff2':>10s}{'error2':>15s}"\
                f" {'higher':>10s}{'error3':>18s}")
dydx = 1.0+np.cos(x)
for i in range(nx):
    print(f"{x[i]:>6.3f} {df1[i]:>10.7f} {df1[i]-dydx[i]:>10.7f}"\
                        f" {df2[i]:>10.7f}{df2[i]-dydx[i]:>15.11f}"\
                        f" {df3[i]:>10.7f}{df3[i]-dydx[i]:>18.14f}")
print("")

f1(x) = 1+x+sin(x)
     x      diff1     error1      diff2         error2     higher            error3
 0.000  1.9999833 -0.0000167  2.0000000 -0.00000000033  2.0000000 -0.00000000000002
 0.200  1.9800502 -0.0000163  1.9800666 -0.00000000033  1.9800666 -0.00000000000002
 0.400  1.9210456 -0.0000154  1.9210610 -0.00000000031  1.9210610 -0.00000000000002
 0.600  1.8253219 -0.0000138  1.8253356 -0.00000000028  1.8253356  0.00000000000001
 0.800  1.6966951 -0.0000116  1.6967067 -0.00000000023  1.6967067  0.00000000000000
 1.000  1.5402933 -0.0000090  1.5403023 -0.00000000018  1.5403023 -0.00000000000002
In [9]:
df = pd.DataFrame({
    '$x$':x,
    'diff1':df1,
    'error1':df1-dydx,
    'diff2':df2,
    'error2':df2-dydx,
    'higher':df3,
    'error3':df3-dydx,
})
df.style.format({"$x$":"{:.1f}", "error1":"{:.6f}", "error2":"{:.11f}", "error3":"{:.14f}"})
Out[9]:
  $x$ diff1 error1 diff2 error2 higher error3
0 0.0 1.999983 -0.000017 2.000000 -0.00000000033 2.000000 -0.00000000000002
1 0.2 1.980050 -0.000016 1.980067 -0.00000000033 1.980067 -0.00000000000002
2 0.4 1.921046 -0.000015 1.921061 -0.00000000031 1.921061 -0.00000000000002
3 0.6 1.825322 -0.000014 1.825336 -0.00000000028 1.825336 0.00000000000001
4 0.8 1.696695 -0.000012 1.696707 -0.00000000023 1.696707 0.00000000000000
5 1.0 1.540293 -0.000009 1.540302 -0.00000000018 1.540302 -0.00000000000002
In [10]:
# df.to_latex("table.tex", index=False, escape=False)
In [11]:
nxx = 101
xx = np.linspace(0.0, 10, nxx)
yy = 1.0+np.cos(xx)
nd = 20
# 各微分法での近似値を個別にリスト内包で計算
d1_vals = [my_diff1(func1, xi, h) for xi in xx[::nd]]
d2_vals = [my_diff2(func1, xi, h) for xi in xx[::nd+1]]
d3_vals = [seven_point_central_difference(func1, xi, h) for xi in xx[::nd+2]]
fig = plt.figure()
plt.plot(xx, func1(xx), label='$f(x) = 1+x+\sin(x)$')
plt.plot(xx, 1.0 + np.cos(xx), label=r"$f'(x)=1+\cos(x)$")
plt.plot(xx[::nd], d1_vals, 'o', label='Central difference')
plt.plot(xx[::nd+1], d2_vals, 'v', label='Richardson extrapolation')
plt.plot(xx[::nd+2], d3_vals, 's', label='7-point central diff')
plt.legend(ncol=1, loc='upper right', fancybox=False, frameon=True)
plt.xlabel('$x$')
plt.ylabel("$f'(x)$")
plt.title('Numerical Differentiation Methods')
plt.grid(True)
fig.savefig('p2_ex02_1.pdf')
plt.show()
No description has been provided for this image