このサイトは更新終了し, こちらに移行しています。

数値解析入門

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

筑波大学農林工学系 奈佐原顕郎

Euler法


dx/dt = 2x をEuler法で解こう(t =0のときx =1とする)。課題6-1と同じ解析である:
  1. /* euler1.c
  2. * 2004/02/15 K. Nishida
  3. */
  4. # include <stdio.h>
  5. main()
  6. {
  7. double t=0;
  8. double dt=0.1;
  9. double x=1;
  10. double dx;
  11. printf("%lf %lf\n", t, x);
  12. for (t=dt; t<=1; t=t+dt)
  13. {dx=2*x*dt;
  14. x=x+dx;
  15. printf("%lf %lf\n", t, x);
  16. }
  17. }
注意: 左端の数字(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)を適当に設定して解いてみよう。以下にその例を示す。
  1. /* euler2.c
  2. * spring motion
  3. * 2004/02/15 K. Nishida
  4. */
  5. # include <stdio.h>
  6. main()
  7. {
  8. double t=0;
  9. double dt=0.01;
  10. double x=1, v=0;
  11. double dx, dv;
  12. double k=1;
  13. double m=1;
  14. printf("%lf %lf %lf\n", t, x, v);
  15. for (t=dt; t<=10; t=t+dt)
  16. {dx=v*dt;
  17. dv=-1*x*(k/m)*dt;
  18. x=x+dx;
  19. v=v+dv;
  20. printf("%lf %lf %lf\n", t, x, v);
  21. }
  22. }
注意: 左端の数字(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 前問を解析的に解いて, その解析解と, 前問の数値解をいっしょにグラフに表示してくらべてみよ。

(ブラウザの「戻る」ボタンで戻ってください。)