数値積分の実行

ベルレ法

シミュレションではってかします めるためにはからった利用しますN  あるとき
d2ri   F-i
dt2  = mi    (i = 1,2, ...,N )
(1)

となり
dri
dt--= vi   (i = 1,2, ...,N )
(2)

けますこれらのN  ありますがなのでi  しま

t + Δt  でのr (t + Δt )  すると
                      dr(t)   Δt2-d2r-(t)
r (t + Δt ) = r(t) + Δt   dt  +  2    dt2  + ⋅⋅⋅
(3)

この(1)(2) 使うと
                            Δt2 F (t)
r(t + Δt ) = r(t) + Δtv (t) +---------+  O((Δt )3)
                             2   m
(4)

となりr(t - Δt)
                               2
r(t - Δt ) = r(t) - Δtv (t) + Δt-F-(t)+  O((Δt )3)
                             2   m
(5)

となります(4) (5) ると
                                   2F-(t)-         4
r(t + Δt) + r(t - Δt) = 2r (t) + Δt   m   + O ((Δt) )
(6)

                                          4
r(t + Δt ) - r (t - Δt ) = 2Δtv (t) + O ((Δt) )
(7)

となりますこのでは     4
(Δt )   していますこの 2 から
                                    F (t)
r(t + Δt ) = 2r(t) - r(t - Δt ) + Δt2----
                                     m
(8)

v (t) = --1- {r(t + Δt) - r(t - Δt)}
       2 Δt
(9)

られますこのベルレばれています

ベルレ法ではられるがずれているので(8) (9) 以下のように します
                                2
                             Δt--F-(t)
r (t + Δt ) = r(t) + Δtv (t) + 2  m
(10)

                   Δt--
v(t + Δt) = v (t) + 2m  {F (t) + F (t + Δt)}
(11)

このベルレばれています

ベルレ法のプログラミング

られたベルレF(t)  F (t + Δt )  っているためそのまま することはできませんをするときは
                    Δt
v ′(t + Δt) = v(t) + ---F (t)
                    2m
(12)

v(t + Δt) = v ′(t + Δt ) + Δt-F (t + Δt)
                         2m
(13)

つにけています

ベルレいたときのプログラムによりm =  1  として以下のような れになります

     v′(t + Δt ) = v(t) + Δ2tF (t)      速 度 の 更V [1]
   r(t + Δt ) = r(t) + Δtv ′(t + Δt)    位 置 の 更V
                  ↓
   r (t + Δt )か ら F (t + Δt )を 計Z    力 の 計Z
                  ↓
v(t + Δt ) = v′(t + Δt ) + Δ2tF (t + Δt ) 速 度 の 更V[2]

この流れをC言語風に書くと、位置をx、速度をv、力をfとして

v+=f*dt*0.5;  速度の更新[1]
x+=v*dt; 位置の更新
xf 力の計算
v+=f*dt*0.5;速度の更新[2]

となります。

ベルレ法のサンプルプログラム

サンプルとして下運ベルレいてみますt = 0  x =  0,v = 0,f = 9.8  t > 0  f = 9.8
    1
x = -f t2
    2
(14)

v = ft
(15)

けますこのベルレめたx, v  するプログラム以下のようになり ます

#include <stdio.h>
#include <math.h>
main()
{
  int i = 0;
  double x,v,f;
  double t,dt;
  dt = 0.1;
  x = 0.0;
  v = 0.0;
  f = 9.8;
  for(i = 1; i < 30; i++){
    t = i*dt;
    v += f*dt*0.5;
    x += v*dt;
    f = 9.8;
    v += f*dt*0.5;
    printf("%f %8.4f %8.4f %8.4f %8.4f %8.4f\n",t,x,t*t*4.9,v,t*9.8,f);
  }
}

このプログラムの実行結果は以下のようになり、ベルレ法で計算した値と理論値 が一致していることがわかります。

virlet0.c virlet1.c
  時間     ベルレx  理論値x  ベルレv  理論値v    力
0.100000   0.0490   0.0490   0.9800   0.9800   9.8000
0.200000   0.1960   0.1960   1.9600   1.9600   9.8000
0.300000   0.4410   0.4410   2.9400   2.9400   9.8000
0.400000   0.7840   0.7840   3.9200   3.9200   9.8000

中略

2.600000  33.1240  33.1240  25.4800  25.4800   9.8000
2.700000  35.7210  35.7210  26.4600  26.4600   9.8000
2.800000  38.4160  38.4160  27.4400  27.4400   9.8000
2.900000  41.2090  41.2090  28.4200  28.4200   9.8000

演習

粒子の初期位置を-1.0、初期速度を0.0、初期状態での力を1.0とし、 粒子に働く力がcos(t)で変化する場合の粒子運動を、 ベルレ法を使って数値計算し理論値と比較し、 そのグラフを作成せよ。 また、ベルレ法で使用する時間刻み幅を変更したとき、 計算値と理論値の差はどのようになるか調べよ。