Atlas / Mathematics / The Computation Thread

Field · Emerged 1922 – 1977

Numerical Methods for PDEs

How do you turn the continuous equations of physics into a finite computation whose answer converges to the truth?

5 chapters5 min read6 turning points 1 open problem

Branched from: Numerical Analysis + Differential Equations
Branched into: Not yet surveyed past here
Figures: Lewis Fry Richardson, Richard Courant, Kurt Friedrichs, Hans Lewy, John von Neumann, Jule Charney, Ray Clough, Radii Fedorenko, Achi Brandt, Frans Pretorius

In brief

The laws of fluids, heat, elasticity, electromagnetism and gravity are partial differential equations: they relate how a quantity changes in space to how it changes in time. Almost none can be solved by formula for a realistic shape or a realistic weather map. Numerical methods replace the continuous equation by a finite one, on a grid of points or a mesh of small pieces, and solve that instead.

The first serious attempt, Lewis Fry Richardson's hand-computed weather forecast of 1922, failed badly. Six years later Courant, Friedrichs and Lewy found one reason such computations fail: the time step must be small enough for the grid to keep up with the physics. The ENIAC weather forecast of 1950 and the finite element method, which grew from engineering in the 1950s and 1960s, made the approach practical. Today it designs aircraft, forecasts the weather and predicted the gravitational waves that LIGO detected in 2015.

Key ideas

DiscretisationEnters 1922: Replace a function by its values at finitely many points, and derivatives by differences between neighbouring values. The differential equation becomes a large system of ordinary equations.
The CFL conditionEnters 1928: For an explicit time-stepping method, the time step must be short enough that information does not travel more than about one grid cell per step. Otherwise the computation becomes unstable and explodes.
Stability and convergenceEnters 1928: A method whose small errors stay small is stable. Lax's equivalence theorem says that, for a consistent method on a linear problem, stability is exactly what is needed for the computed answer to converge to the true one as the grid is refined.
Finite elementsEnters 1943 – 1960: Divide the region into small simple pieces, such as triangles, and approximate the solution on each by a simple polynomial, matched at the corners. The method fits complicated shapes, which grids of squares do not.
MultigridEnters 1964 – 1977: Solve on a hierarchy of grids, coarse and fine. Errors that are smooth on a fine grid are rough on a coarse one, where they can be removed cheaply. The total work grows only in proportion to the number of unknowns.

Draws on other domains

↙ Physics
Elasticity and Continuum Mechanics
The finite element method came from this problem

Chapter I

The Forecast Factory

Differential equations describe how a state changes from moment to moment. For the atmosphere, the state is the wind, pressure, temperature and moisture at every point, and the equations are the partial differential equations of fluid motion. In principle, today's weather determines tomorrow's. In practice nobody could solve the equations.

Lewis Fry Richardson decided to compute instead. He divided the atmosphere into boxes, replaced derivatives by differences between neighbouring boxes, and stepped the equations forward in time by arithmetic. During the First World War, between shifts with an ambulance unit in France, he computed a six-hour forecast for two points in central Europe. It took about six weeks. The result was a change in pressure of 145 hectopascals in six hours, when the real pressure barely changed. He published the whole calculation in 1922, with the failure stated plainly, and imagined a hall of 64,000 people computing the world's weather in time. Much later, Peter Lynch recomputed the forecast and showed that the method was sound. The starting data were unbalanced, and a small smoothing would have given a sensible answer.

Chapter II

Stability

In 1928 Richard Courant, Kurt Friedrichs and Hans Lewy studied grid equations for a different reason: to prove that certain partial differential equations have solutions, by showing that grid solutions converge as the grid is refined. They found that for wave-like equations this happens only if the time step is small enough compared with the grid spacing. If a wave can cross more than one cell in a single step, the grid cannot follow it.

The condition became essential when electronic computers arrived. In 1950 Jule Charney, working with John von Neumann and Ragnar Fjørtoft, ran the first computer weather forecast on the ENIAC. Charney had simplified the equations to remove the fast gravity waves that had wrecked Richardson's forecast, and chose time steps that respected the CFL condition. Each 24-hour forecast took about 24 hours to compute, but the results were realistic. Von Neumann also developed a way to test a scheme's stability by following each wave pattern separately, and in 1956 Peter Lax and Robert Richtmyer proved that, for linear problems, a consistent scheme converges exactly when it is stable.

Chapter III

Elements and Grids

Grids of squares suit the atmosphere but not an aircraft wing. In the 1950s engineers, among them Ray Clough, working with a Boeing team, and John Argyris in London, broke structures into small triangular and rectangular elements, each with simple behaviour, joined at their corners. Clough named it the finite element method in 1960. Courant had proposed the same idea in 1943, as a way to solve variational problems, and mathematicians later proved that the engineers' method converges and how fast. It is now the standard tool for anything with a complicated shape.

Solving the resulting equations was the next bottleneck. A fine three-dimensional mesh has millions of unknowns, and simple iterative methods need many thousands of sweeps to converge, because smooth errors fade slowly. Radii Fedorenko saw in 1964 that a coarser grid could remove those smooth errors cheaply, and Achi Brandt turned the idea into the general multigrid method in 1977. Together with the linear algebra of conjugate gradients, it made problems with billions of unknowns solvable.

Chapter IV

A Closer Look: The CFL Condition

The simplest wave-like equation says that a pattern $u$ moves to the right at speed $c$ without changing shape:

\frac{\partial u}{\partial t} + c\,\frac{\partial u}{\partial x} = 0.

Put it on a grid with spacing $\Delta x$ and time step $\Delta t$ . The upwind scheme updates each value from itself and its left-hand neighbour:

u_j^{\text{new}} = u_j - \nu\,(u_j - u_{j-1}), \qquad \nu = \frac{c\,\Delta t}{\Delta x}.

The number $\nu$ , the Courant number, is the number of cells the pattern moves in one step. If $\nu \le 1$ , the new value is $(1 - \nu)u_j + \nu u_{j-1}$ , a weighted average of two old values, so it can never exceed the largest of them. If $\nu > 1$ , the weight $1 - \nu$ is negative and nothing holds the values in check.

To see what happens, add a tiny sawtooth ripple of size 0.001 that alternates in sign from cell to cell, like rounding error. For this ripple $u_{j-1} = -u_j$ , so each step multiplies it by $1 - 2\nu$ :

Courant number $\nu$	Factor per step	Ripple after 20 steps
0.5	0	0
1.0	−1	0.001
1.5	−2	about 1,000
2.0	−3	about 3.5 million

At $\nu = 2$ the ripple, started at one thousandth, is 3.5 million after 20 steps. No real wave did that. The computation has exploded.

Now take a weather model with a 10-kilometre grid and winds of up to 100 metres per second. The condition $\nu \le 1$ requires $\Delta t \le 10{,}000 / 100 = 100$ seconds, so at least 864 steps for a one-day forecast. Sound waves, at about 340 metres per second, would force steps under 30 seconds. That is why Charney filtered out the fast waves, and why modern models treat them with implicit methods that are not bound by the condition.

Chapter V

From Grids to Spacetime

The same methods now simulate almost every physical system. The hardest test came from general relativity. For decades simulations of two orbiting black holes crashed after a fraction of an orbit, as small errors in Einstein's equations grew without limit. The instability was not in the grid but in the way the equations were written. In 2005 Frans Pretorius found a formulation that stayed stable and followed two black holes through their merger. Ten years later, LIGO's first detection was matched against such simulations. The underlying mathematics is still incomplete, and for the equations of gas dynamics in more than one dimension it is not known whether the simulations converge at all.

Turning points6

1922CRISIS
Richardson's forecast by hand
Lewis Fry Richardson, working between shifts as an ambulance driver in the First World War, computes a six-hour forecast for two points over central Europe from the equations of fluid motion. The calculation takes him about six weeks. It predicts a change in surface pressure of 145 hectopascals in six hours, when the observed pressure barely changed. He publishes the method anyway, together with a dream of a "forecast factory" of 64,000 human computers.
Lewis Fry Richardson
› Sources (2)
- Richardson, L. F. (1922). Weather Prediction by Numerical Process. Cambridge University Press.
- Lynch, P. (2006). The Emergence of Numerical Weather Prediction: Richardson's Dream. Cambridge University Press.
1928PROOF
The Courant–Friedrichs–Lewy condition
Richard Courant, Kurt Friedrichs and Hans Lewy, in Göttingen, use difference equations on grids to prove existence theorems for partial differential equations. Along the way they find that for wave-like equations the grid solution converges only if the time step is small enough relative to the spacing. The condition was a side result of pure mathematics. After 1945 it became the first rule of scientific computing.
Richard Courant, Kurt Friedrichs, Hans Lewy
› Sources (1)
- Courant, R., Friedrichs, K. & Lewy, H. (1928). Über die partiellen Differenzengleichungen der mathematischen Physik. Mathematische Annalen 100: 32–74.
1950REFORMULATION
The first computer forecast
Jule Charney, Ragnar Fjørtoft and John von Neumann run a simplified model of the atmosphere on the ENIAC in Aberdeen, Maryland. Charney had filtered out the fast gravity waves that had ruined Richardson's attempt. Four 24-hour forecasts are produced, each taking about a day of machine time, and some are realistic. Routine computer forecasting began in Sweden and the United States within five years.
John von Neumann, Jule Charney
› Sources (1)
- Charney, J. G., Fjørtoft, R. & von Neumann, J. (1950). Numerical integration of the barotropic vorticity equation. Tellus 2(4): 237–254.
1943 – 1960REFORMULATIONContested
The finite element method
Richard Courant in 1943 approximates a twisting problem by simple functions on triangles. In the 1950s aircraft engineers, among them John Argyris in London and a Boeing team including Ray Clough, independently break wings into small elements to compute stresses. Clough names the method in 1960. Mathematicians later showed it is a systematic way of solving variational problems, with error bounds, and it became the standard tool for engineering analysis.
What is disputed
Who invented the finite element method depends on what counts. Courant's 1943 paper contains the mathematical idea but was not followed up. The engineers of the 1950s developed the practical method without knowing of it. Argyris, Clough, Turner and Zienkiewicz have all been called its founder, and the balance of credit between the mathematical and engineering lines is still discussed.
Richard Courant, Ray Clough
› Sources (2)
- Courant, R. (1943). Variational methods for the solution of problems of equilibrium and vibrations. Bulletin of the American Mathematical Society 49(1): 1–23.
- Turner, M. J., Clough, R. W., Martin, H. C. & Topp, L. J. (1956). Stiffness and deflection analysis of complex structures. Journal of the Aeronautical Sciences 23(9): 805–823.
1964 – 1977REFORMULATION
Multigrid
Standard iterative methods remove rough errors quickly but smooth errors very slowly. Radii Fedorenko in Moscow shows in 1964 that switching to a coarser grid, where smooth errors look rough, fixes this. Achi Brandt develops the idea into a general, practical method in 1977. Multigrid solves the equations of a grid with $N$ points in a number of operations proportional to $N$ , the best possible.
Radii Fedorenko, Achi Brandt
› Sources (2)
- Fedorenko, R. P. (1964). The speed of convergence of one iterative process. USSR Computational Mathematics and Mathematical Physics 4(3): 227–235.
- Brandt, A. (1977). Multi-level adaptive solutions to boundary-value problems. Mathematics of Computation 31(138): 333–390.
2005 – 2006REFORMULATION
Two black holes merge on a computer
For forty years, computer simulations of two orbiting black holes crashed before the holes could merge. The equations of general relativity have many formulations, and most let small errors grow. In 2005 Frans Pretorius, using a reformulation based on harmonic coordinates, simulates a full orbit, merger and ringdown. Within months two other groups succeeded with a different method, "moving punctures".
Frans Pretorius
› Sources (2)
- Pretorius, F. (2005). Evolution of binary black-hole spacetimes. Physical Review Letters 95(12): 121101.
- Campanelli, M., Lousto, C. O., Marronetti, P. & Zlochower, Y. (2006). Accurate evolutions of orbiting black-hole binaries without excision. Physical Review Letters 96(11): 111101.

Applications

Where it is used

Weather and climate↗ Physics · Geophysical Fluid Dynamics
Numerical weather prediction
Every modern forecast solves the equations of fluid motion on a grid covering the globe, with spacing of about ten kilometres. Forecast skill has improved by about a day per decade: a forecast for six days ahead is now as good as a forecast for five days ahead was ten years earlier.
› Sources (1)
- Bauer, P., Thorpe, A. & Brunet, G. (2015). The quiet revolution of numerical weather prediction. Nature 525: 47–55.
Gravitational waves↗ Physics · General Relativity
Templates for LIGO
The signal LIGO detected in September 2015 was identified by comparing it with waveforms of merging black holes computed by numerical relativity. The match gave the masses of the two holes, about 36 and 29 times the Sun's.
› Sources (1)
- Abbott, B. P. et al. (2016). Observation of gravitational waves from a binary black hole merger. Physical Review Letters 116(6): 061102.
Engineering↗ Physics · Elasticity and Continuum Mechanics
Designing on the computer
Aircraft, bridges, engines and car bodies are tested by finite element simulation long before anything is built. Crash tests, for example, are now mostly run on computers, with physical tests to confirm the result.
› Sources (1)
- Zienkiewicz, O. C., Taylor, R. L. & Zhu, J. Z. (2013). The Finite Element Method: Its Basis and Fundamentals, 7th edition. Butterworth-Heinemann.

Open problems

Where the map runs out

Open

Convergence for shock waves in several dimensions

Open as of 2026. Largely settled in one space dimension for data with small total variation, open in two and three.

The equations of gas dynamics form shock waves, and their solutions must be understood in a weak sense. In one space dimension, for data with small total variation, some numerical schemes are known to converge to the right solution. In two or three dimensions nobody knows whether they do, or even whether the equations have a unique physically admissible solution for general data.

Why it is hard: The one-dimensional theory relies on tools that do not extend to higher dimensions. Worse, examples have been found of initial data for the compressible Euler equations with infinitely many admissible weak solutions. If the equations do not pick out one answer, it is unclear what a numerical scheme should converge to.
What resolving it unlocks: A theory would say when simulations of explosions, supersonic flight and supernovae, run every day, are converging to a true solution and when they are only producing plausible pictures.

› Sources (2)

Dafermos, C. M. (2016). Hyperbolic Conservation Laws in Continuum Physics, 4th edition. Springer.
Chiodaroli, E., De Lellis, C. & Kreml, O. (2015). Global ill-posedness of the isentropic system of gas dynamics. Communications on Pure and Applied Mathematics 68(7): 1157–1190.