2024.7.1 solve_ivpのデータ評価点自動設定１【scipy】
solve_ivpはscipyが提供する常微分方程式の初期値問題を解くための関数である。おそらく「IVP」は initial value problemの略。
一般的なオイラー法やルンゲクッタ法では、時間の刻み幅は一定であるのに対して、solve_ivpではデータの評価を行う時刻を必要に応じて自動設定する。そのことについて調べてみた。
code:ivp1.py
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
def func(t, x):
if t < 5:
u = 0
else:
u = 1
return -x + u
ts, yf = 0, 10 # 時間区間
sol = solve_ivp(func, trange, x_init) # t_evalを与えないので、評価時刻は自動設定
t_hist = sol.t
plt.grid()
plt.plot(t_hist, x_hist, '*-')
plt.show()
しかしながら、上のような通常の使い方では関数funcの中で生成される値（t や u）は保持されない。これらをグローバル変数に保存してみよう。
code:ivp2.py
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
def func(t, x):
global tu_hist, t_hist
if t < 5:
u = 0
else:
u = 1
tu_hist.append(t)
u_hist.append(u)
return -x + u
ts = 0
tf = 10
u_hist, tu_hist = [], []
sol = solve_ivp(func, trange, x_init)
t_hist = sol.t
fig = plt.figure()
ax = fig.add_subplot(2,1,1)
ax.set_title('x')
ax.grid()
ax.plot(t_hist, x_hist, '*-')
ax = fig.add_subplot(2,1,2)
ax.set_title('u')
ax.grid()
ax.plot(tu_hist, u_hist, '*-')
plt.show()
結果、上のグラフは状態x、下のグラフは入力uである。
https://scrapbox.io/files/66821558759463001cefb307.png
入力が意図した結果と違う。どういうことであろうか。
状態と入力が生成された時刻をそれぞれプロットしてみる。
code:(続き)ivp2.py
fig = plt.figure()
ax = fig.add_subplot(2,1,1)
ax.set_title('time for x')
ax.grid()
ax.plot(t_hist, '*-')
ax = fig.add_subplot(2,1,2)
ax.set_title('time for u')
ax.grid()
ax.plot(tu_hist, '*-')
plt.show()
https://scrapbox.io/files/6682162c827481001d89154a.png
これらのグラフの横軸はデータ点数である、縦軸は上が x のデータ点に対する時刻、下が入力 u に対する時刻である。
状態のデータ点が16であるのに対して、入力のデータ点は140であることに注意、これだけの回数funcは呼び出されている。
time for x の縦軸より、状態の動きに精度が求められる領域では、時間刻みが細かくなっている。
time for uの縦軸より、funcの呼出しは時系列順ではなく、必要に応じて時間の経過速度が変化したり、場合によっては戻っていることがわかる。