微分方程式

Sageを使って微分方程式を解く方法を紹介します。微分方程式の性質から解析的な解が求まらない場合等に備え、 数値的に微分方程式を方法についても取り上げます。

desolveを使った微分方程式の解法

最初にSageで微分方程式を解くときに一般的に使用されるdesolve関数の使い方について説明します。

例として以下のような2次の微分方程式をSageで表してみます。 $$ \ddot z + \dot z + z = u(t) $$ ここで外部からの変位u(t)を、$u(t) = 1, t \ge 0$のステップ関数とすると上記の微分方程式は 以下のように定義されます。

desolveでは、変数tをvar関数で、変数tの関数ｚをfunction関数で定義し、 微分方程式deに$\ddot z + \dot z + a = 1$の関係式をセットし、 初期条件icsに$[ t_0, z(t_0), \dot z(t_0)]$を与えて解きます。

求まった解$z(t)$とその微分$\dot z(t)$をプロットしてみます。

$z(t)$は、最初は振動していますが、t=10になると定常状態に収まって１となっています。

伝達関数H(s)を使った解法

次に、ラプラス変換と伝達関数を使って同じ微分方程式を解いてみましょう。 $$ \ddot z + \dot z + z = u(t) $$ の両辺をラプラス変換します。zのラブラス変換をZ, uのラブラス変換をUとすると以下の様になります。 $$ s^2 Z + s Z + Z = U $$ ここで伝達関数$H(s) = \frac{Z}{U}$を計算します。 $$ H(s) = \frac{1}{s^2 + s + 1} $$ 単位ステップ応答のラブラス変換は、$\frac{1}{s}$なので、$Z(s)$は以下の様に求めまります。 $$ Z(s) = H(s) U(s) = H(s) \frac{1}{s} = \frac{1}{s^3 + s^2 + s} $$ 最後に$Z(s)$をラプラス逆変換すると$z(t)$の解析解が求まります。

上記の処理をSageを使って行ってみます。

結果のプロット

求まった解をプロットするとdesolveの同じ結果と同じ結果を得ることができました。

微分方程式の数値解法

desolve_system_rk4を使った解法

desolve_system_rk4は、ルンゲクッタ法を使って効率よく微分方程式の数値解をもとめる関数です。

例題の微分方程式を$\ddot z, \dot z$について見てみます。 $$ \ddot z + \dot z + z = u(t) $$ 以下のような等価な連立微分方式を得ることができます。 $$ \left\{\begin{array}{l} \dot z = \dot z\\ \ddot z = - z - \dot z + u(t) \end{array} \right. $$

これをdesolve_system_rk4関数を使って定義すると以下のようになります。 初期条件icsは$[ t_0, z(t_0), \dot z(t_0)]$として解きます。

連立微分方程式の行列とベクトルを使った表現

先の連立微分方程式を行列とベクトルで表現すると以下のようになります。 $$ \left[\begin{array}{c} \dot z\\ \ddot z \end{array}\right] = \left[\begin{array}{cc} 0 & 1\\ -1 & -1 \end{array}\right] \left[\begin{array}{c} z\\ \dot z \end{array}\right] + \left[\begin{array}{c} 0\\ 1 \end{array}\right] u(t), \left[\begin{array}{c} z\\ \dot z \end{array}\right] = \left[\begin{array}{c} 0\\ 0 \end{array}\right] , u(t) = 1, t \ge 0 $$

ここでベクトル$x(t)$を以下のように定義します。 $$ x(t) = \left[\begin{array}{c} z(t)\\ \dot z (t) \end{array}\right] $$ $x(t)$で連立微分方程式を表すと以下のように簡単になります。 $$ \dot x(t) = \left[\begin{array}{cc} 0 & 1\\ -1 & -1 \end{array}\right] x(t) + \left[\begin{array}{c} 0\\ 1 \end{array}\right] u(t), x(0) = \left[\begin{array}{c} 0\\ 0 \end{array}\right] , u(t) = 1, t \ge 0 $$

オイラー法を使った解法

オイラー法は、以下のように行列とベクトルで表現された連立一次常微分方程式を差分近似式を使って解く方法です。 $$ \dot x(t) = A x(t) + B u(t), x(0) = x_0 $$ これを差分近似で表現すると以下のようになります。 $$ x(t + \Delta t) \simeq x(t) + (A x(t) + B u(t)) \Delta t $$

それでは、オイラー法を使って先の微分方程式を解いてみましょう。 A, Bを以下のように定義し、刻み幅$\Delta t$ = 0.05, 終了時刻$t_f$ = 20で計算します。 $$ A = \left[\begin{array}{cc} 0 & 1\\ -1 & -1 \end{array}\right] $$ $$ B = \left[\begin{array}{c} 0\\ 1 \end{array}\right] $$