Numerical ODE Solvers

Integrate differential equations step by step. Learn Euler's method, modified Euler/Heun (RK2), classical RK4, local and global error, stability regions and stiff equations, multistep (Adams) and predictor–corrector methods, adaptive step size, and solving systems — with Euler vs RK4 accuracy tests.

Most real ODEs (non-linear vehicle dynamics, chemical kinetics, circuit simulators) have no closed-form solution. Numerical solvers such as MATLAB's ode45 or SciPy's solve_ivp are built from the methods in this lesson.

An analogy

Euler's method is like navigating by walking straight in the direction you currently face for a fixed distance, then re-checking; RK4 looks ahead at several points before committing to each step.

Euler's method

Uses the slope at the start of each step. Simple, but first-order accurate and can be unstable for stiff or oscillatory problems.

Runge–Kutta methods

Evaluate slopes at intermediate points and average them. RK4 balances accuracy and cost and is the classical default for smooth, non-stiff problems.

Systems and higher order

Convert to first-order systems (state vectors) and apply the same formulas component-wise — pendulums, vehicle dynamics, circuits.

Stiffness, adaptivity, multistep

Stiff systems (fast and slow time scales) force explicit methods to tiny steps; implicit methods (backward Euler, BDF — ode15s) remain stable. Adaptive solvers (Dormand–Prince RK45) estimate error and change h. Adams methods reuse past steps.

Formulas

Worked example: One Euler step

Given: y′ = −2y + t, y(0) = 1, h = 0.1

Find: y(0.1).

Answer: 0.8 (error 0.023)

Worked example: Stability limit

Given: y′ = −50y with Euler

Find: Maximum stable h.

Answer: h ≤ 0.04 s (beyond this the solution oscillates and grows)

Euler vs RK4 convergence test

import math

# y' = -2y + t, y(0) = 1 ; exact y = 0.5t - 0.25 + 1.25 e^(-2t)
f = lambda t, y: -2 * y + t
exact = lambda t: 0.5 * t - 0.25 + 1.25 * math.exp(-2 * t)

def euler(h):
    t, y = 0.0, 1.0
    for _ in range(round(1 / h)):
        y += h * f(t, y); t += h
    return y

def rk4(h):
    t, y = 0.0, 1.0
    for _ in range(round(1 / h)):
        k1 = f(t, y); k2 = f(t + h/2, y + h*k1/2)
        k3 = f(t + h/2, y + h*k2/2); k4 = f(t + h, y + h*k3)
        y += h * (k1 + 2*k2 + 2*k3 + k4) / 6; t += h
    return y

for h in (0.1, 0.05):
    print(f"h={h}: Euler err {abs(euler(h) - exact(1)):.2e}   RK4 err {abs(rk4(h) - exact(1)):.2e}")
# h=0.1: Euler err 3.50e-02   RK4 err 5.33e-06
# h=0.05: Euler err 1.72e-02   RK4 err 3.06e-07

Halving h halves Euler's error (first order) but divides RK4's by ≈ 17 (≈ 2⁴, fourth order). This 'halve-the-step' test verifies both the method and the code.

Common mistakes

Applications

Key takeaways

Going further