Симуляция N-тела на алгоритме скачка Фортрана

Question

Я использую простой алгоритм «перепрыгнуть», я стремлюсь симулировать орбиты Земли Юпитера вокруг Солнца. Я не могу вывести их на орбиту, хотя и уверен, что математика верна. Похоже, что гравитация действует слишком еженедельно, и pl anet просто уплывает от Солнца, что интересно, если я отрегулирую ньютоновское ускорение из-за члена гравитации, умножив его на 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

Answer 1

Ваша самая фундаментальная проблема заключалась в том, что 1 AU = 1,5e11 м, а не 15e11. Тогда вы делали такие вещи, как сброс tcount каждую поездку через l oop. Установите его до начала основного l oop, а затем сбрасывайте только при выводе строки вывода. Оно должно быть обновлено как tcount=tcount+dt, и тогда вы, вероятно, захотите распечатать r(1,2) , r(2,2), r(1,3) , r(2,3), чтобы вы могли построить положения Юпитера и Земли. Также вам может понадобиться go на большее время, чтобы вы могли увидеть несколько полных орбит Земли, и, наконец, поставить тест в нижней части l oop, чтобы он выходил при t>tend. Сделав эти изменения, я получил вывод, который выглядел так: