数値積分の実行

ベルレ法

分子動力学シミュレーションでは粒子を運動方程式に従って動かします。実際の粒子の位置や速度を求めるためには運動方程式から作った差分式を利用します。粒子が $N$ 個あるとき、運動方程式は

$\begin{equation*} \frac{d^2\boldsymbol{r}_i}{dt^2} = \frac{\boldsymbol{F}_i}{m_i}\;\;\;\;(i=1,2,\ldots,N) \end{equation*}$ (1)

となり、速度は

$\begin{equation*} \frac{d\boldsymbol{r}_i}{dt} = \boldsymbol{v}_i\;\;\;\;(i=1,2,\ldots,N) \end{equation*}$ (2)

と書けます。これらの式は $N$ 個ありますが、全て同じ形なので以降は添え字の $i$ を省略します。

時刻 $t+\Delta t$ での位置 $\boldsymbol{r}(t+\Delta t)$ をテーラー展開すると

$\begin{equation*} \boldsymbol{r}(t+\Delta t) = \boldsymbol{r}(t) + \Delta t\frac{d\boldsymbol{r}(t)}{dt} + \frac{\Delta t^2}{2}\frac{d^2\boldsymbol{r}(t)}{dt^2} + \cdots \end{equation*}$ (3)

この式と式(1)、式(2)を使うと

$\begin{equation*} \boldsymbol{r}(t+\Delta t) = \boldsymbol{r}(t) + \Delta t\boldsymbol{v}(t)+ \frac{\Delta t^2}{2}\frac{\boldsymbol{F}(t)}{m} + O((\Delta t)^3) \end{equation*}$ (4)

となり、同様に $\boldsymbol{r}(t-\Delta t)$ は

$\begin{equation*} \boldsymbol{r}(t-\Delta t) = \boldsymbol{r}(t) - \Delta t\boldsymbol{v}(t)+ \frac{\Delta t^2}{2}\frac{\boldsymbol{F}(t)}{m} + O((\Delta t)^3) \end{equation*}$ (5)

となります。式(4)と式(5)の和と差を取ると

$\begin{equation*} \boldsymbol{r}(t+\Delta t) + \boldsymbol{r}(t-\Delta t) = 2\boldsymbol{r}(t) + \Delta t^2\frac{\boldsymbol{F}(t)}{m} + O((\Delta t)^4) \end{equation*}$ (6)

$\begin{equation*} \boldsymbol{r}(t+\Delta t) - \boldsymbol{r}(t-\Delta t) = 2\Delta t\boldsymbol{v}(t) + O((\Delta t)^3) \end{equation*}$ (7)

となります。この式では $(\Delta t)^3$ 以上の項は無視しています。この2式から

$\begin{equation*} \boldsymbol{r}(t+\Delta t) = 2\boldsymbol{r}(t) - \boldsymbol{r}(t-\Delta t)+ \Delta t^2\frac{\boldsymbol{F}(t)}{m} \end{equation*}$ (8)

$\begin{equation*} \boldsymbol{v}(t) = \frac{1}{2\Delta t}\left\{\boldsymbol{r}(t+\Delta t) - \boldsymbol{r}(t-\Delta t)\right\} \end{equation*}$ (9)

が得られます。この式はベルレの式と呼ばれています。

ベルレの式では得られる位置と速度の時間がずれているので、式 (8)と式(9)を以下のように変形します。

$\begin{equation*} \boldsymbol{r}(t+\Delta t) = \boldsymbol{r}(t) + \Delta t\boldsymbol{v}(t) + \frac{\Delta t^2}{2}\frac{\boldsymbol{F}(t)}{m} \end{equation*}$ (10)

$\begin{equation*} \boldsymbol{v}(t +\Delta t) = \boldsymbol{v}(t)+\frac{\Delta t}{2m}\left\{\boldsymbol{F}(t)+\boldsymbol{F}(t+\Delta t)\right\} \end{equation*}$ (11)

この式は速度ベルレの式と呼ばれています。

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

上で得られた速度ベルレの式は、速度の計算に $\boldsymbol{F}(t)$ と $\boldsymbol{F}(t+\Delta t)$ が入っているためそのまま計算することはできません。実際に計算をするときは、式を

$\begin{equation*} \boldsymbol{v}(t+\textstyle\frac{\Delta t}{2}\displaystyle) = \boldsymbol{v}(t)+\frac{\Delta t}{2m}\boldsymbol{F}(t) \end{equation*}$ (12)

と

$\begin{equation*} \boldsymbol{v}(t+\Delta t) = \boldsymbol{v}(t+\textstyle\frac{\Delta t}{2}\displaystyle)+\frac{\Delta t}{2m}\boldsymbol{F}(t+\Delta t) \end{equation*}$ (13)

の2つに分けて、力の計算の前後で計算を行います。

ベルレ法を用いたときのプログラムを作る場合は、無次元化により $m=1$ として、以下のような流れになります。

$\boldsymbol{v}(t+\textstyle\frac{\Delta t}{2}\displaystyle) = \boldsymbol{v}(t)+\frac{\Delta t}{2}\boldsymbol{F}(t)$	速度の更新[1], 式(12)
$\boldsymbol{r}(t+\Delta t) = \boldsymbol{r}(t)+\Delta t\boldsymbol{v}(t+\textstyle\frac{\Delta t}{2}\displaystyle)$	位置の更新, 式(10)
$\boldsymbol{r}(t+\Delta t)\;\;\;\rightarrow\;\;\;\boldsymbol{F}(t+\Delta t)$	力の計算
$\boldsymbol{v}(t+\Delta t) = \boldsymbol{v}(t+\textstyle\frac{\Delta t}{2}\displaystyle)+\frac{\Delta t}{2}\boldsymbol{F}(t+\Delta t)$	速度の更新[2], 式(13)

となります。

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

v	+=	fdt0.5;	速度の更新[1], 式(12)
x	+=	v*dt;	位置の更新, 式(10)
x	→	f	力の計算
v	+=	fdt0.5;	速度の更新[2], 式(13)

となります。

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

サンプルとして落下運動をベルレ法で解いてみます。 $t=0$ で、 $x=0, v=0, f=9.8$ 、 $t>0$ で一定の力 $f=9.8$ が働く場合の物体の運動は

$\begin{equation*} x = \frac{1}{2}ft^2 \end{equation*}$ (14)

$\begin{equation*} v = ft \end{equation*}$ (15)

と書けます。この結果とベルレ法で求めた $x, v$ の結果を比較するプログラムは以下のようになります。

C言語

#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);
  }
}

Python

dt = 0.1
x = 0.0
v = 0.0
f = 9.8
for i in range(1,30):
  t = i*dt
  v += f*dt*0.5
  x += v*dt
  f = 9.8
  v += f*dt*0.5
  print(f'{t:f} {x:8.4f} {t*t*4.9:8.4f} {v:8.4f} {t*9.8:8.4f} {f:8.4f}')

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

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)で変化する場合の粒子運動を、ベルレ法を使って数値計算し理論値と比較し、そのグラフを作成せよ。また、ベルレ法で使用する時間刻み幅を変更したとき、計算値と理論値の差はどのようになるか調べよ。

このページではプログラミング言語として C 言語と Python を使用していますが、 Python で作ったプログラムの動作速度は C 言語に比べて遅いので、 MD シミュレーションのプログラムは C 言語のものを使うことを想定しています。

Python のプログラムは C 言語のプログラムを理解することを目的に作られており、なるべく C 言語に近い書き方になっています。そのため、 Python のプログラムでは、効率の悪い書き方になっている箇所もありますので注意してください。

Koishi's Page