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

Question

I have a set of two ODEs like so: $$ \ddot r = \dot \theta^2r \ - \frac{GM}{r^2}$$ $$ \dot r = -\frac {r \ddot \theta} {2 \dot \theta}$$

$\theta$ and $r$ are functions of time, 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 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 $\theta$ 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 achieveable.

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? I understand that this is very much a computer science problem but I thought people on the maths forum would be more likely to know.

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

Answer 1

I am not certain to well understand the question. Helping that the equations below will help. $$ \dot r = -\frac {r \ddot \theta} {2 \dot \theta}\quad\implies\quad \frac{2}{r}\frac{dr}{dt}=-\frac{\frac{d^2\theta}{dt^2}}{\frac{d\theta}{dt}}$$ $$\ln(r^2)=-\ln\left|\frac{d\theta}{dt}\right|+\text{constant}$$ $$\boxed{\frac{d\theta}{dt}=\frac{c_1}{r^2}}$$ This is already a simpler differential equation than the origonal one.

Without boundary condition one cannot say if the constant $c_1$ can be immediately determined or not. $$ \ddot r = \dot \theta^2r \ - \frac{GM}{r^2}\quad\implies\quad \frac{d^2r}{dt^2}=\left( \frac{d\theta}{dt}\right)^2 r- \frac{GM}{r^2}$$ $$\frac{d^2r}{dt^2}=\left( \frac{c_1}{r^2}\right)^2 r- \frac{GM}{r^2}$$ $$\frac{d^2r}{dt^2}=\frac{c_1^2}{r^3} - \frac{GM}{r^2}$$ $$2\frac{d^2r}{dt^2}\frac{dr}{dt}=2\frac{c_1^2}{r^3}\frac{dr}{dt} - \frac{2GM}{r^2}\frac{dr}{dt}$$ $$\left(\frac{dr}{dt}\right)^2= -\frac{c_1^2}{r^2} +\frac{2GM}{r} +\text{constant}$$ $$\boxed{\frac{dr}{dt}=\pm\sqrt{-\frac{c_1^2}{r^2} +\frac{2GM}{r}+c_2}}$$ This first order ODE contains only $r$. $$t=\pm\int \frac{dr}{\sqrt{-\frac{c_1^2}{r^2} +\frac{2GM}{r}+c_2}}$$ $$t=\pm\frac{1}{c_2}\left(\sqrt{-\frac{c_1^2}{r^2} +\frac{2GM}{r}+c_2}- \frac{GM}{\sqrt{c_2}}\tanh^{-1}\left(\frac{\frac{GM}{\sqrt{c_2}}+\sqrt{c_2}\:r }{\sqrt{-\frac{c_1^2}{r^2} +\frac{2GM}{r}+c_2}} \right) \right)+c_4$$ This is the solution $t$ as a function of $r$.

Without boundary condition one cannot say if these complicated equations could be simplified.

Answer 2

The aim to prepare a system for numerical integration is to isolate the highest order derivatives and then produce an explicit first-order system. This is straight-forward for the given system. $$ \ddot r = \dot\theta^2r-\frac{GM}{r^2}\\ \ddot\theta = -\frac{2\dot r\dot\theta}{r} $$

Note that you can integrate the second equation to the conservation of angular momentum or area speed, $r^2\dot\theta=L=const.$, which allows to eliminate the $\theta$ in the first equation and proceed on towards the Kepler laws.

Now introduce first derivatives for the stepping down of the second order derivatives (either impulse variables or just the speeds), so that in the end \begin{align} \dot r &= v\\ \dot \theta &= \omega\\ \dot v &= r\omega^2-\frac{GM}{r^2}\\ \dot\omega &= -\frac{2v\omega}{r}\\ \end{align}