Numerical solution of 2BP

Equation of motion

The equation of motion for the two-body problem (2BP) is given by:

\[ \ddot{\mathbf{r}} = -\mu \frac{\mathbf{r}}{r^3} \]

where, \(\mathbf{r}\) is the position vector, \(r = |\mathbf{r}|\) is the magnitude of the position vector, \(\mu = G(M + m)\) is the standard gravitational parameter. For Earth, \(\mu \approx 398600 \, \text{km}^3/\text{s}^2\).

RK4 Method

The Runge-Kutta 4th order (RK4) method is a numerical technique used to solve ordinary differential equations (ODEs). The method calculates the solution at the next time step using a weighted average of slopes (derivatives) evaluated at intermediate points.

For a system of ODEs \(\dot{\mathbf{y}} = \mathbf{f}(t, \mathbf{y})\), the RK4 method updates the solution as follows:

1. Compute intermediate slopes: \[ \begin{aligned} k_1 &= \mathbf{f}(t_n, \mathbf{y}_n), \\ k_2 &= \mathbf{f}\left(t_n + \frac{h}{2}, \mathbf{y}_n + \frac{h}{2} k_1\right), \\ k_3 &= \mathbf{f}\left(t_n + \frac{h}{2}, \mathbf{y}_n + \frac{h}{2} k_2\right), \\ k_4 &= \mathbf{f}(t_n + h, \mathbf{y}_n + h k_3). \end{aligned} \]

2. Update the solution: \[ \mathbf{y}_{n+1} = \mathbf{y}_n + \frac{h}{6} (k_1 + 2k_2 + 2k_3 + k_4). \]

Here, \(h\) is the time step, and \(t_n\) and \(\mathbf{y}_n\) are the current time and state vector, respectively.

Implementation of RK4 for 2BP

Here is a Python code that implements the RK4 method for the two-body problem:

# Converted from MATLAB code by Suyog Soman (MTech, AFM 2023-2025) with GPT-4o

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# Standard gravitational parameter for Earth (km^3/s^2)
mu = 398600.0

def f(t, y):
    """
    Computes the derivatives for the two-body problem.
    
    Parameters:
        t : float
            Current time (not used explicitly as the problem is time invariant).
        y : array_like
            State vector [rx, ry, rz, vx, vy, vz].
    
    Returns:
        dydt : ndarray
            Derivatives [vx, vy, vz, ax, ay, az].
    """
    # Unpack the state vector
    rx, ry, rz, vx, vy, vz = y
    # Compute the distance from the origin
    r = np.sqrt(rx**2 + ry**2 + rz**2)
    # Compute the acceleration components
    ax = -mu * rx / r**3
    ay = -mu * ry / r**3
    az = -mu * rz / r**3
    return np.array([vx, vy, vz, ax, ay, az])

def rk4_step(t, y, h):
    """
    Performs one step of the RK4 method.
    
    Parameters:
        t : float
            Current time.
        y : array_like
            Current state vector.
        h : float
            Time step.
    
    Returns:
        y_next : ndarray
            Updated state vector after time h.
    """
    k1 = f(t, y)
    k2 = f(t + h/2, y + h/2 * k1)
    k3 = f(t + h/2, y + h/2 * k2)
    k4 = f(t + h, y + h * k3)
    
    return y + (h/6) * (k1 + 2*k2 + 2*k3 + k4)

def simulate_orbit(y0, t0, tf, h):
    """
    Simulates the orbit using the RK4 integration method for the two-body problem.
    
    Parameters:
        y0 : array_like
            Initial state vector [rx, ry, rz, vx, vy, vz].
        t0 : float
            Initial time.
        tf : float
            Final time.
        h : float
            Time step.
    
    Returns:
        t_values : ndarray
            Array of time points.
        states : ndarray
            Array of state vectors for each time point.
    """
    # Calculate the number of time steps
    num_steps = int(np.ceil((tf - t0) / h))
    # Preallocate arrays for time and state
    t_values = np.linspace(t0, t0 + num_steps * h, num_steps + 1)
    states = np.zeros((num_steps + 1, len(y0)))
    states[0] = y0
    
    t = t0
    y = y0.copy()
    for i in range(1, num_steps + 1):
        y = rk4_step(t, y, h)
        t += h
        states[i] = y
    return t_values, states

if __name__ == "__main__":
    # Initial conditions: 
    r0 = 7000.0  # km (distance from Earth's center)
    v0 = np.sqrt(mu / r0)  # circular orbital speed (km/s)
    
    # Initial state: starting at (r0, 0, 0) with velocity (0, v0, 0)
    y0 = np.array([r0, 0.0, 0.0, 0.0, v0, 0.1])
    t0 = 0.0          # initial time in seconds
    tf = 7200.0       # simulate for 5400 seconds (2 hours)
    h = 10.0          # time step in seconds
    
    # Run the simulation
    t_values, states = simulate_orbit(y0, t0, tf, h)
    
    # Extract positions for plotting
    rx_vals = states[:, 0]
    ry_vals = states[:, 1]
    rz_vals = states[:, 2]
    
    # Plot the orbit in 3D
    fig = plt.figure(figsize=(10,8))
    ax = fig.add_subplot(111, projection='3d')
    ax.plot(rx_vals, ry_vals, rz_vals, label='Orbit')
    ax.scatter([0], [0], [0], color='red', label='Earth')  # Earth's position at the origin
    
    ax.set_xlabel('rx [km]')
    ax.set_ylabel('ry [km]')
    ax.set_zlabel('rz [km]')
    ax.set_title('Orbit using RK4 Method')
    ax.legend()
    plt.show()