システムとモデリング

modelica, Julia, Design Structure Matrix, SysML, 他モデリング全般について。

微分方程式を行列で差分化して解く?

ストラング先生の計算理工学に、行列で微分方程式を差分化して解く方法が載っていました。 浅学のためこういった解法は初めてであまり理解もできていませんが、練習してみたいと思います。

世界標準MIT教科書 ストラング:計算理工学

世界標準MIT教科書 ストラング:計算理工学

例として下記のような微分方程式を解いてみます。

$$ \begin{equation} \frac{dy}{dx}=3x^{2}, y(0) = 0 \end{equation} $$

この方程式の解は y = x^{3}ですが、これを行列で差分化して確認してみます。

yを前進差分で差分化する行列は例えば以下のようなものになります。 $$ A= \begin{bmatrix} 1 & 0 & 0 & 0 \\ -1 & 1 & 0 & 0 \\ 0 & -1 & 1 & 0 \\ 0 & 0 & -1 & 1 \end{bmatrix} $$

これにyをかけると以下のようになります。このように前進差分を取ることができます。 $$ \begin{bmatrix} 1 & 0 & 0 & 0 \\ -1 & 1 & 0 & 0 \\ 0 & -1 & 1 & 0 \\ 0 & 0 & -1 & 1 \end{bmatrix} \begin{bmatrix} y(0) \\ y(1) \\ y(2) \\ y(3) \end{bmatrix}  = \begin{bmatrix} y(0) \\ y(1) - y(0) \\ y(2) - y(1) \\ y(3) - y(2) \end{bmatrix} $$

今回は右辺が 3x^{2}なので $$ \boldsymbol{y} = \begin{bmatrix} y(0) \\ y(1) \\ y(2) \\ y(3) \end{bmatrix} $$ とすると 行列を使った方程式は下記のようになります。 $$ \begin{bmatrix} 1 & 0 & 0 & 0 \\ -1 & 1 & 0 & 0 \\ 0 & -1 & 1 & 0 \\ 0 & 0 & -1 & 1 \end{bmatrix} \boldsymbol{y}  = \begin{bmatrix} 0 \\ 3 \\ 12 \\ 27 \end{bmatrix} $$

これを解くと $$ \boldsymbol{y}  = \begin{bmatrix} 0 \\ 3 \\ 15 \\ 42 \end{bmatrix} $$ が得られます。ただしこの値は y=x^{3}とはまったく一致しません。 おそらくステップ数が4しかないために誤差が大きかったものと思われます。

次はステップ数を増やして解法を確認してみます。

今回はここまでになります。

続き

otepipi.hatenablog.com