数値微分

モジュールのインポート

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 func1(x):
    y = 1.0 + x + np.sin(x)
    return y

関数 derivative

nx, dx, h = 6, 0.2, 0.01
print("f1(x) = 1+x+sin(x)")
print(f"{'x':>8s} {'derivative':>10s} {'error':>10s}")
for i in range(nx):
    x = i*dx
    df3 = derivative(func1, x, h) # 中心差分公式（scipyの関数）
    dydx = 1.0+np.cos(x)
    print(f"{x:>8.3f} {df3:>10.7f} {df3-dydx:>10.7f}")

f1(x) = 1+x+sin(x)
       x derivative      error
   0.000  1.9999833 -0.0000167
   0.200  1.9800502 -0.0000163
   0.400  1.9210456 -0.0000154
   0.600  1.8253219 -0.0000138
   0.800  1.6966951 -0.0000116
   1.000  1.5402933 -0.0000090

プロットデータの計算

nxx = 101
xx = np.linspace(0.0, 10, nxx)
xx

array([ 0. ,  0.1,  0.2,  0.3,  0.4,  0.5,  0.6,  0.7,  0.8,  0.9,  1. ,
        1.1,  1.2,  1.3,  1.4,  1.5,  1.6,  1.7,  1.8,  1.9,  2. ,  2.1,
        2.2,  2.3,  2.4,  2.5,  2.6,  2.7,  2.8,  2.9,  3. ,  3.1,  3.2,
        3.3,  3.4,  3.5,  3.6,  3.7,  3.8,  3.9,  4. ,  4.1,  4.2,  4.3,
        4.4,  4.5,  4.6,  4.7,  4.8,  4.9,  5. ,  5.1,  5.2,  5.3,  5.4,
        5.5,  5.6,  5.7,  5.8,  5.9,  6. ,  6.1,  6.2,  6.3,  6.4,  6.5,
        6.6,  6.7,  6.8,  6.9,  7. ,  7.1,  7.2,  7.3,  7.4,  7.5,  7.6,
        7.7,  7.8,  7.9,  8. ,  8.1,  8.2,  8.3,  8.4,  8.5,  8.6,  8.7,
        8.8,  8.9,  9. ,  9.1,  9.2,  9.3,  9.4,  9.5,  9.6,  9.7,  9.8,
        9.9, 10. ])

yy = 1.0+np.cos(xx)
yy

array([2.00000000e+00, 1.99500417e+00, 1.98006658e+00, 1.95533649e+00,
       1.92106099e+00, 1.87758256e+00, 1.82533561e+00, 1.76484219e+00,
       1.69670671e+00, 1.62160997e+00, 1.54030231e+00, 1.45359612e+00,
       1.36235775e+00, 1.26749883e+00, 1.16996714e+00, 1.07073720e+00,
       9.70800478e-01, 8.71155506e-01, 7.72797905e-01, 6.76710433e-01,
       5.83853163e-01, 4.95153895e-01, 4.11498883e-01, 3.33723979e-01,
       2.62606284e-01, 1.98856384e-01, 1.43111247e-01, 9.59278580e-02,
       5.77776593e-02, 2.90418349e-02, 1.00075034e-02, 8.64849727e-04,
       1.70522421e-03, 1.25202301e-02, 3.32018074e-02, 6.35433127e-02,
       1.03241584e-01, 1.51899968e-01, 2.09032288e-01, 2.74067696e-01,
       3.46356379e-01, 4.25176053e-01, 5.09739179e-01, 5.99200828e-01,
       6.92667130e-01, 7.89204201e-01, 8.87847473e-01, 9.87611337e-01,
       1.08749898e+00, 1.18651237e+00, 1.28366219e+00, 1.37797774e+00,
       1.46851667e+00, 1.55437434e+00, 1.63469288e+00, 1.70866977e+00,
       1.77556588e+00, 1.83471278e+00, 1.88551952e+00, 1.92747843e+00,
       1.96017029e+00, 1.98326844e+00, 1.99654210e+00, 1.99985864e+00,
       1.99318492e+00, 1.97658763e+00, 1.95023259e+00, 1.91438315e+00,
       1.86939749e+00, 1.81572510e+00, 1.75390225e+00, 1.68454667e+00,
       1.60835131e+00, 1.52607752e+00, 1.43854733e+00, 1.34663532e+00,
       1.25125984e+00, 1.15337386e+00, 1.05395542e+00, 9.53997874e-01,
       8.54499966e-01, 7.56455846e-01, 6.60845139e-01, 5.68623155e-01,
       4.80711346e-01, 3.97988097e-01, 3.21279953e-01, 2.51353354e-01,
       1.88906986e-01, 1.34564791e-01, 8.88697381e-02, 5.22783979e-02,
       2.51563786e-02, 7.77467455e-03, 3.06957965e-04, 2.82784380e-03,
       1.53121442e-02, 3.76351202e-02, 6.95737279e-02, 1.10808847e-01,
       1.60928471e-01])

nx = 6
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.3189555  0.5449927  1.9924057  0.7759048  0.1609425 

プロット

fig = plt.figure() # グラフ領域の作成
plt.plot(xx,func1(xx))
plt.plot(xx,yy)
plt.plot(xp3,yy3,'s')
fig.savefig('p1_ex02_1.pdf')

別の関数で計算

def func2(x):
    y =  x * np.sin(x)
    return y

nx = 6
h = 0.01
xp3 = np.linspace(0.4, 10, nx)
yy3 = np.empty(nx)
yy3[:] = [derivative(func2, 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:  0.7578171 -0.8478285 -2.8196704  5.9903515 -0.8361716 -8.9345694 

fig = plt.figure() # グラフ領域の作成
plt.plot(xx,func2(xx))
plt.plot(xx,np.sin(xx)+xx*np.cos(xx))
plt.plot(xp3,yy3,'s')
fig.savefig('p1_ex02_2.pdf')