Numerical Methods

The ideas behind numerical methods

Many engineering problems have no neat formula: the equation cannot be rearranged, the integral has no closed form, the system is too large to solve by hand. Numerical methods replace the exact answer by a sequence of simple calculations that get closer to it, and they come with three questions you should always ask: how fast does it converge (the error against the effort), when does it fail (bad starting points, flat slopes, instability), and how do I know it is right (a check against a known answer or a second method).

Source of error	What it is	How to reduce it
Truncation (method) error	The method is an approximation: a finite number of strips or steps	A smaller step, or a higher-order method
Rounding error	Floating-point numbers are approximate (see Numbers and Floating Point)	Stable formulas, double precision, fewer operations
Model error	The equations do not describe reality exactly	Better physics; compare with measurements
Data error	The inputs are measured with uncertainty	Better data; carry the uncertainty through

Solving equations: bisection and Newton-Raphson

To solve f(x) = 0 there are two workhorses. Bisection needs an interval [a, b] where f changes sign: halve it, keep the half that still changes sign, and repeat. It always works and gains one binary digit per step. Newton-Raphson starts from a guess and follows the tangent to the axis:

x_next = x − f(x) ÷ f′(x)

def newton(f, df, x, tol=1e-12, max_iter=50):
    for _ in range(max_iter):                 # always bound the loop
        step = f(x) / df(x)
        x -= step
        if abs(step) < tol * max(1.0, abs(x)):
            return x
    raise RuntimeError("did not converge")

newton(lambda x: x*x - 2, lambda x: 2*x, 2.0)       # 1.4142135623730951

Integrating: trapezoid and Simpson

To find an area (energy from a power log, distance from a speed log, work from a force curve), split it into strips. The trapezoid rule joins the points by straight lines; Simpson’s rule fits a parabola through each pair of strips. For n strips of width h = (b − a) ÷ n:

Trapezoid: h × [ ½ f₀ + f₁ + f₂ + … + ½ f_n ]    Simpson: (h ÷ 3) × [ f₀ + 4 f₁ + 2 f₂ + 4 f₃ + … + f_n ] (n even)

Strips n	Trapezoid	Error	Simpson	Error
4	1.8961	1.04e-1	2.0046	4.56e-3
8	1.9742	2.58e-2	2.0003	2.69e-4
16	1.9936	6.43e-3	2.0000	1.66e-5
32	1.9984	1.61e-3	2.0000	1.03e-6

The integral of sin(x) from 0 to π is exactly 2. Doubling the strips cuts the trapezoid error by about four (it falls as h²) and the Simpson error by about sixteen (it falls as h⁴): Simpson reaches 1.03e-6 with 32 strips where the trapezoid is at 1.61e-3.

def trapezoid(f, a, b, n):
    h = (b - a) / n
    total = 0.5 * (f(a) + f(b))
    for i in range(1, n):
        total += f(a + i * h)
    return total * h

Differentiating numerically

The slope of f is estimated from nearby values: the forward difference (f(x + h) − f(x)) ÷ h, or the more accurate central difference (f(x + h) − f(x − h)) ÷ (2h). Two errors fight: the truncation error shrinks with h, but the rounding error from subtracting nearly equal numbers grows as h shrinks (cancellation, see Numbers and Floating Point).

Differential equations: Euler and Runge-Kutta

A differential equation dy/dt = f(t, y) describes how a quantity changes; to find y(t) we step forward in time. Euler’s method takes one step along the slope: y_next = y + h × f(t, y). It is simple and first-order accurate (the error falls with h). Runge-Kutta 4 (RK4) evaluates the slope four times per step and is fourth-order accurate, far better for the work.

import numpy as np
from scipy.integrate import solve_ivp

k, Ta = 0.1, 20.0
sol = solve_ivp(lambda t, T: -k * (T - Ta), (0, 60), [90.0], rtol=1e-8, atol=1e-10, dense_output=True)
sol.sol(10.0)          # temperature at t = 10 s: 45.75 C (exact: 20 + 70 e^-1)

Interpolation: reading between the points

When a table gives points but you need a value in between, draw a straight line between the neighboring points:

y = y₀ + (x − x₀) × (y₁ − y₀) ÷ (x₁ − x₀)

If a table gives 4.0 at 20 and 6.5 at 30, then at 24 the value is 4.0 + (24 − 20) × (6.5 − 4.0) ÷ (30 − 20) = 5.00. The error of linear interpolation is at most h² ÷ 8 × max|f″|: for a sine table at 10° spacing that is (0.1745)² ÷ 8 = 0.0038 (0.4 % of full scale).

• Never extrapolate far outside the table: the line (or polynomial) can run away.
• Higher-degree polynomials through many points oscillate (Runge’s phenomenon) and are rarely a good idea; use piecewise linear or cubic splines (scipy.interpolate.CubicSpline), which are smooth and well behaved.
• Table lookups on a microcontroller: a table with linear interpolation is the standard way to linearize a thermistor or a sensor curve (see Embedded Systems).
• Interpolation passes through the points; fitting does not: if the points are noisy, fit a curve instead (next).

def interp(x, x0, y0, x1, y1):
    return y0 + (x - x0) * (y1 - y0) / (x1 - x0)

interp(24, 20, 4.0, 30, 6.5)          # 5.0
np.interp(24, [20, 30], [4.0, 6.5])   # the same, for tables: xp must be increasing

Curve fitting by least squares

To summarize noisy measurements with a model, choose the model’s parameters that minimize the sum of the squared residuals (the vertical distances from the points to the curve). For a straight line y = m x + c there is a closed form:

m = (n Σxy − Σx Σy) ÷ (n Σx² − (Σx)²)    c = (Σy − m Σx) ÷ n

• Polynomials: np.polyfit(x, y, degree). A degree that is too high chases the noise (overfitting): with n points a polynomial of degree n − 1 passes through all of them and is useless between them. Choose the lowest degree that leaves random residuals.
• Exponentials and power laws become straight lines after taking logarithms: fit ln y = ln A − t ÷ τ, then convert back. This weights small values heavily, so for a final answer use nonlinear least squares.
• Nonlinear least squares (scipy.optimize.curve_fit, least_squares) iterates from a starting guess (Gauss-Newton or Levenberg-Marquardt); give a sensible initial guess and check that it converged.
• Report the uncertainty of the parameters (the covariance matrix pcov from curve_fit), not only their values.
• Outliers pull a least-squares line strongly: remove errors you can justify, or use a robust fit (RANSAC, Huber loss).

Linear systems

Many problems reduce to solving A x = b: circuit meshes, statics, finite-element models, calibration. Gaussian elimination reduces the matrix to a triangle and back-substitutes; partial pivoting (swapping in the largest entry) keeps it stable. In practice call the library (numpy.linalg.solve), never invert the matrix.

• Condition number: how much the answer can change for a small change in the data. A matrix with a condition number of 10⁸ loses about eight of the sixteen digits of a double: np.linalg.cond(A). Nearly parallel equations are badly conditioned.
• More equations than unknowns (an over-determined system, as in a fit) are solved in the least-squares sense with np.linalg.lstsq.
• Large sparse systems (finite elements, networks) need sparse matrices and iterative solvers (scipy.sparse.linalg).

Optimization

Finding the minimum of a function (a cost, an error, a loss) is the other big numerical task. Gradient descent steps downhill: x ← x − α f′(x), with a learning rate α. The choice of α matters: too small is slow, too large overshoots and diverges.

• Local minima: a method that goes downhill finds the nearest valley, not necessarily the lowest; try several starting points.
• Use a library: scipy.optimize.minimize offers BFGS, Nelder-Mead (no derivatives) and constrained methods; linear programming has linprog; global search has differential_evolution.
• Constraints and units: scale the variables so that they are of order one, or the optimizer will struggle.

Monte Carlo methods

A Monte Carlo method uses random numbers to estimate a quantity: the average of many random trials. The error falls only as 1 ÷ √N, so ten times the accuracy costs a hundred times the samples, but the method works for problems too complicated for anything else (tolerance stack-ups, reliability, particle transport).

• Seed the random generator (np.random.default_rng(42)) so that results are reproducible, and report the seed.
• Quote the uncertainty: the standard error of the mean falls as s ÷ √N.
• Use independent, well-tested generators (NumPy’s PCG64), not rand() from an old C library.

Use the libraries, and check them

Task	Python (SciPy / NumPy)	MATLAB / Octave	C / C++
Roots	scipy.optimize.root_scalar, fsolve	fzero, fsolve	GSL, a hand-written Newton
Integration	scipy.integrate.quad, np.trapezoid	integral, trapz	GSL, Boost
ODEs	solve_ivp	ode45, ode15s	Boost.odeint, SUNDIALS
Interpolation	np.interp, CubicSpline	interp1, spline	GSL, your own table
Fitting	polyfit, curve_fit	polyfit, lsqcurvefit	Eigen, Ceres
Linear algebra	np.linalg, scipy.linalg	the \ operator	LAPACK, Eigen
Optimization	scipy.optimize.minimize	fminsearch, fmincon	NLopt, Ceres

• Libraries are better than your own code, usually: they handle edge cases and error estimates. Write your own only to learn, or for a microcontroller.
• Check against a case with a known answer (an integral you can do by hand, a circuit you can solve) before you trust a new routine or a new model.
• Change the step or the tolerance and see whether the answer changes in the digits you report: if it does, you have not converged.
• Plot the result. A curve that is the wrong shape is found faster by eye than by any test.

Related

• Programming Languages: the overview of this guide.
• Numbers and Floating Point: the rounding errors behind these methods.
• Filters, PID and Signal Code: filters, PID and the Fourier transform.
• Python for Engineers: NumPy and SciPy in practice.
• Structures and Strength of Materials and Robotics Tools: where these methods are applied.

FAQ

Which numerical integration method should I use?

For smooth functions Simpson’s rule or SciPy’s quad, which adapts itself. For measured samples at equal spacing the trapezoid rule is robust and simple.

Why does my simulation blow up?

Usually the time step is too large for the fastest process (a stability limit), or the model has an error in a unit or a sign. Halve the step and see whether the result changes; try an adaptive or implicit solver.

How small should the step size be?

Small enough that halving it changes the answer by less than the accuracy you need, and not so small that rounding error takes over. Compare two step sizes and look at the difference.

Is least squares the best way to fit data?

It is the standard when the errors are random and similar in size. With outliers or errors in both variables use a robust fit or total least squares; always look at the residuals.

References

• Press, W. H., Teukolsky, S. A., Vetterling, W. T. and Flannery, B. P., Numerical Recipes, Cambridge University Press.
• Burden, R. L. and Faires, J. D., Numerical Analysis; Chapra, S. C. and Canale, R. P., Numerical Methods for Engineers.
• SciPy documentation (optimize, integrate, interpolate, linalg); Higham, N., Accuracy and Stability of Numerical Algorithms.

Last reviewed: October 2026. Error figures are computed in the page builder with double-precision arithmetic.