Navigation

長時間の計算

　長時間の計算をしてみよう． \(|v_0|=1.0\) の初期条件であれば円軌道を回り続けるはずである．

長時間の計算を行う前に

tmax = 300000.0d0，n=3000000（時間刻みは0.1）まで計算を行ってみよう．

もちろん，先のプログラムで上記のようにパラメータを変更すれば計算をすることはできる．しかし，いくつか注意がある． これまでのプログラムでは毎ステップごとに時間，位置，速度を出力している．それをそのまま3000000ステップ分を出力すると300Mbを超える容量のファイルになってしまう． また実はファイルへの出力は意外と計算時間を要する．私の環境では，この計算に120秒程度の時間がかかった． これくらいなら待てると思うかもしれないが，たかが二体の計算でこの時間は粒子の数が増えた時のことを考えると恐ろしい． ここで興味があるのは，長時間の間，円軌道が維持されるのかどうかである．毎ステップごとに出力しなくてもよいだろう． ということで，出力部分に以下のような変更を行う（前ページの102~107行目あたり）．

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
integer                     i,dout
integer,parameter::         n=3000000
integer,parameter::         noutmax=1000
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

if(n .gt. noutmax) then
    dout = int(n/noutmax)
else
    dout = 1
end if

open(unit=10,file='outputRK.dat',status='replace')
write(10,*)'t_rk, x_rk, y_rk, E_rk'
do i = 1,n
    if(mod(i,dout) .eq. 0) then
        write(10,*)t(i),',',x(i),',',y(i),',',E(i)
    end if
end do
close(10)

この変更では，データの最大出力数（noutmax=1000）を設定する． 全計算ステップがこの数を超える場合は，この出力数内に収まるように出力間隔を調整する．

if(n .gt. noutmax) then
    dout = int(n/noutmax)
else
    dout = 1
end if

まずは，とうとう条件分岐if文の登場． 眺めていればわかると思うが，ifの後の括弧内の条件が真の時， つまりnがnoutmaxよりも大きければ（.gt.はgreater thanの意，n > noutmaxと書いてもよい，比較演算子という） その下にかかれた処理（dout = int(n/noutmax)）を実行せよ． そうではない（偽）ならば，else以下の処理（dout = 1）を実行せよ．である． 最後にend ifを書いて条件分岐の処理を終える．

n=300000はnoutmax=1000よりも大きいので，dout=int(300000/1000)=300を計算する．

Note

「if (~) then」の次に「else if (*) then」を加えて，条件を重ねていくことができる．

Note

比較演算子には，「.lt., <（より小さい）」　「.le., <=（以下）」「.eq.,==（等しい）」「.ne.,/=（等しくない）」「.ge.,>=（以上）」がある．

Note

論理演算子には「.and.」「.or.」などがある．条件の結合などが行える．「 \(0 \lt n \leq 10\)」は「if(n .gt. 0.0d0 .and. n .le. 10.0d0) then ...」と書ける．

Note

int()は実数を整数に切り捨てで変換する．ex. int(10.0/3.0)=int(3.33)=3

do i = 1,n
    if(mod(i,dout) .eq. 0) then
        write(10,*)t(i),',',x(i),',',y(i),',',E(i)
    end if
end do

次に出力部は，doループの中にif文をいれて「mod(i,dout) .eq. 0」が真の時だけ出力をしないさいとする．

mod(i,dout)は余り関数で，iをdoutで割った時の余りを返す． 今，dout=300なので，iが300の倍数のときだけmod(i,dout)はゼロになる．この時if文の条件が真になり出力の処理が実行される． つまり，iが300の倍数のときだけ出力される．

このようにすることで，どんなにステップ数を増やしても出力数は常にnoutmaxを維持する．

さて計算

tmax = 300000.0d0，n=3000000で計算した軌道である．

X vs Y　of celestial body 300k sec

Note

データの出力間隔が広い時に，折れ線をつかってプロットすると軌道ではないところを線で結んでしまうことがある．そのためここでは点でプロットした．

軌道に厚みが出てしまっている．半径1の円を周回するはずが，徐々に半径が小さくなっている． これは明らかにおかしい．

時間刻みを小さくしてみよう．

X vs Y　of celestial body 300k sec

赤い点がtmax = 300000.0d0，n=30000000(時間刻み0.01)で計算した結果である． 半径1の周回軌道を維持しているようである．時間刻み0.1の場合には，長時間の計算による計算誤差の蓄積によって不自然な軌道が現れていたのだ． ただし当然ながら計算ステップを増やした分だけ計算時間は増える．約10秒もかかった．

万有引力の作用する（太陽を中心に固定した）系の全エネルギーEは，

\[E = \frac{1}{2}m_e v^2 - G\frac{m_s m_e}{r}\]
\[E^* = E/{\left(\frac{Gm_e m_s}{a}\right)}=\frac{1}{2}{v^*}^2 - \frac{1}{r^*}\]

である．以下のような計算を追加して時間変化をプロットした．

E = 0.5d0*(u**2+v**2) - 1.0d0/sqrt(x**2+y**2)
Total energy for long time

時間刻み0.1の場合，時間とともにエネルギーが減少してしまっている．それに比べれば時間刻み0.01は一定を維持しているように見える．

Note

ルンゲ－クッタ法は長時間の計算を行うとエネルギーが減少してしまうことが知られている．〔 ハミルトン力学系のためのシンプレクティック数値積分法，吉田春夫 〕

さらに長時間

時間刻み0.01でも，さらに長時間計算すればいずれはエネルギーが減少するのだろう，と予測して， tmax = 1000000.0d0，n=10000000(時間刻み0.01)を行おうとすると，コンパイルのところでエラーがでる．

$ gfortran.exe celestial_RK.f90
/tmp/ccwZ2sLL.s: Assembler messages:
/tmp/ccwZ2sLL.s:1058: Error: value of 4800005200 too large for field of 4 bytes at 5136
/tmp/ccwZ2sLL.s:1176: Error: value of 4800005793 too large for field of 4 bytes at 5729

数を変えてみるとtmax = 380000.0d0，n=3800000まではコンパイルできるが， tmax = 390000.0d0，n=3900000だとエラーがでる．ただしこの場合はエラーメッセージが異なる．

*/libgfortran.dll.a(d001123.o):(.text+0x2): 再配置がオーバーフローしないように切り詰められました: *
良くわからないが，数が大きくなるとエラーだし，メモリの確保に問題があるのかもしれない．
これまでのコードでは，xyuvEの全時間分データを配列で記憶している．ステップ数（n）が増えるとメモリの確保が厳しそうである．
3900000*7(配列の数)*8(倍精度実数の情報量,byte)=218Mbyteだが，このあたりに何かあるのか？
さらに，n=3800000では，コンパイルができても実行はできない．実行までできるのはn=3500000あたりまでのようである．
明確な原因はわからないが，メモリの使い過ぎは明らかなので，配列の使用を減らす努力をしよう．

コードをよく眺めていると，配列を使う必要がなかったような気がしてくる． これまでのコードでは，位置と速度などの情報を配列に記憶した後，それを再び呼び出してエネルギーを計算したりファイルに出力したりしていた． 最初のdoループですべて済むのではないか？doループの中でエネルギーを計算し，ファイルに出力すれば良さそうだ．

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
!+++++++++++++++++++++++++++++++++++++++++++++++
program celestial

use globals
use subfunc

implicit none

integer                     i,dout
integer,parameter::         n=10**9
integer,parameter::         noutmax=10000
double precision            tmax
double precision            t,x,y,u,v
double precision            E,E0,dE
double precision            kx1,kx2,kx3,kx4,ku1,ku2,ku3,ku4
double precision            ky1,ky2,ky3,ky4,kv1,kv2,kv3,kv4
double precision            start_time,finish_time

!-----CPU time
call CPU_TIME(start_time)

x = 1.0d0
y = 0.0d0
u = 0.0d0
v = 1.0d0
t = 0.0d0
tmax = 10.0d6

delta_t = (tmax-t)/dble(n)

open(unit=10,file='outputRK.dat',status='replace')
write(10,*)'t_rk, x_rk, y_rk, E_rk, dE_rk'

if(n .gt. noutmax) then
    dout = int(n/noutmax)
else
    dout = 1
end if

do i = 1, n-1

    E = 0.5d0*(u**2+v**2) - 1.0d0/dsqrt(x**2+y**2)
    if(i .eq. 1) then
        E0 = E
    end if
    dE = abs((E - E0)/E0)

    if(mod(i-1,dout) .eq. 0) then
        write(10,*)t,',',x,',',y,',',E,',',dE
    end if

    kx1 = fx(x,y,u)*delta_t
    ky1 = fy(x,y,v)*delta_t
    ku1 = fu(x,y,u)*delta_t
    kv1 = fv(x,y,v)*delta_t
    kx2 = fx(x+0.5d0*kx1,y+0.5d0*ky1,u+0.5d0*ku1)*delta_t
    ky2 = fy(x+0.5d0*kx1,y+0.5d0*ky1,v+0.5d0*kv1)*delta_t
    ku2 = fu(x+0.5d0*kx1,y+0.5d0*ky1,u+0.5d0*ku1)*delta_t
    kv2 = fv(x+0.5d0*kx1,y+0.5d0*ky1,v+0.5d0*kv1)*delta_t
    kx3 = fx(x+0.5d0*kx2,y+0.5d0*ky2,u+0.5d0*ku2)*delta_t
    ky3 = fy(x+0.5d0*kx2,y+0.5d0*ky2,v+0.5d0*kv2)*delta_t
    ku3 = fu(x+0.5d0*kx2,y+0.5d0*ky2,u+0.5d0*ku2)*delta_t
    kv3 = fv(x+0.5d0*kx2,y+0.5d0*ky2,v+0.5d0*kv2)*delta_t
    kx4 = fx(x+kx3,y+ky3,u+ku3)*delta_t
    ky4 = fy(x+kx3,y+ky3,v+kv3)*delta_t
    ku4 = fu(x+kx3,y+ky3,u+ku3)*delta_t
    kv4 = fv(x+kx3,y+ky3,v+kv3)*delta_t

    x = x + (kx1+2.0d0*kx2+2.0d0*kx3+kx4)/6.0d0
    u = u + (ku1+2.0d0*ku2+2.0d0*ku3+ku4)/6.0d0
    y = y + (ky1+2.0d0*ky2+2.0d0*ky3+ky4)/6.0d0
    v = v + (kv1+2.0d0*kv2+2.0d0*kv3+kv4)/6.0d0
    t = t + delta_t

end do

close(10)

!-----Output CPU time
call CPU_TIME(finish_time)
write(*,*) 'CPU time:',finish_time - start_time

end program celestial
!+++++++++++++++++++++++++++++++++++++++++++++++
ということで，配列を一つも使わずに書くことができる．
まず更新前の位置と速度の情報をつかってエネルギーの計算，ファイル出力を行うようにした（初期条件での出力ができるように）．また計算時間がかかるようになってきたので，およその計算時間を表示するようにした．
call CPU_TIME(start_time)

!main process!

call CPU_TIME(finish_time)
write(*,*) 'CPU time:',finish_time - start_time

Note

CPU_TIME()は処理系の時間を返してくれる組み込みサブルーチンである．上記のようにプログラムのはじめとおわりで2回呼び出し時刻を記録し，その差を計算することで計算にかかったおおよその時間を知ることができる．

さてtmax = 10000000.0d0，n=1000000000（時間刻みは0.01）の計算が無事に実行できた．軌道と全エネルギーは，

torajectory of celestial body for 10M
Total energy for 10M
軌道は見事に円を維持している．すばらしい．
ただし全エネルギーは，やはり一方的に減少を続けている．これは何とかならないだろうか？
ところで調子に乗ってさらに10倍長く，tmax = 100000000.0d0，n=10000000000の計算を試みたが，またもやコンパイルが通らない．
integer,parameter::         n=10**10
                                1
エラー: Arithmetic overflow at (1)

整数型（4byte）の扱える数の範囲は-2,147,483,648 ～ 2,147,483,647（≒2×10 9 ）である． これを超えてしまったためである．

Table Of Contents

Previous topic

万有引力

Next topic

例題：ファン・デル・ポール振動子

This Page

Navigation

© Copyright 2014, Kentaro TANAKA. Created using Sphinx 1.2.1.