The Reality of ODEs in Practice

Most DEs arising in science and engineering cannot be solved analytically:

• The force \(F(x,x')\) is nonlinear (pendulum, van der Pol, …)
• Parameters are time-varying or data-driven
• The model is a large coupled system

Even when an analytic solution exists, it may involve complicated integrals or infinite series requiring numerical evaluation anyway.

Discretization

Divide \([t_0, T]\) into \(N\) equal subintervals of length

\[h = \frac{T - t_0}{N} \qquad \text{(step size)}.\]

Grid: \(t_n = t_0 + nh\), \(n = 0,1,\ldots,N\).

Derivation

Integrate \(x'(t)=f(t,x(t))\) from \(t_n\) to \(t_{n+1}\):

\[\int_{t_n}^{t_{n+1}} x'(t)\,dt = \int_{t_n}^{t_{n+1}} f(t,x(t))\,dt. \tag{6.4}\]

Left side: \(x(t_{n+1})-x(t_n)\) (Fundamental Theorem).

Right side: left-hand rectangle rule — approximate \(f\) by its value at \(t_n\):

\[x(t_{n+1}) = x(t_n) + h\,f(t_n,x(t_n)) + \underbrace{O(h^2)}_{\text{local error}}.\]

Geometric Picture — The Animation

At each step: follow the slope field at the current point for distance \(h\).

Euler method: \(x'=t-x\), \(x(0)=1\), \(h=0.25\). Green = Euler step (local tangent); orange = local error; blue = exact solution; red = numerical approximation.

Error Analysis

Example 6.3 — Euler Convergence

Modified Euler (RK2) — Predictor–Corrector

Idea: use the trapezoidal rule instead of the rectangle rule for the integral in (6.4):

\[x(t_{n+1})-x(t_n) = \frac{h}{2}\bigl[f(t_n,x(t_n))+f(t_{n+1},x(t_{n+1}))\bigr]+O(h^3).\]

This gives an implicit equation. Fix it with a predictor–corrector:

The Runge–Kutta Method (RK4)

The workhorse of numerical ODEs. Evaluates \(f\) at four points per step:

Method Comparison

Error Sources and Stiff Equations

Three sources of error in any finite difference scheme:

• Truncation error — from approximating a continuous problem by a discrete one. This is the \(O(h^p)\) analysis we have been doing.

• Stability error — the discrete algorithm itself may be unstable. An early error can be amplified exponentially as it propagates forward. The step size \(h\) must stay below a problem-dependent threshold.

• Roundoff error — computers represent real numbers with finite precision. The actually computed \(\overline{X}_n\) differs from the exact discrete \(X_n\); reducing \(h\) too far eventually makes roundoff dominate.

To help develop our intuition for these numerical methods, let’s view this visual web app.

Extending to Systems (§6.3)

For the 2D system \(x'=f(t,x,y)\), \(y'=g(t,x,y)\), \(x(t_0)=x_0\), \(y(t_0)=y_0\):

Euler: \[X_{n+1}=X_n+hf(t_n,X_n,Y_n), \quad Y_{n+1}=Y_n+hg(t_n,X_n,Y_n).\]

Modified Euler — predictors then correctors (same structure, applied to both components simultaneously).

RK4 — compute eight slopes \(k_{11},k_{21},k_{12},k_{22},k_{13},k_{23},k_{14},k_{24}\) (one pair per component at each of the four evaluation points), then:

\[X_{n+1}=X_n+\frac{h}{6}(k_{11}+2k_{12}+2k_{13}+k_{14}),\] \[Y_{n+1}=Y_n+\frac{h}{6}(k_{21}+2k_{22}+2k_{23}+k_{24}).\]

Why Use solve_ivp?

Rather than coding methods from scratch, use scipy.integrate.solve_ivp:

• Adaptive step size: automatically tightens \(h\) where the solution changes rapidly, relaxes it where it is smooth.
• Multiple solvers: 'RK45' (default), 'RK23', 'DOP853', 'Radau' (stiff), 'BDF' (stiff).
• Dense output: continuous interpolant sol.sol(t) for plotting.
• Event detection: find zero-crossings (e.g., when a population hits zero).

from scipy.integrate import solve_ivp

sol = solve_ivp(
    fun    = lambda t, y: [f(t, y[0])],  # rhs: (t, y-vector) → dy/dt vector
    t_span = (t0, T),                     # integration interval
    y0     = [x0],                        # initial condition (list/array)
    method = 'RK45',                      # solver
    dense_output = True,                  # enable sol.sol(t) interpolant
    rtol   = 1e-6,                        # relative tolerance
    atol   = 1e-8,                        # absolute tolerance
)

solve_ivp — Single Equation

sol = solve_ivp(
    fun          = lambda t, y: [t - y[0]],
    t_span       = (0, 2),
    y0           = [1.0],
    method       = 'RK45',
    dense_output = True,
    rtol=1e-6, atol=1e-8,
)
print(f"Steps: {len(sol.t)}   x(2) ≈ {sol.y[0,-1]:.8f}   exact: {exact_ex(2):.8f}")

Steps: 11   x(2) ≈ 1.27067069   exact: 1.27067057

solve_ivp — System Example (Lotka–Volterra)

a, b, c, d = 0.6, 0.5, 0.3, 0.4
sol_lv = solve_ivp(
    fun = lambda t, z: [a*z[0]-b*z[0]*z[1],   # x' = ax - bxy
                        -c*z[1]+d*z[0]*z[1]],  # y' = -cy + dxy
    t_span = (0, 60),
    y0     = [c/d + 0.4, a/b],
    method = 'RK45',
    dense_output = True,
    rtol=1e-9, atol=1e-11, max_step=0.05,
)

Solver Selection Guide

Key Takeaways

• The Euler method (6.5) \(X_{n+1}=X_n+hf(t_n,X_n)\) follows the slope field step by step. It is first-order: global error \(O(h)\), so halving \(h\) halves the error.

• The modified Euler / RK2 method uses a predictor–corrector pair (6.7)–(6.8). Global error \(O(h^2)\): halving \(h\) quarters the error.

• Runge–Kutta RK4 uses a weighted average of four slope evaluations \(k_1,\ldots,k_4\). Global error \(O(h^4)\): halving \(h\) reduces error by a factor of 16. RK4 is the standard workhorse of numerical ODE solving.

• All three methods extend to systems by treating the state as a vector and applying the formula component-wise. Error orders are unchanged.

• Error sources: truncation (Taylor series), stability (scheme may amplify errors), and roundoff (finite-precision arithmetic). Stiff equations require implicit solvers.

• scipy.integrate.solve_ivp provides high-quality adaptive RK solvers with automatic step-size control. Use dense_output=True for a continuous interpolant, method='Radau' for stiff problems, and always tighten rtol/atol for accuracy-sensitive work.