Runge-Kutta Method

Differential Equations

The Runge-Kutta methods are a family of iterative methods for solving ordinary differential equations. The fourth-order Runge-Kutta method (RK4) is one of the most widely used numerical methods for ODEs due to its balance between accuracy and computational efficiency.

Mathematical Formulation

Given an initial value problem:

\frac{d y}{d t} = f (t, y), y (t_{0}) = y_{0}

The RK4 method approximates the solution at $t_{n + 1} = t_{n} + h$ using four evaluations of the function $f$ :

k_{1} k_{2} k_{3} k_{4} y_{n + 1} = f (t_{n}, y_{n}) = f (t_{n} + \frac{h}{2}, y_{n} + \frac{h}{2} k_{1}) = f (t_{n} + \frac{h}{2}, y_{n} + \frac{h}{2} k_{2}) = f (t_{n} + h, y_{n} + h k_{3}) = y_{n} + \frac{h}{6} (k_{1} + 2 k_{2} + 2 k_{3} + k_{4})

The method is fourth-order accurate, meaning the local truncation error is $O (h^{5})$ and the global error is $O (h^{4})$ .

Derivation Intuition

RK4 can be thought of as a weighted average of four slope estimates:

$k_{1}$ : slope at the beginning of the interval
$k_{2}$ : slope at the midpoint using $k_{1}$
$k_{3}$ : slope at the midpoint using $k_{2}$
$k_{4}$ : slope at the end using $k_{3}$

The final step uses Simpson's rule-like weighting: $\frac{1}{6} (k_{1} + 2 k_{2} + 2 k_{3} + k_{4})$ .

Python Implementation

The following code implements RK4 for solving ODEs:

import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl

def set_publication_style():
    """Set publication-quality matplotlib style."""
    mpl.rcParams.update({
        'font.family': 'serif',
        'font.size': 12,
        'axes.labelsize': 14,
        'axes.titlesize': 16,
        'axes.linewidth': 1.2,
        'axes.labelpad': 8,
        'axes.titlepad': 10,
        'xtick.labelsize': 12,
        'ytick.labelsize': 12,
        'xtick.direction': 'in',
        'ytick.direction': 'in',
        'xtick.top': True,
        'ytick.right': True,
        'xtick.major.size': 6,
        'ytick.major.size': 6,
        'xtick.major.width': 1.2,
        'ytick.major.width': 1.2,
        'legend.fontsize': 12,
        'legend.frameon': False,
        'lines.linewidth': 2,
        'lines.markersize': 6,
        'figure.dpi': 100,
        'savefig.dpi': 300,
        'savefig.bbox': 'tight'
    })

set_publication_style()

def rk4(f, y0, t0, t_end, dt, *args):
    """
    Runge-Kutta 4th order method for solving ODEs.
    
    Parameters:
    -----------
    f : callable
        Function f(t, y, *args) that returns dy/dt
    y0 : float or array
        Initial value y(t0)
    t0 : float
        Start time
    t_end : float
        End time
    dt : float
        Time step size
    *args : additional arguments
        Additional arguments to pass to f
    
    Returns:
    --------
    t : array
        Time points
    y : array
        Solution values at each time point
    """
    n_steps = int((t_end - t0) / dt)
    t = np.linspace(t0, t_end, n_steps + 1)
    y = np.zeros(n_steps + 1)
    y[0] = y0

    for i in range(n_steps):
        t_i = t[i]
        y_i = y[i]
        
        # Compute the four k values
        k1 = f(t_i, y_i, *args)
        k2 = f(t_i + dt / 2, y_i + dt * k1 / 2, *args)
        k3 = f(t_i + dt / 2, y_i + dt * k2 / 2, *args)
        k4 = f(t_i + dt, y_i + dt * k3, *args)
        
        # Update solution
        y[i + 1] = y_i + dt * (k1 + 2 * k2 + 2 * k3 + k4) / 6

    return t, y

# Example: Exponential growth
def exponential_growth(t, y, k):
    """dy/dt = k*y"""
    return k * y

# Parameters
k = 1.0  # Growth rate
y0 = 1.0  # Initial value
t0 = 0.0  # Start time
t_end = 5.0  # End time
dt = 0.1  # Time step

# Solve using RK4
t_rk4, y_rk4 = rk4(exponential_growth, y0, t0, t_end, dt, k)

# Analytical solution for comparison
t_analytical = np.linspace(t0, t_end, 1000)
y_analytical = y0 * np.exp(k * t_analytical)

# Plot results
plt.figure(figsize=(10, 6))
plt.plot(t_rk4, y_rk4, 'o-', label='RK4 Solution', markersize=4)
plt.plot(t_analytical, y_analytical, '--', label='Analytical Solution', linewidth=2)
plt.xlabel('Time (t)')
plt.ylabel('y(t)')
plt.title('Exponential Growth: RK4 vs Analytical Solution')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.savefig('figures/runge_kutta_comparison.png', dpi=300, bbox_inches='tight')
plt.show()

# Calculate error
error = np.abs(y_rk4 - y0 * np.exp(k * t_rk4))
print(f"Maximum error: {np.max(error):.2e}")
print(f"Mean error: {np.mean(error):.2e}")

Comparison with Other Methods

RK4 offers a good balance between accuracy and computational cost:

Euler's Method: First-order ( $O (h)$ ), requires 1 function evaluation per step
RK2 (Midpoint): Second-order ( $O (h^{2})$ ), requires 2 function evaluations per step
RK4: Fourth-order ( $O (h^{4})$ ), requires 4 function evaluations per step

For many problems, RK4 provides excellent accuracy with reasonable computational cost, making it the method of choice for many applications.

Advantages and Limitations

Advantages:

High accuracy (fourth-order)
Stable for a wide range of problems
Self-starting (doesn't require previous solution values)
Efficient for non-stiff problems

Limitations:

Not optimal for stiff ODEs (may require very small step sizes)
Fixed step size in basic implementation (adaptive versions exist)
Requires four function evaluations per step