B..2.1.2 solve_ivp() の利用

scipy には、ODEPACK 由来でない、より新しい関数群も用意されている。 例えば solve_ivp(func, t_span, y0, [,method, t_eval,$\cdots$]]) を使って初期値問題が解ける。

• func は $f$ を計算する関数であるが、引数並びに注意が必要である。 solve_ivp() のときは func(t, x[, param]) という順であり、 odeint() のときの func(x, t [, param]) とは違っている (t と x が逆)。
• t_span は、 [t0, T] のような、初期時刻と最終時刻からなるリストである。

  区間をどのように分割するかは、ユーザーが指定するのではなく、 自動的に決められる。 そのかわり、以下のサンプル・プログラムでは、 rtol=, atol= で許容される誤差を指定する。 内部で、いわゆる刻み幅の自動調節が行われる。

詳しいことは scipy.integrate.solve_ivp -- ScyPy v1.11.3 Manual を見よ。

malthus3.py

# malthus3.py --- dx/dt=a*x, x(0)=x0 を複数の初期値について解いて解曲線を描く

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

def malthus(t, x, a): # 順番に注意 f(t,x,パラメーター)
    return a * x

a=1
t0=0; T=1; n=10
t=np.linspace(t0,T,n+1)

for i in range(0,10):
    x0=[i*0.2] # スカラーでない
    sol=solve_ivp(malthus, [t0,T], x0, args=(a,),
                  dense_output=True, rtol=1e-6,atol=1e-8)
    x=sol.sol(t)
    # 解曲線描画
    plt.plot(t,x.T, label='x0='+'{:6.1f}'.format(x0[0]))

plt.ylim(0, 5)  # 縦軸の範囲指定
    
# タイトル、凡例の位置、横軸と縦軸の説明、グリッド
plt.title('Malthus: dx/dt=ax, x(0)=x0; a='+str(a))
plt.legend(loc='upper left')
plt.xlabel('t')
plt.ylabel('x')
plt.grid()

plt.show()

図 11: malthus3.py で複数の初期値に対する解の解曲線を描いてみる