INFORMATION

リンク先
定数係数の2階微分方程式解法ツール(陽、陰オイラー法、RK4、解析解の比較ができます)を公開しました。(2017/3/12)

常微分方程式の解法

陽解法、陰解法

・使用例(一階微分)

 

まず一番簡単な例として、以下の式を解く場合を考えます。

$$\frac{dx}{dt}=x、x(0)=1$$

まずは解析解を求めます。

$$\frac{1}{x}dx=dt$$

両辺を積分して

$$ln(x)=t+C$$ $$x=C\cdot exp(t)$$

初期条件:x(0)=1から、C=1。よって解析解は

$$x=exp(t)$$

続いて数値計算で解くためには式を離散化する必要があります。

1time step先の計算を、現在値と、傾きから求める陽解法(Explicit Euler法)と、1time step先の値自身を使う陰解法(Implicit Euler)をご紹介します。

各使用例を以下ご紹介します。

<陽解法の場合> $$\frac{x_{i+1}-x_i}{\Delta t}=x_i、x(0)=1$$

式をまとめ直すと、

$$x_{i+1}=(1+\Delta t)x_i、x(0)=1$$ <陰解法の場合> $$\frac{x_{i+1}-x_i}{\Delta t}=x_{i+1}、x(0)=1$$

式をまとめ直すと、

$$x_{i+1}=\frac{x_i}{1-\Delta t}、x(0)=1$$

上式をVBA、Python、Scilabで求めたものは、それぞれ以下をご参照ください。

VBAPythonScilab

・使用例(一階微分&非線形)

<陰解法+Newton-Raphson法の場合>

以下の非線形の一階微分方程式を解くことを考えます。

$$\frac{dy}{dt}=-\left(y^4-1\right)、y(0)=2$$

右辺が非線形になっているため、一度にパッと逆行列を求めることができません。

上記の問題に対してNerton Raphson法を適用した以下の一般化式で求めます。

$$\left(\boldsymbol{I}-h\boldsymbol{J}(t_{n+1},\boldsymbol{Y}_{n+1,i})\Delta\boldsymbol{Y}_{n+1,i+1}\right)=-\left(\boldsymbol{Y}_{n+1,i}-hf\left(t_{n+1},\boldsymbol{Y}_{n+1,i}\right)-\boldsymbol{Y}_{n}\right)$$

以下の著書を参考にしています。

Computational Methods in Engineering, 著者: S.P. Venkateshan、Prasanna Swaminathan

ブログでも紹介したのですが、実用的で分かりやすくて良い本です。

・使用例(二階微分)

もう一例、今度は2階微分です。

$$\frac{d^2x}{dt^2}=-x、x(0)=1、\frac{dx}{dt}=0$$

まずは解析解を求めます。

$$\frac{d^2x}{dt^2}+x=0$$

ラプラス変換すると、

$${\mathfrak{L}}\left[ \frac{dx^2}{dt^2} \right]+{\mathfrak{L}}\left[x\right] = {s^2}F( s ) - sf( 0 ) - f( 0 )+F(s)=0$$

右辺を整理しなおして、

$$ \left({s^2}+1 \right) F( s ) - sf( 0 ) - f( 0 )=0$$

初期条件:x(0)=f(0)=1、x'(0)=f'(0)=0より、

$$ \left({s^2}+1 \right) F( s ) = 1 $$ $$ F( s ) = \frac{1}{\left({s^2}+1 \right)} $$

逆ラプラス変換を行い、以下が解となります。

$$ x( t ) = sin(t) $$

なお、上式の変形はラプラス変換の微分則を使っております。ラプラス変換の詳細については、こちらControlをご参照ください。

ラプラス変換の微分則↓

$${\mathfrak{L}}\left[ {\frac{d^2}{dt^2}f( t )} \right] = {s^2}F( s ) - sf( 0 ) - f( 0 )$$

逆ラプラス変換(三角関数)↓

  $$sin(at)=\frac{a}{s^2+a^2}$$

続いて数値計算で解くためには式を離散化する必要があります。

<陽解法の場合>

u=\(\frac{dx}{dt}\)とすると、\(\frac{d^2x}{dt^2}=\frac{du}{dt}\)

$$\begin{eqnarray} \begin{cases} \frac{x_{i+1}-x_i}{\Delta t}=u_i & \\ \frac{u_{i+1}-u_i}{\Delta t}=-x_i & \end{cases} \end{eqnarray}$$

式をまとめ直すと、

\begin{align*} \begin{bmatrix} x_{i+1} & \\ u_{i+1} & \end{bmatrix} &= \begin{bmatrix} 1 & \Delta t \\ -\Delta t & 1 \end{bmatrix} \begin{bmatrix} x_i & \\ u_i & \end{bmatrix} \end{align*} <陰解法の場合>

u=\(\frac{dx}{dt}\)とすると、\(\frac{d^2x}{dt^2}=\frac{du}{dt}\)

$$\begin{eqnarray} \begin{cases} \frac{x_{i+1}-x_i}{\Delta t}=u_{i+1} & \\ \frac{u_{i+1}-u_i}{\Delta t}=-x_{i+1} & \end{cases} \end{eqnarray}$$

式をまとめ直すと、

\begin{align*} \begin{bmatrix} 1 & -\Delta t \\ \Delta t & 1 \end{bmatrix} \begin{bmatrix} x_{i+1} & \\ u_{i+1} & \end{bmatrix} &=\begin{bmatrix} x_i & \\ u_i & \end{bmatrix} \end{align*}

上式をVBA、Python、Scilabで求めたものは、それぞれ以下をご参照ください。

VBAPythonScilab

・使用例(二階微分&非線形)

<陰解法+Newton-Raphson法の場合>

例としてVan der Poles式を解く方法をご紹介します。

陽解法、陰解法、4次のRunge-Kutta法の精度、安定性比較

・WPF、Silverlightで作成したツール

・各手法で解いた場合の精度、安定性を簡単に比較できるように、ツールを作成しました。

ODEPlot

・係数、初期値、計算刻みを変更し、各手法の精度、安定性を比較してみてください。

高速自動微分

当サイトに不具合、ご意見等ございましたらCEsolutionにお知らせください。

Page Top