pypy.com/

Python、Unity、FX自動化などを勉強しています。あと、コーラと車も好きです。そこらへんについて、たまに記事を書きます。

c言語で微分方程式を解く(修正オイラー法)

前回の記事では微分方程式を解く方法として、オイラー法を紹介しましたが、さらに精度の良い方法として修正オイラー法があるので今回はそれについて書きたいと思います。
前回の記事は、以下のリンクよりアクセスできるので、気になる方はご覧ください。

taiboman.com - taiboman リソースおよび情報

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

修正オイラー法の考え方

前回と同じように、
 \frac{dy}{dx} = f(x,y) \tag{1}
について考えます。
f:id:SaidaTaisei:20200119153219p:plain
上図について、赤線の傾きは(1)式より、f(x,y) となります。
ここまでは、オイラー法と同様で、修正オイラー法になると増えるのは緑線の傾きです。
赤線の傾きが f(x,y) であるから、x=x+h での y の値は
\begin{align} 
y(x+h)=y(x)+hf(x,y)
\end{align} 
と書くことができます。
ここで、k_1=hf(x,y) とすると、
y(x+h)=y(x)+k_1 \tag{2}
となります。
(2)式をもとにすると、x=x+h における傾き(緑線の傾き)は、
\begin{align} 
f(x+h,y(x)+k_1)
\end{align} 
この緑線の傾きをもとに、x=x+hy の値を計算すると、
\begin{align} 
y(x+h)=y(x)+hf(x+h,y(x)+k_1)
\end{align} 
ここで、k_2=hf(x+h,y(x)+k_1)とすると
y(x+h)=y(x)+k_2 \tag{3}
上の図のように、一般的には、求めたい値は赤線と緑線の間に存在するため、
修正オイラー法では、
y(x+h)=y(x)+ \frac{1}{2}(k_1+k_2) \tag{4}
とします。
これを続けることで、(1)式の解を得ることができます。

実際にコードを書いてみる

それでは、コードを書いていきます。
今回は、前回のオイラー法と同じように、
\begin{align} 
\frac{dy}{dx} = xy \tag{5}
\end{align}
を初期値(x_0,y_0)=(0.0,4.0)
分割数100で解きたいと思います。
なお、(5)式の厳密解は
\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++){
    double k1 = h*f(x,y);
    double k2 = h*f(x+h,y+k1);
    y = y + (k1+k2)/2.0;
    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);
}

出た値をグラフにしてみると以下のようになります。
f:id:SaidaTaisei:20200119211352p:plain
見るとわかるように、グラフがほとんど重なっており、精度がオイラー法より良くなったことがわかります。
実際にx=3.0の時の値を見ると厳密解が 360.06853 で求めた値が 359.03017 であるため、そこまでの精度を求めない場合には充分であると思います。
今回は以上で終わります。