数値微分

モジュールのインポート

import numpy as np
from scipy.misc import derivative # 数値微分
import matplotlib.pyplot as plt

import scienceplots
#Warning : As of version 2.0.0, you need to add import scienceplots before setting the style (plt.style.use('science')).
plt.style.use(['science', 'notebook'])

中心差分公式

def my_diff1(func, x, h): # 通常の中心差分公式
    fp = func(x+h)
    fm = func(x-h)
    d1 = (fp-fm)/2.0/h
    return d1

リチャードソンの外挿

def my_diff2(func, x, h): # リチャードソンの外挿の式
    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

微分する関数

def func1(x):
    y = 1.0 + x + np.sin(x)
    return y

nx = 6
dx = 0.2
h = 0.01
print("f1(x) = 1+x+sin(x)")
print(f"{'x':>8s} {'diff1':>10s} {'error1':>10s}"\
      f"{'diff2':>10s} {'error2':>10s} {'derivative':>10s} {'error3':>10s}")
for i in range(nx):
    x = i*dx
    df1 = my_diff1(func1, x, h)# 中心差分公式
    df2 = my_diff2(func1, x ,h)# リチャードソンの外挿の式
    df3 = derivative(func1, x, h)# 中心差分公式（scipyの関数）
    dydx = 1.0+np.cos(x)
    print(f"{x:>8.3f} {df1:>10.7f} {df1-dydx:>10.7f}"\
          f"{df2:>10.7f} {df2-dydx:>10.7f} {df3:>10.7f} {df3-dydx:>10.7f}")
print("")

f1(x) = 1+x+sin(x)
       x      diff1     error1     diff2     error2 derivative     error3
   0.000  1.9999833 -0.0000167 2.0000000 -0.0000000  1.9999833 -0.0000167
   0.200  1.9800502 -0.0000163 1.9800666 -0.0000000  1.9800502 -0.0000163
   0.400  1.9210456 -0.0000154 1.9210610 -0.0000000  1.9210456 -0.0000154
   0.600  1.8253219 -0.0000138 1.8253356 -0.0000000  1.8253219 -0.0000138
   0.800  1.6966951 -0.0000116 1.6967067 -0.0000000  1.6966951 -0.0000116
   1.000  1.5402933 -0.0000090 1.5403023 -0.0000000  1.5402933 -0.0000090

nx = 7
xp = np.linspace(0.0, 10, nx)
yy1 = np.empty(nx)
yy1[:] = [my_diff1(func1, xp[i], h) for i in range(nx)]
print(f"{'diff1:':>11s}",end=' ')
[print(f"{yy1[i]:10.7f}",end=' ') for i in range(nx)]
print("")

diff1:  1.9999833  0.9042780  0.0183424  1.2836575  1.9273522  0.5388036  0.1609425

xp2 = np.linspace(0.2, 10, nx)
yy2 = np.empty(nx)
yy2[:] = [my_diff2(func1, xp2[i], h) for i in range(nx)]
print(f"{'diff2:':>11s}",end=' ')
[print(f"{yy2[i]:10.7f}",end=' ') for i in range(nx)]
print("")

diff2:  1.9800666  0.7404685  0.0523729  1.3779777  1.9003827  0.5094812  0.1609285

xp3 = np.linspace(0.4, 10, nx)
yy3 = np.empty(nx)
yy3[:] = [derivative(func1, xp3[i], h) for i in range(nx)]
print(f"{'derivative:':>11s}",end=' ')
[print(f"{yy3[i]:10.7f}",end=' ') for i in range(nx)]
print("")

derivative:  1.9210456  0.5838601  0.1032565  1.4685089  1.8693830  0.4807200  0.1609425

nxx = 101
xx = np.linspace(0.0, 10, nxx)
yy = 1.0+np.cos(xx)

fig = plt.figure() # グラフ領域の作成
plt.plot(xx,func1(xx))
plt.plot(xx,yy)
plt.plot(xp,yy1,'o')
plt.plot(xp2,yy2,'v')
plt.plot(xp3,yy3,'s')
fig.savefig('p2_ex02_1.pdf')