INFORMATION
-
リンク先
- 定数係数の２階微分方程式解法ツール（陽、陰オイラー法、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で求めたものは、それぞれ以下をご参照ください。
VBA、Python、Scilab・使用例（一階微分&非線形）
<陰解法+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
ブログでも紹介したのですが、実用的で分かりやすくて良い本です。
・使用例（二階微分）
もう一例、今度は２階微分です。
$$\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で求めたものは、それぞれ以下をご参照ください。
VBA、Python、Scilab・使用例（二階微分&非線形）
<陰解法+Newton-Raphson法の場合>例としてVan der Poles式を解く方法をご紹介します。
陽解法、陰解法、４次のRunge-Kutta法の精度、安定性比較
・WPF、Silverlightで作成したツール
・各手法で解いた場合の精度、安定性を簡単に比較できるように、ツールを作成しました。
・係数、初期値、計算刻みを変更し、各手法の精度、安定性を比較してみてください。
Runge-Kutta法
・使用例（一階微分）
当サイトに不具合、ご意見等ございましたらCEsolutionにお知らせください。