とらりもんHOME  Index  Search  Changes  Login

とらりもん - C言語による微分方程式の解析 (Euler法) Diff

  • Added parts are displayed like this.
  • Deleted parts are displayed like this.

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;
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

{{attach_view(euler1.png)}}

上の図で, 緑色の線が解析解(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()
{
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);
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

{{attach_view(11b.png)}}

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

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

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

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

{{attach_view(11c.png)}}

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