ODE Lab II: Central Potential Problem

ODE Lab II: Central Potential Problem#

In this lab, we will implement two numerical integration methods—Runge-Kutta 4th order (RK4) and Leapfrog/Verlet—to study a 2D central potential problem. This problem models a body in an inverse-square law force field, such as planetary motion under Newtonian gravity.

The equation governing motion in a central potential is:

(403)#\[\begin{align} \frac{d^2 \mathbf{r}}{dt^2} = \mathbf{a} = - \nabla \frac{1}{|\mathbf{r}|} = - \frac{\mathbf{r}}{|\mathbf{r}|^3} \end{align}\]

This describes how a particle moves under a gravitational or electrostatic potential following an inverse-square law.

In Cartesian coordinates, this system is rewritten as:

(404)#\[\begin{align} \frac{dx}{dt} &= u, \quad \frac{du}{dt} = -\frac{x}{(x^2 + y^2)^{3/2}} \\ \frac{dy}{dt} &= v, \quad \frac{dv}{dt} = -\frac{y}{(x^2 + y^2)^{3/2}} \end{align}\]

where:

  • \((x, y)\) are the position coordinates,

  • \((u, v)\) are the velocity components,

  • The acceleration components are computed using the inverse-square law.

We will solve this system numerically using RK4 and Leapfrog/Verlet methods.

import numpy as np
from matplotlib import pyplot as plt

Fourth-Order Runge-Kutta Method (RK4)#

The RK4 method estimates the solution using four intermediate slopes:

(405)#\[\begin{align} k_1 &= \Delta t \, f(X_n, t_n) \\ k_2 &= \Delta t \, f\left(X_n + \frac{1}{2}k_1, t_n + \frac{1}{2}\Delta t\right) \\ k_3 &= \Delta t \, f\left(X_n + \frac{1}{2}k_2, t_n + \frac{1}{2}\Delta t\right) \\ k_4 &= \Delta t \, f\left(X_n + k_3, t_n + \Delta t\right) \end{align}\]

The update step is:

(406)#\[\begin{align} X_{n+1} = X_n + \frac{1}{6} k_1 + \frac{1}{3} k_2 + \frac{1}{3} k_3 + \frac{1}{6} k_4. \end{align}\]

This method is fourth-order accurate, meaning it introduces an error on the order of \(\mathcal{O}(\Delta t^4)\) per step.

# We improve the implementation from last lecture slightly

def RK4(f, x0, t0, dt, n):
    T = [t0]
    X = [x0]

    for i in range(n):
        k1 = dt * np.array(f(X[-1]         ))
        ...

        T.append(T[-1] + dt)
        X.append(...)

    return np.array(T), np.array(X)

Defining the Acceleration Function#

We define a function that computes the acceleration based on the inverse-square law:

def va(rv):
    x, y, u, v = rv
    ...
    return np.array([u, v, ...])

Now we apply RK4 to solve the system:

T, RV = RK4(va, np.array([1, 0, 0, 0.75]), 0, 0.1, 1000)

X_RK4 = RV[:,0]
Y_RK4 = RV[:,1]

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[4], line 1
----> 1 T, RV = RK4(va, np.array([1, 0, 0, 0.75]), 0, 0.1, 1000)
      3 X_RK4 = RV[:,0]
      4 Y_RK4 = RV[:,1]

Cell In[2], line 8, in RK4(f, x0, t0, dt, n)
      5 X = [x0]
      7 for i in range(n):
----> 8     k1 = dt * np.array(f(X[-1]         ))
      9     ...
     11     T.append(T[-1] + dt)

TypeError: unsupported operand type(s) for *: 'float' and 'ellipsis'

plt.scatter([0],[0], color='k')
plt.plot(X_RK4, Y_RK4)
plt.xlim(-1.2,1.2)
plt.ylim(-1.2,1.2)
plt.gca().set_aspect('equal')

Leapfrog/Verlet Method#

Leapfrog integrates the equations in a staggered manner:

(407)#\[\begin{align} v_{n+1/2} &= v_n + \frac{1}{2} a(r_{n}) \Delta t, \\ r_{n+1} &= r_n + v_{n+1/2} \Delta t, \\ v_{n+1} &= v_{n+1/2} + \frac{1}{2} a(r_{n+1}) \Delta t. \end{align}\]
# Leapfrog/Verlet Method
def leapfrog(a, r0, v0, t0, dt, n):
    
    T = [t0]
    X = [r0[0]]
    Y = [r0[1]]
    U = [v0[0]]
    V = [v0[1]]

    for _ in range(n):
        ...

    return np.array(T), np.array([X, Y]).T, np.array([U, V]).T

With leapfrog(), the “force” term is easier to implement:

def a(x, y):
    ...

Running the leapfrog integration:

T, R, V = leapfrog(a, (1, 0), (0, 0.75), 0, 0.1, 1000)

X_lf = R[:,0]
Y_lf = R[:,1]

plt.scatter([0],[0], color='k')
plt.plot(X_lf,  Y_lf)
plt.xlim(-1.2,1.2)
plt.ylim(-1.2,1.2)
plt.gca().set_aspect('equal')

Comparing RK4 and Leapfrog#

Now integrate this up to 2,000 dynamic time and compare between RK4 and leapfrog.

T, RV = RK4(va, np.array([1, 0, 0, 0.75]), 0, 0.1, 20_000)
X_RK4 = RV[:,0]
Y_RK4 = RV[:,1]

T, R, V = leapfrog(a, (1, 0), (0, 0.75), 0, 0.1, 20_000)
X_lf = R[:,0]
Y_lf = R[:,1]

plt.scatter([0],[0], color='k')
plt.plot(X_lf,  Y_lf)
plt.plot(X_RK4, Y_RK4)
plt.xlim(-1.2,1.2)
plt.ylim(-1.2,1.2)
plt.gca().set_aspect('equal')

Energy Computation and Conservation Analysis#

To analyze the stability of our integrators, we compute the total energy:

(408)#\[\begin{align} E = \frac{1}{2} (u^2 + v^2) - \frac{1}{\sqrt{x^2 + y^2}} \end{align}\]
def energy(R, V):
    ...

E_rk4 = energy(RV[:,:2], RV[:,2:])
E_lf  = energy(R, V)

plt.plot(T, E_rk4, label='RK4 Energy', linestyle='--')
plt.plot(T, E_lf,  label='Leapfrog Energy')
plt.xlabel('Time')
plt.ylabel('Total Energy')
plt.ylim(-10,10)
plt.legend()

Conclusion#

  • RK4 provides high accuracy but does not conserve energy perfectly.

  • Leapfrog maintains energy conservation better, making it suitable for long-term orbital simulations.

Contents