我正在尝试使用 Fortran 95 来实现 N 体问题。在将其推广到 N 体之前，我首先尝试计算单个轨道体（下面我的代码中的体 2）围绕一个固定体（体 1 在下面的代码）。问题是输出不正确，但我的方程式找不到错误！这就是我想要实现的：

$V_{2,x} \leftarrow V_{2,x} - \frac{ G \times M_1 \times (x_2-x_1)}{r^3} \times dt$

$V_{2,y} \leftarrow V_{2,y} - \frac{ G \times M_1 \times (y_2-y_1)}{r^3} \times dt$

$x_2 \leftarrow x_2 + V_{2,x} \times dt$

$y_2 \leftarrow y_2 + V_{2,y} \times dt$

这是我的代码：

program gravitation
implicit none

integer :: k, total
real, dimension(2) :: M, x, y, Vx, Vy
real, parameter :: G = 6.6738E-11, pi = 3.1415926
real :: time, dt, r

dt = 1000. ! dt = 1000 seconds
total = 5000 ! total number of points to be output 

M(1)=1.9891E+09 ! Sun's mass
x(1)=0.
y(1)=0.
Vx(1)=0.
Vy(1)=0.

M(2)=5.9736E+03 ! Earth's mass
x(2)=1.4960E+11 ! Earth's average orbit as initial point
y(2)=0.
Vx(2)=0.
Vy(2)=2.9786E+04 ! Earth's average orbital speed

open(unit=101,file="a.txt")
write(101,*) "time, x, y"
write(101,"(E10.4,2(A1,1x,E10.4))") 0., ",", x(2), ",", y(2)

do k=1,total
  r = sqrt( (x(2)-x(1))**2 + (y(2)-y(1))**2 )
  Vx(2) = Vx(2) - G * M(1) * dt * (x(2)-x(1)) / r**3 ! i.e., Vx2 = Vx2 - G.M.dt.(x2-x1)/r^3
  Vy(2) = Vy(2) - G * M(1) * dt * (y(2)-y(1)) / r**3
    x(2) = x(2) + Vx(2) * dt
    y(2) = y(2) + Vy(2) * dt
  write(101,"(E10.4,2(A1,1x,E10.4))") k*dt, ",", x(2), ",", y(2)
enddo

end program gravitation

输出文件显示第二个body的x位置完全没有变化（卡在初始值0.1496E+12），而y位置发生了变化。

有谁看到这里有什么问题？非常感谢。