Navigation

多体系

　金星や火星あるいは月を加えた系を計算したくなる． 地球と太陽の二つしかないときは二体問題だが，月が増えれば三体問題．三体以上は多体と言い，解析的に軌道を得ることができないことが知られている．

数が増えると

これまでのコードの書き方ではいくつかの問題が生じる． 一つは計算の主要部分，

fu_temp = fu(x,y)
fv_temp = fv(x,y)
x = x + fx(u,v)*delta_t + 0.5d0*fu_temp*delta_t**2
y = y + fy(u,v)*delta_t + 0.5d0*fv_temp*delta_t**2
u = u + 0.5d0*(fu_temp + fu(x,y))*delta_t
v = v + 0.5d0*(fv_temp + fv(x,y))*delta_t

これまでの二体問題では，（太陽を中心に固定していたので） 地球の位置と速度の情報（x,y,u,v）だけを更新していればよかったが， これにもう一つ天体が加わると，その分だけ変数が（2次元であれば）４つ増えることになる． さらにもう一つ増やすと．．．と天体を増やすたびに変数がどんどん増えていく． おまけに変数ごとに上記のような計算式を書いて，となると非常に面倒である．

次に，大切な相互作用の計算部分でも，

r = sqrt(x**2+y**2)
xx = -y / r**3

これまでは，太陽との万有引力だけを計算していればよかったが， 新たに天体を加えたら，その天体との相互作用も計算しなくてはいけない．

そして， 運動方程式（再掲）も，

(1)\[\frac{d^2 \boldsymbol{x}_i}{dt^2} = \sum_{j \neq i}Gm_j\frac{\boldsymbol{x}_j- \boldsymbol{x}_i}{|\boldsymbol{x}_j-\boldsymbol{x}_i|^3}\]

太陽と地球だけだった場合は，無次元化により方程式の係数をすっきりさせることができたが， それぞれ異なる質量をもっている天体が存在する場合はそうはいかない． 無次元化はするが，それぞれの質量に関する情報も必要である．

ということで，思い切って，どっと変更．

  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
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
!+++++++++++++++++++++++++++++++++++++++++++++++
module globals

double precision,parameter::    pi = acos(-1.0d0)
double precision,parameter::    G = 6.67384d-11

integer,parameter::             ndim = 2
integer,parameter::             nbody = 5

end module globals
!+++++++++++++++++++++++++++++++++++++++++++++++
module var_trajectory

use globals
implicit none

double precision                coord(ndim,nbody)
double precision                veloc(ndim,nbody)

end module var_trajectory
!+++++++++++++++++++++++++++++++++++++++++++++++
module param_planet

use globals
implicit none

double precision                mass(nbody)

end module param_planet
!+++++++++++++++++++++++++++++++++++++++++++++++
module cal_ene

use globals
use param_planet
implicit none

contains
    subroutine cal_kin(v,e)
        double precision,intent(in)::     v(ndim,nbody)
        double precision,intent(out)::    e
        integer     i,j
        double precision    v2

        e = 0.0d0
        do i=1,nbody
            v2 = 0.0d0
            do j=1,ndim
                v2 = v2 + v(j,i)**2
            end do
            e = e + 0.5d0*mass(i)*v2
        end do
    end subroutine cal_kin
    !++++++++++++++++++++++++
    subroutine cal_pot(c,e)
        double precision,intent(in)::   c(ndim,nbody)
        double precision,intent(out)::    e
        integer     i,j,k
        double precision    rij(ndim)
        double precision    rij2,rij_abs

        e = 0.0d0
        do i=1,nbody-1
            do j=i+1,nbody
                rij(:) = c(:,j)-c(:,i)
                rij2 = 0.0d0
                do k=1,ndim
                    rij2 = rij2 + rij(k)**2
                end do
                rij_abs = sqrt(rij2)
                e = e + mass(i)*mass(j)/rij_abs
            end do
        end do
    end subroutine cal_pot

end module cal_ene
!+++++++++++++++++++++++++++++++++++++++++++++++
module cal_force

use globals
use param_planet
implicit none

contains
    subroutine force(c,a)
        double precision,intent(in)::       c(ndim,nbody)
        double precision,intent(out)::      a(ndim,nbody)
        integer     i,j,k
        double precision    rij(ndim)
        double precision    rij2,rij_abs,inv_rij_abs
        double precision    ftmp
        double precision    fij(ndim),f(ndim,nbody)

        f = 0.0d0
        do i = 1,nbody-1
            do j = i+1,nbody
                rij(:) = c(:,j)-c(:,i)
                rij2 = 0.0d0
                do k=1,ndim
                    rij2 = rij2 + rij(k)**2
                end do
                rij_abs = sqrt(rij2)
                inv_rij_abs = 1.0d0/rij_abs
                ftmp = mass(i)*mass(j)*(inv_rij_abs**3)
                fij(:) = ftmp*rij(:)
                f(:,i) = f(:,i) + fij(:)
                f(:,j) = f(:,j) - fij(:)
            end do
        end do
        do i = 1,nbody
            a(:,i) = f(:,i) / mass(i)
        end do
    end subroutine force

end module cal_force
!+++++++++++++++++++++++++++++++++++++++++++++++
program celestial

use globals
use var_trajectory
use param_planet
use cal_ene
use cal_force

implicit none

integer                     i,dout
integer,parameter::         n=10**6
integer,parameter::         noutmax=10000
integer,parameter::         nlive=10
double precision            tmax
double precision            t,dt,dt2
double precision            cnor,vnor,mnor,tnor
double precision            E,E0,Ev,Ep,dE,Eold,sum_dE
double precision            start_time,finish_time
double precision            accel(ndim,nbody)
integer                     nliveout

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

!----- mass
mass(1) = 1.989d30          ! mass of sun[kg]
mass(2) = 5.972d24          ! mass of earth[kg]
mass(3) = 4.869d24          ! venus
mass(4) = 6.419d23          ! mars
mass(5) = 7.348d22          ! mass of moon[kg]

!----- initial conditions
coord(1,1) = 0.0d0
coord(2,1) = 0.0d0
coord(1,2) = 1.496d11       !average radius of revolution of earth[m]
coord(2,2) = 0.0d0
coord(1,3) = 1.082d11
coord(2,3) = 0.0d0
coord(1,4) = 2.279d11
coord(2,4) = 0.0d0
coord(1,5) = 3.844d8        !average radius of revolution of moon[m]
coord(1,5) = coord(1,5) + coord(1,2)
coord(2,5) = 0.0d0

veloc(1,1) = 0.0d0
veloc(2,1) = 0.0d0
veloc(1,2) = 0.0d0
veloc(2,2) = 2.978d4        !average velocity of earth[m/s]
veloc(1,3) = 0.0d0
veloc(2,3) = 3.502d4
veloc(1,4) = 0.0d0
veloc(2,4) = 2.413d4
veloc(1,5) = 0.0d0
veloc(2,5) = 1.022d3        !average velocity of moon[m/s]
veloc(2,5) = veloc(2,5) + veloc(2,2)


!----- time
t = 0.0d0
tmax = 3.2d8              !365day*24hour*60min*60sec = 31536000sec = 3.2d7

!----- normalizse
cnor = sqrt((coord(1,2)-coord(1,1))**2 + (coord(2,2)-coord(2,1))**2)
mnor = mass(1)
tnor = sqrt(cnor**3/(G*mnor))
vnor = cnor/tnor

coord = coord/cnor
veloc = veloc/vnor
mass = mass/mnor
t = t/tnor
tmax = tmax/tnor

!------ main
dt = (tmax-t)/dble(n)
dt2 = 0.5d0*dt

open(unit=10,file='outputNvV.dat',status='replace')
write(10,*)'t_NvV, E_NvV, dE_NvV, x1_NvV, y1_NvV, x2_NvV, y2_NvV, &
    &  x3_NvV, y3_NvV, x4_NvV, y4_NvV, x5_NvV, y5_NvV'

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

do i = 1, n-1

    if(mod(i,nliveout) .eq. 0) then
        write(*,*)'current step:',i,'/total step:',n
    end if

    if(i .eq. 1) then
        call cal_kin(veloc,Ev)
        call cal_pot(coord,Ep)
        E0 = Ev-Ep
        Eold = E0
    end if

    call cal_kin(veloc,Ev)
    call cal_pot(coord,Ep)
    E = Ev-Ep

    dE = abs((E - E0)/E0)
    sum_dE = sum_dE + abs(E-Eold)

    if(mod(i-1,dout) .eq. 0) then
        write(10,*)t,',',E,',',dE,',',coord(1,1),',',coord(2,1),',',coord(1,2),',',coord(2,2),',',coord(1,3) &
        & ,',',coord(2,3),',',coord(1,4),',',coord(2,4),',',coord(1,5),',',coord(2,5)
    end if

    if(i .eq. 1) then
        call force(coord,accel)
    end if

    veloc = veloc + dt2 * accel
    coord = coord + dt * veloc
    call force(coord,accel)
    veloc = veloc + dt2 * accel

    t = t + dt
    Eold = E

end do

close(10)

write(*,*) 'Time step:',dt
write(*,*) 'Average error of energy update:',sum_dE/dble(n)

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

end program celestial
!+++++++++++++++++++++++++++++++++++++++++++++++

貼り付けるとかなり長いコードになってしまったが，中身はそれほど変わっていない．

大切な主要部分は書き方を変更．

if(i .eq. 1) then
    call force(coord,accel)
end if

veloc = veloc + dt2 * accel
coord = coord + dt * veloc
call force(coord,accel)
veloc = veloc + dt2 * accel

以前の書き方よりもすっきりしているように見える． ただし，座標（coord），速度（veloc）と加速度（accel）の変数はすべて配列で記述している． モジュール「var_trajectory」の中で定義されているように，ndim=2,nbody=5，2×5の大きさをもつ配列である． coord(1,1)は1番目の天体のx座標，veloc(2,5)は5番目の天体のy方向速度を格納する．

スキームはこれまでと同じく速度ベルレ法である．見た目が少々違うが，

\[\begin{split}\dot x &= u, \dot u &= \Phi (x)\end{split}\]

に対して，

次のようなアルゴリズムであり，中身は変わっていない．

\[\begin{split}&u^* = v(t) + \Phi (x(t))\frac{\Delta t}{2}\\ &x(t+\Delta t) = x(t) + u^* \Delta t\\ &u(t+\Delta t) = u^* + \Phi (x(t+\Delta t))\frac{\Delta t}{2}\end{split}\]

まず，半歩進めた仮の速度v* を計算する． 最初のステップ（i=1）では，計算に必要なaccelの値がないので，この時のみif文内で計算している． その後は，1ステップ前でaccelの計算が済んでいる． この速度を使って位置の更新． この位置を使って速度の更新である．

加速度の計算 \(\Phi (x)\) は，サブルーチンを作って行うことにした． このためのサブルーチン「force」はモジュール「cal_force」のなかに（contain）含まれている．

module cal_force

use globals
use param_planet
implicit none

contains
    subroutine force(c,a)
        double precision,intent(in)::       c(ndim,nbody)
        double precision,intent(out)::      a(ndim,nbody)
        integer     i,j,k
        double precision    rij(ndim)
        double precision    rij2,rij_abs,inv_rij_abs
        double precision    ftmp
        double precision    fij(ndim),f(ndim,nbody)

        f = 0.0d0
        do i = 1,nbody-1
            do j = i+1,nbody
                rij(:) = c(:,j)-c(:,i)
                rij2 = 0.0d0
                do k=1,ndim
                    rij2 = rij2 + rij(k)**2
                end do
                rij_abs = sqrt(rij2)
                inv_rij_abs = 1.0d0/rij_abs
                ftmp = mass(i)*mass(j)*(inv_rij_abs**3)
                fij(:) = ftmp*rij(:)
                f(:,i) = f(:,i) + fij(:)
                f(:,j) = f(:,j) - fij(:)
            end do
        end do
        do i = 1,nbody
            a(:,i) = f(:,i) / mass(i)
        end do
    end subroutine force

end module cal_force

系内にあるすべての天体間の組合せに対して距離を計算し，万有引力を計算している．

f(:,i) = f(:,i) + fij(:)
f(:,j) = f(:,j) - fij(:)

天体jから天体iへの作用fijを計算したら，iに作用する力f(:,i)に加える，とともに， 天体iから天体jへの作用（-fij，反作用）としてjに作用する力f(:,j)にも加えている． その他の部分はよく見れば何を計算しているか分かるだろう．

Note

たとえばdoループを使って，二天体間の距離の二乗rij2を計算しているが，この時ループの直前で「rij2=0.0d0」と値をゼロにするのを忘れないように．これを忘れると，その他の天体間のための計算に用いたrij2の値を引き継いでしまう．

Note

これまでの関数（function）による書き方と似ている．関数は，単一の値を返すのに対して，サブルーチンは引数に値を入れて返す．この場合はc=coordの値を受け取りa=accelに結果を入れて返している．

Note

天体の質量，初期位置，初速度はWikipediaの情報を元に決めた．Wikipediaは理科年表の値などを載せているようである．初期位置はx軸上の平均公転半径に相当する場所，初期速度は平均軌道速度とした．

さて，以下が結果である．

X vs Y　of many-body

なんてことはない，それぞれの天体の円（にしか見えない）軌道が描かれる．

よく見ないと分からないが，例えば月と地球の軌道を見比べれば，月が地球の軌道の外側と内側を交互に移動しながら周回する様子も分かる．

あるいは，たとえば金星の速度を20%増しにして長時間の計算をしてみると，金星だけでなく火星の軌道も徐々にずれてくる．

X vs Y　of many-body
delta_E of many body sim

一応，エネルギーの初期値からの変動も載せる． 時間刻みは \(\Delta t = 0.637\) （赤）から \(\Delta t = 0.637^{-5}\) （黒）まで 速度ベルレ法なので，初期値からの変動は小さく，また刻みを小さくすればするほど小さくなる．

Table Of Contents

Previous topic

エネルギーの保存

Next topic

可視化（アニメ）

This Page

Navigation

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