How to numerically approximate solutions to a set of two second order two-variable PDEs using a python script (boundary value problem)?

Question

I have a set of two PDEs like so (pseudo-code):

r'' = (s')^2 * r - G*M*r^-2
r' = -(r*s'')/(2*s')

s and r are functions of time, ' the derivative, G and M are constants.

They were obtained from the Euler-Lagrange equations for the moon moving around the earth in 2D polar coordinates. They may even be wrong; their correctness is not the point of this question.

I have tried to go about solving them using the scipy library but it requires me to first separate the equations into first order ODEs such that each equation contains either r or s but not both. I frankly don't know how to do that and whilst this may be possible for these two equations, for the more complex versions I plan on working with it is unlikely to be achievable.

So my question is: How can I go about solving them directly without any manipulation? Are there any libraries out there that will do that?

Also, this is a boundary value problem. I know values of r and s at the start and end of my given time interval (from the JPL ephemerides) but not r' or s'.

EDIT: After looking at the comments I put together some code with scipy:

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
from skyfield.api import load, Topos

G = 6.67430e-11
M = 5.97219e24

eph = load("de421.bsp")
observer = eph["earth"] + Topos(latitude_degrees=0, longitude_degrees=0)
moon = eph["moon"]

time1 = load.timescale().utc(2000, 1, 1, 0, 0, 0)
moon_position = observer.at(time1).observe(moon)
theta0, _, r0 = moon_position.radec()
time2 = load.timescale().utc(2000, 1, 1, 0, 0, 1)
moon_position = observer.at(time2).observe(moon)
theta1, _, r1 = moon_position.radec()

theta0 = float(theta0.radians)
thetadot = theta0 - float(theta1.radians)
r0 = float(r0.m)
rdot = r0 - float(r1.m)
print(r0, r1.m)
print(theta0, theta1.radians)
print(thetadot, rdot)


def equations(t, y):
    theta, r, thetadot, rdot = y
    dydt = [rdot, (thetadot**2 * r - G * M / r**2), thetadot, -r * (thetadot**2) / (2 * rdot),]
    return dydt


t_eval = np.linspace(0, 10000, 1000000)
sol = solve_ivp(equations, [0, 10000], [r0, theta0, rdot, thetadot])
t_plot = sol.t
y_plot = sol.y

#plt.plot(t_plot, y_plot[0], label="r(t)")
plt.plot(t_plot, y_plot[1], label="theta(t)")
plt.xlabel("Time")
plt.ylabel("Values")
plt.legend()
plt.show()

Unfortunately no matter how long or short my timescale is theta always becomes 0 at the end which can't be correct. Also the way I'm representing the system of equations in the equations function is wrong, how would I do it correctly?

Answer 1

For an attempt at finding an analytical solution (although I note that the question asks for a numerical approximation), you could try using sympy, in particular the dsolve_system function, e.g.,

from sympy import *
from sympy.solvers.ode.systems import dsolve_system

s, r, G, M, t = symbols("s r G M t")

sf = Function(s)(t)
rf = Function(r)(t)

sp = Derivative(sf, t)  # derivative
spp = Derivative(sf, (t, 2))  # second derivative

rp = -(r * spp) / (2 * sp)
rpp = sp**2 * r - G * M / r**2

eqs = [Eq(Derivative(rf, t), rp), Eq(Derivative(rf, (t, 2)), rpp)]

dsolve_system(eqa, t=t)

Unfortunately, it seems that sympy can't solve this system of equations analytically.

Answer 2

In your equations, r is radius and s is polar angle (I’ll use theta below). If you multiply your second equation by 2.r.sdot and rearrange you get

or

Hence you have a first invariant

h is the angular momentum (per unit mass) and it is constant because gravity is a central force, so imparts no torque and hence no change of angular momentum.

Substituting in your first equation you get

If you multiply through by rdot and rearrange you get

or

and hence you have a second invariant

E is the total energy (per unit mass) It will be negative for a bound orbit. It is constant because there are no external forces.

So, really, you have two first-order ODEs:

and

which you could actually solve sequentially.

Obviously you have the two unknown invariants, h and E. However, you can estimate initial guesses for these from approximate motion in a circle:

where T is the period of the orbit.

If you eliminate time t then you should get a polar curve r(theta), which represents an ellipse.

Answer 3

You need to decide on one order of state components and keep to it.

At the moment you unpack

theta, r, thetadot, rdot = y

but you assemble the derivatives in the order

r, rdot, theta, thetadot

This can not work.

The initial values then have a still different order

r0, theta0, rdot, thetadot