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

问：

这个想法是建立一个欧拉积分方法，然后应用耦合常微分方程来描述水星绕太阳公转时的位置和速度。请参阅此处的耦合常微分方程。

我需要帮助的问题是，在绘制数据后，尽管更改了参数，但我只得到了圆形轨道 - 并且轨道周期是错误的，在天文单位中应该是 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()

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