連立常微分方程式の解き方

Question

#####要旨#####

連立線形常微分方程式
dX/dt = AX + b
のAの各要素の値が時間的に変化する場合の、
Pythonの解き方を知りたいです。

#####詳細#####

連立線形常微分方程式

dX/dt = AX + b

　＃X=(x1, x2, ・・・)

において，

行列Aが時間的に変化しない場合は

scipy.integrate.odeintなどで解くことができると理解しております．

行列Aの各要素の値が時間的に変化する場合，

どのように解くのでしょうか．

イメージとしては，

A|t=t1(t1の時のA)　→　［［1,2,3］［4,5,6］［7,8,9］］

A|t=t2　→　［［0,1,2］［3,4,5］［6,7,8］］

・・・

といった具合で時間と共に変化します．

　＃要素の値は予め与えることができます．(計算する前にわかります)

行列Aは疎行列を想定してますが，

3重対角行列でも対称行列でもありません．

なにかライブラリを使って計算できるものでしょうか．
舌足らずかもしれず大変恐縮ですが、よろしくお願いいたします．

Comments

[aokikenichi]

すみません理解が追い付てないのかと思いますが
Xも時間の関数となる変数で、Aも時間の関数となる変数ということでしょうか
これだともう常微分方程式ではないと思いますが

[tmyuki_0716]

理解不足でした、ご指摘の通りだと思います。

dx/dt=ax＋b
の
a=1 t≦t1
a=2 t＞t1
のような問題の連立ver.が解けないか、
と思いました。

[Penpen7]

常微分方程式ではあると思いますが

Answer 1

適当にAに当てはめてやってみました(t依存)。手計算出来ないと思うので、あっているかどうかはわかりません。

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

def A(t):
    return np.array([[np.sin(t), np.cos(t), 0.1],
                     [np.cos(t), -0.5, np.exp(-t)],
                     [np.sin(t), -1, 0]])
def f(x, t):
    return A(t)@x
t0 = np.linspace(0, 10, 1000)
x = odeint(f, [0, 1, 2], t0)

plt.figure()
plt.plot(t0, x[:,0], label="$x_1$")
plt.plot(t0, x[:,1], label="$x_2$")
plt.plot(t0, x[:,2], label="$x_3$")
plt.xlabel("t")
plt.ylabel("x")
plt.legend()
plt.show()

Comments

[tmyuki_0716]

ありがとうございます。
引数funcを、このように定義するのですね！勉強になりました。
結果ですが、Excelを使った計算と殆ど一致することを確認しました。
(微妙に差がありますが、解法の差だと理解してます。)
大変ありがとうございました。

また、コーディングもスマートで、大変勉強になりました。