Euler Integration to plot orbit of Mercury, only plotting circular orbits and incorrect period

The idea is to set up an Euler integration method and then apply coupled ODEs that describe the position and velocity of Mercury as it orbits around the Sun. See here coupled ODEs.

The problem that I need help with, is after plotting the data, I only get circular orbits despite changing the parameters - and the orbital period is wrong, should be 0.24 in astronomical units and it returns 0.17. This period changes with fewer orbits iterated, also.

Please could someone point me in the right direction? I suspect it is a problem with the Euler updater and/or x_s, y_s use.

Please note the plotting code is separate to the theory, they read/write to a file to communicate.

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()

I have tried to use the Euler method on its own and it seems to work. Then I checked the ODEs and the constants and they are correct as far as I can see. I could alter the code to make the period longer, but I would be entirely cheating.

Solution

Your algorithm is fine, except that, if you want to rederive the known orbital period, the initial velocity can't be a free parameter but must be set from the other orbital parameters. Specifically, the instantaneous norm of the velocity for a given Sun-planet distance r is

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

For your setup, vy0 = v, with

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

gives you a period of 0.24080... for any stable delta_t. For other random vy0, you'll see that the orbit turns elliptic.

Note aside, using the Euler method is not appropriate for orbital mechanics. The reason is because this method (along with things like Runge-Kutta) are not conservative; their dissipative nature makes the energy decay. What you want is a so-called symplectic method like Leap-frog.