pythonphysics# 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.

- Python Jinja2 LaTeX Table
- Getting attributes of a class
- How can I print many significant figures in Python?
- How to allow list append() method to return the new list
- Calculate Last Friday of Month in Pandas
- Python type hint for Iterable[str] that isn't str
- How to iterate over a list in chunks
- How to exit the entire application from a Python thread?
- Running shell command and capturing the output
- How do I pass a variable by reference?
- Convert range(r) to list of strings of length 2 in python
- How can I get the start and end dates for each week?
- how to use send_message() in python-telegram-bot
- Python conditional replacement based on element type
- How can I count the number of items in an arbitrary iterable (such as a generator)?
- Find longest consecutive range of numbers in list
- Insert text in braces with asyncpg
- How does one put a link / url to the web-site's home page in Django?
- How to determine if a path is a subdirectory of another?
- Custom Keybindings for Ipython terminal
- FastAPI asynchronous background tasks blocks other requests?
- How to make sure that information from one file is duplicated into several text documents, without specific lines
- Installing a Python environment with Anaconda
- sklearn pipeline model predicting same results for all input
- Brew command not found after installing Anaconda Python
- How to get an XPath from selenium webelement or from lxml?
- Pipe PuTTY console to Python script
- How to align the axes of a figure in matplotlib?
- Persist ParentDocumentRetriever of langchain
- How to reset index in a pandas dataframe?