Numerical integration and ODEs

Integrals and differential equations are the language of engineering, but most of the ones that matter cannot be solved by hand. The discharge of a river is an integral of measured velocities across the channel. The energy in a load cycle is the integral of a recorded power signal. The motion of a vehicle, the temperature of a battery pack and the concentration in a reactor follow differential equations with nonlinear terms that defeat the methods of MATH-201. This chapter covers the methods that compute them anyway: the trapezoidal and Simpson rules, Gauss quadrature, and the Euler and Runge–Kutta methods for ordinary differential equations, along with the idea of stiffness that decides which ODE solver to use.

Numerical integration (quadrature)

The idea is to replace the integrand by an interpolating polynomial (Lesson 4) on each small piece of the interval and integrate that exactly.

The trapezoidal rule

Joining neighboring points with straight lines and adding the areas of the trapezoids gives, for n equal panels of width h,

Its error is about \(\frac{(b - a)h^2}{12}|f''|\): second order, so halving h quarters the error. It is exact for straight lines and works directly on measured data at uneven spacing.

Simpson's rule

Fitting a parabola through each set of three points (an even number of panels) gives Simpson's 1/3 rule:

Its error is about \(\frac{(b - a)h^4}{180}|f^{(4)}|\): fourth order, so halving h divides the error by 16. Surprisingly, it is exact for cubics as well as parabolas.

Richardson extrapolation and adaptive methods

If the error of a rule is known to be proportional to \(h^2\), two estimates with h and h/2 can be combined to cancel the leading error: \(I \approx (4T_{h/2} - T_h)/3\). Applied to the trapezoidal results for sin x with 4 and 8 panels, this gives 2.00027, which is exactly Simpson's rule with 8 panels. Repeating the idea gives Romberg integration. Library routines go further: adaptive quadrature (scipy.integrate.quad, MATLAB integral) estimates the error on each subinterval and refines only where the integrand is difficult, delivering a requested accuracy with as few evaluations as possible.

Gauss quadrature

Newton–Cotes rules (trapezoid, Simpson) use equally spaced points. Gauss quadrature also chooses where to evaluate, picking the points and weights to integrate the highest possible polynomial degree exactly: n points integrate polynomials up to degree 2n − 1 exactly. On [−1, 1] the two-point rule is

Mapped to \([0, \pi]\), two evaluations of sin x give 1.936 and three give 2.0014, better than Simpson with four panels and five evaluations. Gauss quadrature is the standard way finite element programs integrate stiffness matrices over each element, where every evaluation costs real time.

Ordinary differential equations

An initial value problem \(y' = f(t, y)\), \(y(t_0) = y_0\) is solved numerically by stepping forward in time with step h, using the slope f to predict the next value.

Euler's method

The simplest step follows the tangent line:

ENCM-100 used it to simulate a draining tank. It is first order: the global error after reaching a fixed time is proportional to h, so halving the step halves the error and doubles the work. That is too slow for accurate work.

Runge–Kutta methods

Runge–Kutta methods sample the slope at several points within each step and average them. The improved Euler (Heun) method predicts with Euler, then averages the slopes at both ends (second order). The classical fourth-order Runge–Kutta method (RK4) is the most famous:

Four slope evaluations per step buy fourth-order accuracy: halving h divides the error by 16.

Systems and higher-order equations

Every method applies unchanged to systems, with y and f as vectors. Higher-order equations are first rewritten as first-order systems. The mass–spring–damper of MATH-201, \(mx'' + cx' + kx = F(t)\), becomes

This is how vehicle suspensions, robot arms, orbits, chemical reactors and power-system transients are simulated.

Stability and stiff equations

A method can be accurate in principle and still blow up. Applied to the decay equation \(y' = -\lambda y\) (λ > 0), Euler gives \(y_{k+1} = (1 - h\lambda)y_k\), which decays only if \(|1 - h\lambda| < 1\), that is \(h < 2/\lambda\). Larger steps make the numerical solution oscillate and grow without bound, even though the true solution decays to zero.

An equation is stiff when it contains a very fast-decaying component alongside the slow behavior of interest. Explicit methods (Euler, RK4) must then take tiny steps to stay stable, long after the fast component has died away and accuracy no longer needs them. Stiffness is common: chemical kinetics with fast and slow reactions, circuits with very different time constants, and the heat equation on a fine grid.

Implicit methods must solve an equation (often a nonlinear system, by Newton's method) at every step, so each step costs more, but they can take far larger steps on stiff problems. Library solvers make the choice easy: in scipy.integrate.solve_ivp use RK45 (the default, explicit, adaptive) for ordinary problems and Radau or BDF for stiff ones; MATLAB's equivalents are ode45 and ode15s. If an explicit solver is crawling with tiny steps, suspect stiffness.

import numpy as np
from scipy.integrate import simpson, quad, solve_ivp

x = np.arange(0, 17, 2.0)
q = np.array([0, 0.32, 0.81, 1.25, 1.42, 1.30, 0.94, 0.48, 0])
print(f"discharge (Simpson): {simpson(q, x=x):.2f} m3/s")
print(f"integral of sin on [0, pi]: {quad(np.sin, 0, np.pi)[0]:.10f}")

drag = lambda t, v: 9.81 - (0.25/80)*v**2
sol = solve_ivp(drag, [0, 4], [0.0], rtol=1e-8, atol=1e-10)
print(f"skydiver v(4) = {sol.y[0, -1]:.3f} m/s")

stiff = lambda t, y: -50*(y - np.cos(t))
for method in ["RK45", "BDF"]:
    s = solve_ivp(stiff, [0, 10], [0.0], method=method, rtol=1e-6, atol=1e-9)
    print(f"{method}: y(10) = {s.y[0, -1]:.6f}, function evaluations = {s.nfev}")

Output: discharge (Simpson): 13.16 m3/s; integral of sin on [0, pi]: 2.0000000000; skydiver v(4) = 33.875 m/s; RK45: y(10) = −0.849612, function evaluations = 2942; BDF: y(10) = −0.849612, function evaluations = 428.