\usepackagecancel

Return to Python で数値解析

SciPy で数値解析

Python で数値解析するためのライブラリである SciPy の使用例を示します。グラフは Matplotlibpyplot.plot() を使って描きます。三角関数などの数学関数については,SymPy を使わず,全て NumPy の関数を使います。

ライブラリの import と設定

In [1]:
# NumPy を使います
import numpy as np
# Matplotlib でグラフを描きます
import matplotlib.pyplot as plt

# グラフの目盛設定用
from matplotlib.ticker import MultipleLocator

# 以下はグラフを SVG で Notebook にインライン表示させる設定
%config InlineBackend.figure_formats = ['svg']

# mathtext font の設定
plt.rcParams['mathtext.fontset'] = 'cm'

SciPy による数値微分

SciPy には数値微分をしてくれる関数 scipy.misc.derivative がありましたが,非推奨となり,SciPy 1.12.0 からは使えなくなったようです。

解析的な微分の定義は(前方差分の極限として)

dfdxlimh0f(x+h)f(x)h

でした。滑らかな関数であれば,(後方差分の極限として)
dfdx=limh0f(x)f(xh)h
としても良いです。

現実世界では h0 の極限はとれませんから十分小さい値として h を定義し,近似的な数値微分を(前方差分と後方差分の平均である中心差分として)以下のように定義しましょう。

dfdx12{f(x+h)f(x)h+f(x)f(xh)h}=f(x+h)f(xh)2h

関数 f(x) が初等関数(および初等関数の合成関数)である場合は解析的に微分できますが,そうでない場合,例えば f(x) が数値的にしか与えられていないとか,何か解析的に微分できない場合には,数値微分によって近似的な値を求める必要があります。

また,解析的に微分できそうだが,人力ではけっこう大変そうな場合にも,数値微分によって近似的な値を使うほうが便利な場合があるかもしれません。

以下のように,関数 f(x) を数値微分してくれる関数 ndiff() を定義します。

In [2]:
def ndiff(f, x, h=1e-5):
    # h を省略すると h=1e-5 として計算する
    return (f(x + h) - f(x - h))/(2*h)

例として,f(x)=sinxx=3 における数値微分を求めてみます。(答えは cos3 です。)

In [3]:
# h を省略した場合
print('       ndiff(f, 3) = %.15f' % ndiff(np.sin, 3))

for h in [1e-4, 1e-5, 1e-6, 1e-7, 1e-8]:
    print('h = %3.1e, ndiff = %.15f' % (h, ndiff(np.sin, 3, h)))

print(' '*10 + 'cos(%3.1f) = %.15f' % (3, np.cos(3)))
       ndiff(f, 3) = -0.989992496590319
h = 1.0e-04, ndiff = -0.989992494952602
h = 1.0e-05, ndiff = -0.989992496590319
h = 1.0e-06, ndiff = -0.989992496730485
h = 1.0e-07, ndiff = -0.989992494926373
h = 1.0e-08, ndiff = -0.989992490763036
          cos(3.0) = -0.989992496600445

以上の結果からわかるように,h をやみくもに小さくすればするほど,数値微分の精度があがる,というわけではないです。

○練習:関数の数値微分

f(x)=x3ex1 の極大値を求めるための準備として,まず,f(x) とその数値微分 ndiff(f, x)0x5 の範囲でグラフにする。

In [4]:
def f(x):
    return x**3/(np.exp(x) - 1)
In [5]:
# ax を使う際の最初のおまじない
fig, ax = plt.subplots()

# x の範囲
xrange = np.linspace(0.001, 5, 200)

# y = f(x) のグラフ
ax.plot(xrange, f(xrange), label = '$f(x)$')
# y = f'(x) のグラフ
ax.plot(xrange, ndiff(f, xrange), label = '$f^{\prime}(x)$')

# x 軸のラベル
ax.set_xlabel('$x$')
# グリッド
ax.grid()
# 凡例の表示
ax.legend();

上のグラフから,数値微分した微分係数 f(x) の値がゼロになる x は,2<x<3 のあたりであること,そのあたりで確かに f(x) が極大(最大)となっていることが見てとれます。

SciPy による数値積分

scipy.integrate.quad()

SciPy には数値積分してくれる関数 scipy.integrate.quad() が用意されています。それを使ってみましょう。結果には積分値と推定誤差が出力されます。

0πsin(sinx)dx の計算をします。

数値計算の際は,三角関数等は NumPy の関数を使い,円周率 π も NumPy の np.pi を使います。

In [6]:
from scipy.integrate import quad

# quad を使うときは変数名 x 決め打ち
def sinsin(x):
    return np.sin(np.sin(x))

# quad で数値積分
# quad(f, a, b): 被積分関数 f(x) を x = a から b まで数値積分
quad(sinsin, 0, np.pi)
Out[6]:
(1.7864874819500525, 2.4687274069628184e-11)

scipy.integrate.quad() は上記のように,数値積分の値と絶対誤差の推定値をタプル (tuple) で返します。

数値積分値(と誤差)を個別に取り出すには,例えば以下のように。

In [7]:
y, abserr = quad(sinsin, 0, np.pi)

print("積分値は", y)
print("誤差は", abserr)
積分値は 1.7864874819500525
誤差は 2.4687274069628184e-11

参考:シンプソン法による数値積分

シンプソン法で abf(x)dx を求める。

積分区間 [a,b]N(=2n) 等分(偶数等分)し,

h=baN,fi=f(a+ih)

とする。シンプソン法の公式は

abf(x)dxh3(f0+4f1+f2)+h3(f2+4f3+f4)+h3(f4+4f5+f6)+h3(f6+4f7+f8)++h3(fN2+4fN1+fN)=h3(f0+2i=1n1f2i+4i=1nf2i1+f2n)

SciPy では scipy.integrate.quad() を使って数値積分することができるのでわざわざ自作の関数をつくる必要はないのですが,教育的見地からあえて,シンプソン法の公式を使って数値積分の値を返す関数 simpson() を以下のように定義してみます。

In [8]:
def simpson(f, a, b, Ndiv):
    h = (b - a)/Ndiv
    simp = 0
    for i in range(0, Ndiv, 2):
        simp += f(a + i*h) + 4*f(a + (i+1)*h) + f(a + (i+2)*h)
    return simp*h/3

N の値を変えて 0πsin(sinx)dxsimpson() で数値積分して,精度を確認してみます。

In [9]:
def f(x):
    return np.sin(np.sin(x))

for N in [10*2**i for i in range(9)]:
    print('N = %4d, %.15f' %(N, simpson(f, 0, np.pi, N)))
N =   10, 1.786721213744079
N =   20, 1.786501256275661
N =   40, 1.786488331266680
N =   80, 1.786487534856180
N =  160, 1.786487485253951
N =  320, 1.786487482156504
N =  640, 1.786487481962955
N = 1280, 1.786487481950862
N = 2560, 1.786487481950104

以上の結果から,どんなことがわかるかな?

SciPy による方程式の数値解

例としてベッセル関数 J0(x) のゼロ点を求めてみます。SciPy におけるベッセル関数 J0(x)scipy.special.jv(0,x) です。

まずは J0(x)=0 となる x のあたりをつけるためにグラフを描いてみます。

In [10]:
from scipy.special import jv

def J0(x):
    return jv(0,x)

# ax を使う際の最初のおまじない
fig, ax = plt.subplots()

xrange = np.linspace(0, 10, 200)
ax.plot(xrange, J0(xrange))

# x 軸
ax.plot(xrange, 0*xrange)

ax.set_xlabel('x')
ax.set_ylabel('J0(x)')
ax.grid()
ax.set_xlim(0, 10)
ax.set_title('ベッセル関数 J0(x)');

scipy.optimize.root_scalar

方程式の数値解を求めてくれる関数は scipy.optimize.root_scalar です。2<x<3 のあたりで x 軸と交わっていますから,このあたりの解を探します。

In [11]:
from scipy.optimize import root_scalar

# x が 2 ~ 3 の間で j0(x) = 0 の根を探す
root_scalar(J0, bracket = [2, 3])
Out[11]:
      converged: True
           flag: 'converged'
 function_calls: 8
     iterations: 7
           root: 2.404825557695807

scipy.optimize.root_scalar() は結果を上記のように表示します。解 root のみを表示したい場合は以下のようにします。

In [12]:
# 根のみを出力するには...
ans = root_scalar(J0, bracket = [2, 3])
print(ans.root)
2.404825557695807
In [13]:
# ほんとに根になってるか確認。ゼロに近ければよい。
J0(ans.root)
Out[13]:
-1.7918655681235457e-14
In [14]:
# どうしてももっと精度が欲しいというのであれば...
# xtol に,小さい値をいれてみる
ans = root_scalar(J0, bracket = [2, 3], xtol=1e-15)
ans.root
Out[14]:
2.4048255576957724
In [15]:
# ほんとに根になってるか確認。ゼロに近ければよい。
J0(_)
Out[15]:
1.535884772676773e-16

○練習:関数の極大値

  1. f(x)=x3ex1 が極大値をとるときの x の値,およびそのときの極大値を求める。
  2. g(x)=x5ex1 が極大値をとるときの x の値,およびそのときの極大値も同様にして求めてみよ。

まず f(x) の数値微分 ndiff(f, x) をグラフにしてみる。

In [16]:
def f(x):
    return x**3/(np.exp(x) - 1)

# ax を使う際の最初のおまじない
fig, ax = plt.subplots()

xrange = np.linspace(0, 5, 200)
ax.plot(xrange, ndiff(f, xrange), label = '$f^{\prime}(x)$')

# x 軸
ax.plot(xrange, 0*xrange)

ax.set_xlabel('$x$')
ax.set_ylabel(' ')
ax.grid()
ax.legend();

ndiff(f, x)2<x<3 のあたりでゼロになっているので,そのあたりを root_scalar() で探してみる。

In [17]:
def df(x):
    return ndiff(f, x)

ans = root_scalar(df, bracket = [2, 3])
ans.root
Out[17]:
2.8214393721317643
In [18]:
# そのときの f(x) の極大値は...
f(_)
Out[18]:
1.4214354727477367

SciPy による1階常微分方程式の数値解法

scipy.integrate.solve_ivp()

SciPy には常微分方程式を数値的に解いてくれる関数 scipy.integrate.solve_ivp() が用意されています。それを使ってみましょう。

例として,簡単なロジスティック方程式

dxdt=f(x)=(1x)x

を,初期条件 t0=0x0=x(t0)=0.1t=t0=0 から t1=10 まで解きます。

solve_ivp() を使う場合は変数名が決め打ちで t, y です。連立常微分方程式の場合は y がリストになります。

In [19]:
from scipy.integrate import solve_ivp

# scipy.integrate.solve_ivp() は
# dy/dt = f(t, y) の形を解く
# solve_ivp を使う際は変数名 t, y は決め打ち
# solve.ivp() に渡す関数の引数の順番に注意。t が先
# 連立でないときは,1次元リストとして返す

def f(t, y):
    x = y[0]
    return [(1-x)*x]

# t の範囲
tini = 0
tend = 10
# 計算する点は t の範囲を N 等分して N+1 点
N = 20
t_list = np.linspace(tini, tend, N+1)
# y の初期値
y_ini = [0.1]

answer = solve_ivp(f, [tini, tend], y_ini, t_eval = t_list, 
                   # デフォルトの許容誤差は大きめなので
                   # 小さい値に設定
                   rtol = 1.e-12, 
                   atol = 1.e-15)
print(answer)
  message: The solver successfully reached the end of the integration interval.
  success: True
   status: 0
        t: [ 0.000e+00  5.000e-01 ...  9.500e+00  1.000e+01]
        y: [[ 1.000e-01  1.548e-01 ...  9.993e-01  9.996e-01]]
      sol: None
 t_events: None
 y_events: None
     nfev: 1790
     njev: 0
      nlu: 0
In [20]:
# ax を使う際の最初のおまじない
fig, ax = plt.subplots()

ax.scatter(answer.t, answer.y[0], s = 10)

ax.grid()
ax.set_xlabel("$t$")
ax.set_ylabel("$x(t)$")
ax.set_title("ロジスティック方程式の数値解");

参考までに,解析解は

x(t)=exp(t)9+exp(t)

です。(解析解がわかっているときには,わざわざ数値計算する必要はないのですが,練習ということで。)

In [21]:
# 解析解
def xsol(t):
    return np.exp(t)/(9+np.exp(t))

# ax を使う際の最初のおまじない
fig, ax = plt.subplots()

# 解析解
ts = np.linspace(tini, tend, 200)
ax.plot(ts, xsol(ts), c = 'tab:blue', label='解析解')
# 数値解
ax.scatter(answer.t, answer.y[0], zorder = 2,
           c = 'tab:orange', s = 10, label = '数値解')

ax.grid()
ax.set_xlabel("$t$")
ax.set_ylabel("$x(t)$")
ax.set_title("ロジスティック方程式の解")
ax.legend();

SciPy による2階常微分方程式の数値解法

scipy.integrate.solve_ivp()

SciPy には常微分方程式を数値的に解いてくれる関数 scipy.integrate.solve_ivp() が用意されています。この函数は,1階常微分方程式だけでなく,連立1階常微分方程式も,従って高階の常微分方程式も解いてくれるので,それを使ってみましょう。

例として,簡単な減衰+強制振動の方程式

d2xdt2=xadxdt+bcos(t)

を,

  • a,b にいくつか値を入れて
  • 初期条件 t0=0x0=x(t0)=3,v0=v(t0)=0
  • t=t0=0 から t1=20 までを解きます。

2階常微分方程式は,

dxdt=F1(t,x,v)=vdvdt=F2(t,x,v)=xav+bcost

として連立1階微分方程式の形にして解きます。

In [22]:
from scipy.integrate import solve_ivp

# scipy.integrate.solve_ivp() は
# dy/dt = F(t, y) の形を解く

def F(t, y, a, b):
# solve.ivp() に渡す関数の引数の順番に注意。t が先
    x = y[0]
    v = y[1]
    F1 = v
    F2 = -x - a*v + b*np.cos(t)
    return [F1, F2]

tini = 0
tend = 20
t_list = np.linspace(tini, tend, 200)
y_ini = [3, 0]

a = 0.5; b = 0.2
answer1 = solve_ivp(F, [tini, tend], y_ini, t_eval = t_list, 
                    args=(a, b), 
                    # デフォルトの許容誤差は大きめなので
                    # 小さい値に設定
                    rtol = 1.e-12, 
                    atol = 1.e-15)

print(answer1)
  message: The solver successfully reached the end of the integration interval.
  success: True
   status: 0
        t: [ 0.000e+00  1.005e-01 ...  1.990e+01  2.000e+01]
        y: [[ 3.000e+00  2.986e+00 ...  3.669e-01  3.840e-01]
            [ 0.000e+00 -2.740e-01 ...  1.881e-01  1.509e-01]]
      sol: None
 t_events: None
 y_events: None
     nfev: 8834
     njev: 0
      nlu: 0

横軸に t= answer1.t,縦軸に x(t)= answer1.y[0] をとってグラフにしてみます。

In [23]:
# ax を使う際の最初のおまじない
fig, ax = plt.subplots()

# 数値解 x は .y[0] のようにして取り出す
# 数式を label に設定する際は $ で囲んで LaTeX 式に
# a, b の値を凡例に
key = '$a=%3.1f, b=%3.1f$' % (a, b)
ax.plot(answer1.t, answer1.y[0], label=key)

ax.grid()
ax.set_xlabel("$t$")
ax.set_ylabel("$x(t)$")
ax.legend()
ax.set_title("減衰+強制振動")
ax.set_xlim(0, 20);

○練習:減衰+強制振動

減衰+強制振動の方程式を以下の場合に解き,グラフにする。

まず a の影響を調べるため,

  1. a=0,b=0
  2. a=0.2,b=0
  3. a=0.5,b=0
  4. a=0.8,b=0

次に,b の影響を調べるため,

  1. a=0,b=0
  2. a=0,b=0.2
  3. a=0,b=0.5
  4. a=0,b=0.8

○練習:単振り子

2階常微分方程式

d2θdt2=4π2sinθ

を初期条件

θ=45,  dθdt=0at  t=0

のもとで t=0 から t=5 まで数値的に解け。

θ y[0]θ˙ y[1] として

In [24]:
def F(t, y):
    # solve.ivp() に渡す関数の引数の順番に注意。t が先
    theta = y[0]
    dottheta = y[1]
    F1 = dottheta
    F2 = -4*np.pi**2*np.sin(theta)
    return [F1, F2]

tini = 0
tend = 5
t_list = np.linspace(tini, tend, 200)
# 初期条件はラジアンになおして
y_ini = [np.radians(45), 0]

answer2 = solve_ivp(F, [tini, tend], y_ini, t_eval = t_list, 
                    # デフォルトの許容誤差は大きめなので
                    # 小さい値に設定
                    rtol = 1.e-12, 
                    atol = 1.e-15)

一方,θ があまり大きくないとすれば,テイラー展開の1次の近似で sinθθ であるから

d2θdt24π2θ

= に変えて方程式にすると,これは単振動の方程式であるから,同じ初期条件で以下のように解析的に解ける。

θs=θ0cos(2πt)=π4cos(2πt)

また,このとき

θ˙s= π22sin(2πt)

この単振動の解も重ねてグラフにしてみよ。周期にはどのような違いがあるか。

In [25]:
# ax を使う際の最初のおまじない
fig, ax = plt.subplots()

# 数値解 x は .y[0] のようにして取り出す

# 縦軸の値をラジアンから度に変換
theta = np.degrees(1)*(answer2.y[0])
# 数値解
ax.plot(answer2.t, theta, label='数値解')

# 同じ初期条件の単振動解
ax.plot(t_list, 45*np.cos(2*np.pi*t_list), label = '単振動解')

ax.grid()
ax.set_xlabel("$t$")
ax.set_ylabel(r"$\theta$ (度)")
ax.legend();

参考:ルンゲ・クッタ法による2階常微分方程式の数値解法

さて,SymPy で自作のルンゲ・クッタ法 rk2() を使った数値解法の際は,かなり時間がかかりました。

NumPy の関数を使った場合はどうでしょうか。

In [26]:
def rk2(Func, x0, v0, t0, t1, N):
    """ d^2 x/dt^2 = Func(t, x, v) を
        初期条件 t = t0 で x = x0, v = v0 とし,
        [t0, t1] を N 分割して
        4次のルンゲ・クッタ法で解く
    """
    h = (t1 - t0)/N
    T = [t0]
    X = [x0]
    V = [v0]
    t = t0
    x = x0
    v = v0
    for i in range(N):
        k1 = h*v
        m1 = h*Func(t, x, v)
        k2 = h*(v+m1/2)
        m2 = h*Func(t+h/2, x+k1/2, v+m1/2)
        k3 = h*(v+m2/2)
        m3 = h*Func(t+h/2, x+k2/2, v+m2/2)
        k4 = h*(v+m3)
        m4 = h*Func(t+h, x+k3, v+m3)
        x = x + (k1 + 2*k2 + 2*k3 + k4)/6
        v = v + (m1 + 2*m2 + 2*m3 + m4)/6
        t = t + h 
        T.append(t)
        X.append(x)
        V.append(v)
    return [T, X, V]
In [27]:
def F(t, x, v):
    # NumPy の関数・定数を使う
    return -4*np.pi**2*np.sin(x)

x0 = np.radians(45)
v0 = 0
t0 = 0
t1 = 3

for N in [100, 200, 400, 800, 1600, 3200]:
    T2, X2, V2 = rk2(F, x0, v0, t0, t1, N)
    # X2 の最後の要素 X2[-1] を表示
    print("N =%5d, x(t1) =%15.12f" % (N, X2[-1]))
N =  100, x(t1) = 0.591465894161
N =  200, x(t1) = 0.591544091180
N =  400, x(t1) = 0.591548830281
N =  800, x(t1) = 0.591549121191
N = 1600, x(t1) = 0.591549139198
N = 3200, x(t1) = 0.591549140318
In [28]:
%%time

T2, X2, V2 = rk2(F, x0, v0, t0, t1, 1600)
CPU times: user 12.1 ms, sys: 14 µs, total: 12.1 ms
Wall time: 11.9 ms

rk2() は同じコードですので,F(t, x, v) に NumPy の関数を使っただけで,約 600 倍 (!) 程度高速になっていることがわかります。

ですので,解析的に計算をしたいなら SymPy ですが,数値計算を高速にしたいなら NumPy,という具合に使い分けが肝心です。

○練習:単振り子の周期

数値計算で求めた単振り子の周期を調べるため,θ˙= y[1] のグラフを描いてみます。比較のために,単振動の場合の θ˙s も一緒にグラフにします。

In [29]:
# ax を使う際の最初のおまじない
fig, ax = plt.subplots()

# 数値解
ax.plot(answer2.t, answer2.y[1], label='数値解')

# 同じ初期条件の単振動解
ax.plot(t_list, -np.pi**2/2*np.sin(2*np.pi*t_list), label = '単振動解')

ax.xaxis.set_major_locator(MultipleLocator(0.5))
ax.set_xlim(0, 5)

ax.grid()
ax.set_xlabel("$t$")
ax.set_ylabel(r"$\dot{\theta}$")
ax.legend();

上のグラフから,単振動の場合は t=0θ˙=0 から始まり,t=0.5 で折り返し,t=1 で再び θ˙=0 となることから,周期が 1 であることがわかります。

単振り子の数値解の場合は t=1 ではなく,1<t1<1.5 にあることがわかります。この値を数値的に求めてください。

ヒント:以下のように,x を入れると t= 0 から t= x まで数値的に微分方程式を解き,最後の値 θ˙(x) を返す関数を dottheta(x) として定義してみます。

In [30]:
def F(t, y):
    theta = y[0]
    dottheta = y[1]
    F1 = dottheta
    F2 = -4*np.pi**2*np.sin(theta)
    return [F1, F2]

def dottheta(x):
    tini = 0
    tend = x
    t_list = np.linspace(tini, tend)
    # 初期条件 theta0 はラジアンになおして
    y_ini = [np.radians(theta0), 0]

    answer2 = solve_ivp(F, [tini, tend], y_ini, t_eval = t_list, 
                        # デフォルトの許容誤差は大きめなので
                        # 小さい値に設定
                        rtol = 1.e-12, 
                        atol = 1.e-15)
    return (answer2.y[1])[-1]
In [31]:
theta0 = 45
ans1 = root_scalar(dottheta, bracket = [1, 1.5])
t1 = ans1.root
print(t1)
1.0399733431967593

念のため,次に「ひとまわり」して戻ってくる時刻は 2<t2<2.5 にあります。t2 を求め,時間間隔 t2t1t1 が等しいかどうか,確認してください。

In [32]:
ans2 = root_scalar(dottheta, bracket = [2, 2.5])
t2 = ans2.root
print('t1      =', t1)
print('t2 - t1 =', t2 - t1)
t1      = 1.0399733431967593
t2 - t1 = 1.039973343196681

念には念をいれて,0<t<5 の範囲で θ˙(t) が正の側からゼロになる時刻 t1,t2,t3, を求め,時間間隔が一定であるか(つまり,単振り子の運動が確かに周期的であるか)確認してください。

In [33]:
t0 = [0]
for i in range(1, 5):
    ans = root_scalar(dottheta, bracket = [i, i+0.5], xtol=1e-15)
    t0.append(ans.root)
    print('t%1d - t%1d = %15.13f' %(i, i-1, t0[i]-t0[i-1]))
t1 - t0 = 1.0399733431968
t2 - t1 = 1.0399733431967
t3 - t2 = 1.0399733431966
t4 - t3 = 1.0399733431965

時間間隔が(ほぼ)一定であることが確認できたら,これを t=0 における初期条件 θ0=45, θ˙0=0 の場合の周期とします。

同様にして,t=0 における初期条件

  • θ0=1, θ˙0=0 (単振動解との比較も兼ねて)
  • θ0=60, θ˙0=0
  • θ0=70, θ˙0=0
  • θ0=80, θ˙0=0

の場合の単振り子の周期を求めなさい。この結果から,初期条件としての振幅 θ0 と周期の間にはどのような関係があるか,考察してください。