c言語で微分方程式を解く（オイラー法）

お久しぶりです。久しぶりにブログを書きます。

今回は常微分方程式を解くということをテーマで書きます。

普段、unityなんかを使われる方はあまり意識しなくても、力なんかを与えてやると物理エンジンが勝手に計算してくれるのですが、自分でも解けたほうがいいだろうということで、解いていこうと思います。

今回使う手法はオイラー法という一番シンプルで簡単な方法です。

精度はあまりよくありませんが、簡単なものであればこの方法を使うこともできると思います。

それでは書いていきます。

追記
より精度よく計算できる修正オイラー法について記事を書きました。
興味のある方は是非ご覧ください。
c言語で微分方程式を解く（修正オイラー法） - pypy.com/

さらに精度よく計算できるルンゲクッタ法の記事を書きました。
c言語で微分方程式を解く（ルンゲクッタ法） - pypy.com/

二階微分方程式のも書きました。
c言語で二階常微分方程式を解く（オイラー法） - pypy.com/

オイラー法の仕組み
コード

オイラー法の仕組み

$\frac{dy}{dx} = f(x,y)　\tag{1}$

上式について考えます。

グラフが以下のような形になったとします。

この時、あるxとそこからh増加したx+hとそれに対応するy(x)、y(x+h)を考えます。

以下の図のような感じです。

上図の赤線の傾きは(1)より、 $f(x,y)$ であるから、

y成分の増加量は $hf(x,y)$ であり、

$y(x+h)$ は以下の式に近似できる。

$y(x+h)　\approx　y(x)+hf(x,y)\tag{2}$

ここまで、図を用いて説明を書きましたが、もっと単純にテイラー展開を用いた場合、

\begin{align}
y(x+h) & =y(x)+\frac{h}{1!}\frac{dy}{dx}+\frac{h^2}{2!}\frac{d^2y}{dx^2}+\cdots\\\\
& \approx　y(x)+hf(x,y)
\end{align}

となり、(2)式と同じ結果になります。

上の図では、結構大きな誤差があるように見えますが、hを小さくすれば誤差は小さくなります。

これを続けることで、(1)式の解を得ることができます。

コード

それでは、コードを書いていきます。
今回は、
\begin{align}　
\frac{dy}{dx} = xy　\tag{3}
\end{align}
を初期値 $(x_0,y_0)=(0.0,4.0)$
分割数100で解きたいと思います。
なお、(3)式の厳密解は
\begin{align}　
y = 4e^ \frac{x^2}{2}　
\end{align}
です。
厳密解との比較もしたいと思います。

#include <stdio.h>
#include <math.h>

#define N 100 //分割数

double f(double x,double y);
double y_true(double x);

int main(void){
  //初期値
  double x=0.0;
  double y=4.0;
  double h=(3.0-0.0)/(double)N; //刻み幅

  printf("x       y       y_true\n");
  for(double i=0; i<N; i++){
    y = y + h*f(x,y);
    x = x + h;
    printf("%5.5lf, %5.5lf, %5.5lf\n",x,y,y_true(x));
  }
}

double f(double x,double y){
  return x*y;
}
double y_true(double x){
  return 4.0*exp(x*x/2.0);
}