mathematics//numerical methods

Numerical methods are the algorithms that turn continuous mathematics into a finite sequence of arithmetic operations a computer can carry out (stepping a differential equation through time, solving a large linear system, taking a derivative of a program), and their craft is getting an answer without the rounding and the approximation eating it. Every simulator, every state estimator in an autopilot and every training run of a network is a numerical method, and most of their strange failures (a simulated oscillator that gains energy, a filter whose covariance goes negative, a training that suddenly diverges) are numerical before they are physical.


Numerical methods are the algorithms that turn continuous mathematics into a finite sequence of arithmetic operations a computer can carry out (stepping a differential equation through time, solving a large linear system, taking a derivative of a program), and their craft is getting an answer without the rounding and the approximation eating it. Every simulator, every state estimator in an autopilot and every training run of a network is a numerical method, and most of their strange failures (a simulated oscillator that gains energy, a filter whose covariance goes negative, a training that suddenly diverges) are numerical before they are physical.

The members of the family split by the question they answer. To advance a model in time, numerical integration chooses a step and a method, from the one-line Euler method that a microcontroller runs at a kilohertz to the Runge-Kutta method that buys accuracy with four evaluations per step; numerical stability says whether the method stays bounded when the system does, and stiffness is the case where a fast mode forces tiny steps on a slow simulation. To get a derivative, finite differences perturb the input and divide (enough for a handful of parameters), and automatic differentiation transforms the program to get exact gradients for millions. To keep a statistic inside a loop with fixed memory, Welford's algorithm updates a mean and a variance sample by sample without the cancellation of the textbook formula. And before any of these, variable scaling puts the numbers in similar ranges, which on its own fixes a surprising share of problems.

Numerical error is born in rounding, grows with conditioning, enters with the data and travels downstream.

Rounding is about 10−1610^{-16}10−16 in double precision and 10−710^{-7}10−7 in single (floating-point arithmetic); a badly conditioned problem amplifies it by its condition number; a covariance estimated from few samples brings its own distorted eigenvalues; and an ill-conditioned calibration matrix turns sensor noise into a wrong RRR for the Kalman filter, and from there into the control.

The linear algebra underneath is already solved. Behind NumPy's linalg sit BLAS and LAPACK, kernels tuned for decades (OpenBLAS, Intel MKL), and the habits that matter are about using them well: solve a system with a factorization instead of forming an inverse, use least squares routines instead of the normal equations, which square the condition number (linear system of equations). On a microcontroller the equivalent is a small library such as CMSIS-DSP, with matrices of a few dozen entries in single precision.

The simplest method that works is usually right. For a few seconds of a smooth system, a library solver with default settings is enough; a hand-written fixed-step integrator earns its place when it must run on a microcontroller or in a hardware-in-the-loop rig with deadlines (flyswatter rule). Where the hardware has no floating point, scaling becomes part of the design (fixed-point arithmetic).

The modelling side of turning time into steps, and what it changes in a system's behaviour, is discretization; the use of all of this to ask what-if questions is simulation.