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

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$$.

• I suspect there is no way to write this equation as a system of linear first order ODEs. Do you understand how to transform it into a system of first order ODEs? Do you have additional boundary conditions? The value of $r$ and $\theta$ on the boundary will not be enough for a second order problem. Perhaps for example you have some sort of periodicity condition since this is about an orbit? Dec 21, 2023 at 11:53
• I don't know how to transform this into a linear first order system which is why I was hoping for there to be a way to solve it without doing so. In terms of additional boundary conditions, I could probably estimate r' and theta' from the Ephemerides but it may not be that accurate Dec 21, 2023 at 11:58
• These are not PDEs, they are ODEs. Dec 21, 2023 at 17:38
• stackoverflow.com/questions/53813499/… , stackoverflow.com/questions/65084648/… are places you can start Dec 21, 2023 at 17:42
• math.stackexchange.com/questions/4373805/… Dec 21, 2023 at 17:43

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}

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.

• This is not the answer he was looking for. Dec 21, 2023 at 17:38
• Sorry I didn't understood your question. I thought that you were looking for the ODE $$\frac{d^2r}{dt^2}=\frac{c_1^2}{r^3} - \frac{GM}{r^2}$$ which contains only one variable $r$. In addition I gave the analytic solution which is useful to compare with the numerical solution in order to check the accuracy of the numerical results. Dec 22, 2023 at 10:35
• The next intermediary step is to set $r(t)=1/u(\theta(t))$, directly leading to $r=\frac{R}{1+e\cos\theta}$ and the Kepler equation. Dec 22, 2023 at 15:15