欧拉积分绘制水星轨道，仅绘制圆形轨道和不正确的周期

Question

这个想法是建立欧拉积分方法，然后应用耦合常微分方程来描述水星绕太阳运行时的位置和速度。看这里。

我需要帮助的问题是，在绘制数据后，尽管更改了参数，我只得到圆形轨道 - 并且轨道周期是错误的，应该是天文单位的 0.24，它返回 0.17。这个周期也随着迭代轨道的减少而变化。

请有人指出我正确的方向吗？我怀疑这是 Euler 更新程序和/或 x_s、y_s 使用的问题。

请注意，绘图代码与理论是分开的，它们读取/写入文件以进行通信。

import numpy as np
import matplotlib.pyplot as plt #not necessary as no active visulalisation

# Constants for the orbit of Mercury
e = 0.205630
aph = 0.466697
per = 0.3074995
a = 0.387098
n = 1000000
delta_t = 0.000001

#take number of orbits as user input, fewer than 10 for computational time's sake
num_orbits = int(input("Enter the number of orbits to simulate: "))
n = int(num_orbits * per / delta_t)

# Setting the arrays
x = np.zeros(n)
y = np.zeros(n)
vx = np.zeros(n)
vy = np.zeros(n)

# Initial position of Mercury
x0 = -a
y0 = 0
x[0] = x0
y[0] = y0

# Initial speed of Mercury
vx0 = 0
vy0 = 12 #given in script
vx[0] = vx0
vy[0] = vy0

# Position of the Sun
x_s = -e * aph
y_s = 0

# Other constants
M = 1
G = 4 * np.pi**2

#functions for given ODE in script
def funcx_dotdot(x, y):
    r = np.sqrt((x - x_s)**2 + (y - y_s)**2)
    return -G * M * (x - x_s) / (r**3)

def funcy_dotdot(x, y):
    r = np.sqrt((x - x_s)**2 + (y - y_s)**2)
    return -G * M * (y - y_s) / (r**3)

# Euler method to update positions and velocities (for any similar ODE)
def euler_update(x, y, vx, vy, delta_t):
    ax = funcx_dotdot(x, y)
    ay = funcy_dotdot(x, y)
    
    new_vx = vx + delta_t * ax
    new_vy = vy + delta_t * ay
    
    new_x = x + delta_t * new_vx
    new_y = y + delta_t * new_vy
    
    return new_x, new_y, new_vx, new_vy

# Initialize orbital_period
orbital_period = 0

for i in range(1, n):
    x[i], y[i], vx[i], vy[i] = euler_update(x[i-1], y[i-1], vx[i-1], vy[i-1], delta_t)

# Find the time step when Mercury crosses its initial x-position
initial_x = x[0]

#work out orbital period based on when mercury returns to (near) starting position
for i in range(1, n):
    if (x[i] > initial_x and x[i - 1] < initial_x) or (x[i] < initial_x and x[i - 1] > initial_x):
        orbital_period = i * delta_t
        break

print(f"Orbital period of Mercury: {orbital_period} astronomical units of time")

# Create an array to store time values
time = np.arange(0, n * delta_t, delta_t)

#save data in columns separated by tab: time, x, y
data = np.column_stack([time, x, y])
datafile_path = "/Users/danclough/Spyder/xy_orbit.txt"
np.savetxt(datafile_path , data, fmt=['%f','%f', '%f'])

import matplotlib.pyplot as plt
import numpy as np

data = np.loadtxt('xy_orbit.txt')

time = data[:, 0]
x = data[:, 1]
y = data[:, 2]

# Plot the orbit
plt.figure(figsize=(8, 6))
plt.plot(x, y, linestyle='-', marker='.')
plt.scatter((-0.205630*0.466697), 0, color='orange', marker='o', label='Sun')
plt.scatter(x[-1], y[-1], color='red', marker='o', label='Mercury')
plt.title('Orbit of Mercury around the Sun')
plt.xlabel('X Position')
plt.ylabel('Y Position')
plt.axis('equal')
plt.legend()
plt.grid(True)
plt.show()

我尝试过单独使用欧拉方法，它似乎有效。然后我检查了常微分方程和常数，据我所知它们是正确的。我可以改变代码来延长时间，但我完全是在作弊。

Answer 1

你的算法很好，只是，如果你想重新推导已知的轨道周期，初始速度不能是自由参数，而必须根据其他轨道参数来设置。具体来说，对于给定的太阳-行星距离 r，速度的瞬时范数是

v = sqrt(GM(2/r - 1/a))

对于您的设置，

vy0 = v

，以及

r = np.sqrt((x0-x_s)**2.0 + (y0-y_s)**2.0)

为任何稳定的

0.24080...

提供一段

delta_t

的时间。对于其他随机 vy0，您会看到轨道变成椭圆形。

请注意，使用欧拉方法并不适合轨道力学。原因是因为这种方法（以及诸如 Runge-Kutta 之类的方法）并不保守；它们的耗散性质使能量衰减。你想要的是所谓的 symplectic 方法，比如 Leap-frog。