Sageでグラフを再現してみよう：微分方程式（オイラー法）編
この企画は、雑誌や教科書にでているグラフをSageで再現し、 グラフの意味を理解すると共にSageの使い方をマスターすることを目的としています。
今回は、島根大学電子制御システム工学科の吉田和信氏の Octaveによる動的システムシミュレーション入門 の例題2.1の図2.1を題材にします。
例題2.1について
以下のような連立一次常微分方程式をオイラー法を使って解く問題です。 $$ \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, uは以下の様な行列とベクトルで与えられています。刻み幅$\Delta t$ = 0.05, 終了時刻$t_f$ = 20とします。 $$ \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 $$
octave連携の準備
sageからoctaveを使う場合のグラフ処理を助ける関数を定義します。
|
|
octaveを使ったオイラー法の解析例
最初は、octaveを使った解析例（解ex2_1.m）をそのままsageで実行します。ex2_1.mはワークシート上部のDataタグで内容を確認することができます。
ベクトル$x$は、$ z, \dot z $を要素に持ちます。 $$ x(t) = \left[\begin{array}{c} z(t)\\ \dot z (t) \end{array}\right] $$
図の青の線が、$z(t)$、緑の線が$\dot z(t)$に相当し、$z(t)$のピークで$\dot z(t)$が0になっているのが分かります。
|
|
Sageでのオイラー法の解析例
次にSageを使って同じEuler法で解いてみます。ここではSageの特徴である、インタラクティブ機能を使って、 刻み幅dtを0.05, 01, 0.5, 1から選択できるようにしました。刻み幅によってグラフがどのように変化するかをみてください。
オイラー法は、連立一次常微分方程式を数値微分を使って解く方法であり、刻み幅によって解が正しく求まらないことがあることに注意が必要です。
Click to the left again to hide and once more to show the dynamic interactive window |
|||||||||||||||
Sageの関数を使った例
Sageにも連立一次常微分方程式を解く関数がいくつか提供されています。ここでは、 desolve_system_rk4を使った例を示します。
例題を$ z, \dot z $で表すと、 $$ \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 $$
これを展開すると、以下の様になります。 $$ \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)]$を与えて解きます。
|
|
伝達関数H(s)を使った解法
$\ddot z, \dot z, z, u$をまとめると以下の２次の微分方程式になります。 $$ \ddot z + \dot z + z = u(t) $$ 上記の両辺をラプラス変換します。zのラブラス変換をZ, uのラブラス変換をUとすると以下の様になります。 $$ s^2 Z + s Z + s = 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を使って行ってみます。
\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{1}{3} \, {\left(\sqrt{3} \sin\left(\frac{1}{2} \, \sqrt{3} t\right) + 3 \, \cos\left(\frac{1}{2} \, \sqrt{3} t\right)\right)} e^{\left(-\frac{1}{2} \, t\right)} + 1
\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{1}{3} \, {\left(\sqrt{3} \sin\left(\frac{1}{2} \, \sqrt{3} t\right) + 3 \, \cos\left(\frac{1}{2} \, \sqrt{3} t\right)\right)} e^{\left(-\frac{1}{2} \, t\right)} + 1
|
結果のプロット
求まった解をプロットすると図2.1と同じ結果となり、解析解と数値計算の結果が一致することが確認できました。
|
|
微分方程式の解
最後に数式処理システムSageの微分方程式を解く関数desolveを使った例を示します。
desolveでは、変数tをvar関数で、変数tの関数ｚをfunction関数で定義し、 微分方程式deに$\ddot z + \dot z + a = 1$の関係式をセットし、 初期条件icsに$[ t_0, z(t_0), \dot z(t_0)]$を与えて解きます。
解もグラフも先の伝達関数の例と同です。
\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{1}{3} \, {\left(\sqrt{3} \sin\left(\frac{1}{2} \, \sqrt{3} t\right) + 3 \, \cos\left(\frac{1}{2} \, \sqrt{3} t\right)\right)} e^{\left(-\frac{1}{2} \, t\right)} + 1
\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{1}{3} \, {\left(\sqrt{3} \sin\left(\frac{1}{2} \, \sqrt{3} t\right) + 3 \, \cos\left(\frac{1}{2} \, \sqrt{3} t\right)\right)} e^{\left(-\frac{1}{2} \, t\right)} + 1
|
|
|
まとめ
行列とベクトルで表される連立一次常微分方程式をオイラー法で数値微分することによって、 n次（今回は2次）の常微分方程式の解を求めることができる。
さらにSageを使った解析解の求め方を紹介し、微分方程式をSageを使ってどのように解くかを図2.1 を例に説明しました。
是非、この機会にSageを使っていろいろな微分方程式を解いてみてください。
|
|