Project on the Kepler Problem (template)

4.3. Project on the Kepler Problem (template)#

Brenton LeMesurier, lemesurierb@charleston.edu

Version of November 19, 2025

4.3.1. Introduction#

The objective of this project is to compare various methods for solving initial value problems for ordinary differential equations. The main test case will be the \textbf{Kepler Problem} of modeling planetary motion. This is a useful test case because it has the complications of being a system of several nonlinear equations, and yet exact solutions are known, and there are conserved quantities that can be used as accuracy checks.

The initial task described below is to solve one case by one simple method (Euler’s method) and present the results graphically. Once that is working, we will compare with other time step sizes, compute solutions to far longer times, and compare to the results of various other more accurate time discretization methods. Graduate students will have the opportunity to carry the investigation further, such as by reading about and exploring more advanced time-discretization methods, or considering systems with more planets.

4.3.2. Background on the Kepler Problem, modeling a planetary orbit#

The system of two second order nonlinear ODEs

\[ \ddot{x} = - \frac{x}{(x^2 + y^2)^{3/2}}, \quad \ddot{y} = - \frac{y}{(x^2 + y^2)^{3/2}}, \quad t_0 \leq t \leq t_{final} \]

gives a description of the motion of one body around another under gravitational attraction. (Aside: the full description in terms of two bodies moving in three dimensional space can be reduced to this by using center of mass coordinates, and the fact that in these coordinates, the motion in space stays in a certain plane.)

This system has two conserved quantities, energy $H$ and angular momentum $L$,

\[\begin{split} \begin{split} H(t) &= \frac12 (\dot{x}^2 + \dot{y}^2) - \frac{1}{\sqrt{x^2 + y^2}}, = H_0, \; \text{a constant}, \\ L(t) &= \dot{y} x - \dot{x} y, = L_0, \; \text{a constant}. \end{split} \end{split}\]

For negative energy, the solutions move along orbits $(x(t), y(t))$ that are ellipses. These are most easily described using polar coordinates $x = r \cos\theta$, $y = \sin\theta$: $$ r = \frac{L_0^2}{1 + e \cos (\theta - \theta_0)}, $$ (orbit)

with eccentricity $e = \sqrt{1 + 2 H_0 L_0^2}< 1.$

(Aside: for $H_0 > 0$ the orbits are hyperbolas and for $H_0 = 0$ they are parabolas:.)

4.3.3. Initial Tasks#

Rewrite the equation as a system of four first order equations for the quantities $$ z_1(t) = x(t),\; z_2(t) = y(t),\; z_3(t) = \dot{x}(t),\; z_4(t) = \dot{y}(t). $$
Solve these equations with the basic Euler method for the case
$$ x(0) = z_1(0) = 1-e,\; y(0) = z_2(0) = 0,\; \dot{x}(0) = z_3(0) = 0,\; \dot{y}(0) = z_4(0) = \sqrt{\frac{1+e}{1-e}} $$ <br> with eccentricity $e = 0.6$. For now, do this for $0 \leq t \leq 2\pi$ and with time step size $h = 0.0005$. <br> (It can be shown that these initial data give an exact solution of period $2\pi$, so this solution should complete one period.)
Make a parametric graph of $x(t)$ vs $y(t)$ for this computed approximation, and also display the exact elliptical orbit from Equation \eqref{orbit} on this same graph.
Plot graphs of $H(t)$ vs $t$ and $L(t)$ vs $t$. Due to numerical errors, these will not be exactly constant, and the deviations from their initial values give indications of the accuracy of the numerical solution method.

4.3.4. Further Tasks#

Once we have this basic case working with suitable graphical presentation:

Increase the time to about ten orbits ($t_{final} = 20 \pi$). Reduce the time step size $h$ as necessary to get a somewhat reasonable result; however it will still be visibly inaccurate even for very small time step size.
Solve with the other standard methods provided: the explicit trapezoid and explicit midpont methods. Compare accuracy with “equal time cost”, meaning with time step size $h$ twice as large for the latter methods as for Euler’s method (Why is this roughly equal time cost?)

A final option

Solve with the implicit midpoint method, using fixed point iteration to solve the nonlinear equations. It helps to define the increment $\delta u_n = u_{n+1} - u_n$ and write the method for solving $du/dt = f(t,u)$ as

\[\begin{split} \begin{split} \frac{\delta u_n}{h} &= f \left( t + \frac{h}{2}, u_n + \frac{\delta u_n}{2} \right) \\ u_{n+1} &= u_n + \delta u_n \end{split} \end{split}\]

so that the iteration is solving for $\delta u_n$.
For the initial approximation to start the iteration, use the approximate value of $\delta u_n$ given by Euler’s method.
How many iterations? Experiment: try a large number like ten which should give convergence within rounding error, but also try with just three iterations: the results might be interesting.
You could also experiment with a stopping condition based on the successive iterates being very close in the sense of the maximum norm.