fortran - Fortran 蛙跳算法的 N 体仿真

Question

我正在使用一个简单的“跨越式”算法，我的目标是模拟地球的轨道，即木星围绕太阳的轨道。尽管相当确定数学是正确的，但我无法将它们送入轨道。似乎重力每周都在起作用，行星只是漂浮在远离太阳的地方，有趣的是，如果我通过将重力项引起的牛顿加速度乘以 rad2 来调整它，我发现该系统确实产生了相当稳定的轨道，但也有很多大半径。

program physim
  Implicit none
  integer :: i,j,n,day ! Integer variables
  doubleprecision :: G , r(1:3,1:10) , a(1:3, 1:10) , v(1:3, 1:10) , m(1:3), dt, Au, dr(1:3), 
rad2(1:3), t, tcount, tend, tout

! constants
day = 86400
tout = 10*day
tend = 20*day

Au = 15e11
n = 3
G = 6.67e-11
!n = 2
dt = 100
!sun
r(1,1) = 0.
r(2,1) = 0.
r(3,1) = 0.
v(1,1) = 0.
v(2,1) = 0.
v(3,1) = 0.
m(1) = 1.9898e30
!earth
r(1,2) = Au
r(2,2) = 0.
r(3,2) = 0.
v(1,2) = 0.
v(2,2) = 30000
v(3,2) = 0.
m(2) = 6e24
!jupiter
r(1,3) = 5.2*Au
r(2,3) = 0.
r(3,3) = 0.
v(1,3) = 0.
v(2,3) = 13070
v(3,3) = 0.
m(3) = 2e27


do
a = 0
tcount = 0
do i = 1, n
    do j = 1, n
        !calculating acceleration
    if (i==j)cycle
    dr(1:3) = r(1:3, j) - r(1:3, i)

    rad2 = dr(1)**2 + dr(2)**2 + dr(3)**2
    a(1:3, i) = a(1:3, i) + G*m(j)*dr(1:3)/(rad2*sqrt(rad2))


    end do
end do

do i = 1, n
    r(1:3 ,i) = r(1:3, i) + v(1:3, i)*dt
    v(1:3, i) = v(1:3, i) + a(1:3, i)*dt

end do
t = t + dt
tcount = t + dt
if(tcount>tout) then
!write(6,*) a(1,2)
!write(6,*) rad2


write(6,*) a(1,1) , a(2,1), a(3, 2)


end if

end do
end program

标签： fortran

Answer 1

您最基本的问题是 1 AU = 1.5e11 m，而不是 15e11。然后你正在做一些事情，比如重置tcount循环中的每一次旅行。在主循环开始之前设置它，然后仅在打印出一行输出时重置。它应该更新tcount=tcount+dt，然后你可能想要打印出来，r(1,2) , r(2,2), r(1,3) , r(2,3)这样你就可以绘制木星和地球的位置。另外，您可能应该多花些时间，这样您就可以看到几个完整的地球轨道，最后在循环的底部进行测试，以便它会在何时退出t>tend。进行这些更改后，我得到了如下所示的输出：