C言語による微分方程式の解析 (Euler法)

dx/dt = 2x をEuler法で解こう(t =0のときx =1とする)。課題6-1と同じ解析である：

/* euler1.c
* 2004/02/15 K. Nishida
*/
# include <stdio.h>
main()
{
double t=0;
double dt=0.1;
double x=1;
double dx;
printf("%lf %lf\n", t, x);
for (t=dt; t<=1; t=t+dt){
  dx=2*x*dt;
  x=x+dx;
  printf("%lf %lf\n", t, x);
  }
}

注意: 左端の数字（6.など）は打ち込まないで良い。

$ gcc euler1.c -o euler1
$ ./euler1 > euler1.txt
$ gnuplot
GNUPLOT> plot "euler1.txt" w p, exp(2*x) w l

上の図で, 緑色の線が解析解（x=exp(2t)）, 赤点がオイラー法による数値解である。

課題13-1 上を参考にして, dx/dt = 2x (1- x) をEuler法で解け(t =0のときx =0.1とする)。課題5-2と同じ解析である。

課題13-2 以下のLotka-Volterra方程式を解き, 結果をGNUPLOTで表示せよ：

   dS/dt = pS - qWS
   dW/dt = -rW + sWS 

t =0でS =2, W =0.5とする。p =q =r =s =1とする。

Newtonの運動方程式は,

   d2x / dt2 = F/m   (xは質点の座標, Fは力, mは質量) 

であり, 2階の微分方程式なので, 直接的には積分できない。そこで, v=dx/dtという変数を導入し, xとvの2つについての1階の連立微分方程式に変換する：

   dv / dt = F/m
   dx / dt = v 

あとは, 上のLotka-Volterra方程式などと同様にEuler法で解くことができる。 たとえば, バネの振動運動を解いてみよう。バネ定数kのバネに質量mのおもりがついて振動しているときの運動方程式は,

   dv / dt = - x * k/m
   dx / dt = v 

となる。初期条件および, パラメータ(ここではmとk)を適当に設定して解いてみよう。以下にその例を示す。

/* euler2.c
* spring motion
* 2004/02/15 K. Nishida
*/
# include <stdio.h>
main(){
double t=0;
double dt=0.01;
double x=1, v=0;
double dx, dv;
double k=1;
double m=1;
printf("%lf %lf %lf\n", t, x, v);
for (t=dt; t<=10; t=t+dt){
  dx=v*dt;
  dv=-1*x*(k/m)*dt;
  x=x+dx;
  v=v+dv;
  printf("%lf %lf %lf\n", t, x, v);
  }
}

注意: 左端の数字（6.など）は打ち込まないで良い。

$ gcc euler2.c -o euler2
$ ./euler2 > euler2.txt
$ gnuplot
GNUPLOT> plot "euler2.txt" using 1:2 title "x" w l, "euler2.txt" using 1:3 title "v" w l

課題13-3 おもりに空気抵抗が働く場合を考えて, 上で作ったバネの振動のプログラムを改造せよ。 ヒント：空気抵抗は,

   dv / dt = - x * k/m - v * a    (aは抵抗係数というパラメータ) 

という形で運動方程式に入る。aに適当な値を設定して18行目をそのように改造すれば良い。

以下にこの課題の解の例を示す：

課題13-4 前問を解析的に解いて, その解析解と, 前問の数値解をいっしょにグラフに表示してくらべてみよ。