In the Gulf War, on 25 February 1991, a Patriot missile battery at Dhahran failed to engage an incoming Scud, which struck a barracks and killed 28 American soldiers. The US General Accounting Office traced the failure to a software problem: the system counted time in tenths of a second, and a tenth is not a number a binary computer can store exactly. Each tick carried a tiny error. After about 100 hours of continuous running the clock was off by roughly a third of a second, which for a missile moving at the speed of a Scud is more than half a kilometre.
Nothing about that error was large. What mattered was that nobody had asked how a small error would grow over time. That question is the whole of numerical analysis: how to compute answers with finite digits, how big the error is, and which methods keep it small and which let it explode.
Almost every number you see from a weather forecast, a bridge simulation, a flight path, a game physics engine or an AI model came out of a numerical method. Mathematics gives the equation; numerical analysis gives you an answer you can trust, with an honest estimate of how wrong it is.
Who takes it: Mathematics, Computer Science, Physics and engineering students, usually in the second or third year, once they know calculus, linear algebra and differential equations. Mathematics programmes (MIT's Course 18, Cambridge) treat it as a rigorous subject with proofs; engineering programmes teach it as a toolkit.
What it is for: turning continuous mathematics into finite computation, and knowing when to believe the result.
Comes after: Differential Equations (34/03), Linear Algebra (34/04) and Programming Fundamentals (34/09).
Leads to: Partial Differential Equations (34/79), where every real simulation is numerical; also Machine Learning (34/32), which is numerical optimisation at scale.
Honest line: you will not learn to write a production solver. You will see the six ideas every solver is built around: floating point, the two-error trade-off, fast root-finding, conditioning, order of accuracy and stability. Easier first reads: Differential Equations (34/03), which first met Euler's method breaking down, and Calculus (12/11).
Computers keep about sixteen digits, and 0.1 is not one of them
A computer stores a real number as a sign, an exponent and a fixed-length string of binary digits, in the standard first published as IEEE 754 in 1985 (revised in 2008 and 2019). The common "double" format has 53 significant bits, so neighbouring numbers near 1 differ by 2−52, about 2.2 × 10−16. That gap is called machine epsilon. Between any two stored numbers there is a hole, and every answer is rounded into the nearest one.
Proof sketch: why 0.1 is not stored exactly. A fraction has a finite binary expansion only when its denominator, in lowest terms, is a power of 2. One tenth is 1/(2 × 5), and the 5 is the problem: its binary expansion repeats forever, 0.000110011001100…. The real trick is the denominator test. The computer cuts the repetition after 53 significant bits and keeps the nearest value, which is 0.100000000000000005
Notice that 0.5, which is 2−1, is stored perfectly, and so is any whole number up to 253 = 9,007,199,254,740,992. One more than that cannot be stored: the last button shows it silently rounding back down. The Patriot clock did what the slider does, but with only 24 bits and 3.6 million additions: the tiny drift per tick added up to the third of a second. Rounding errors are small; accumulated or amplified, they are not.
The fix is rarely "use more digits". It is to design the computation so that errors cannot pile up, which the next five ideas show.
Every method has two errors that fight
Suppose you want the slope of a curve at a point and can only evaluate the function. The textbook recipe is the forward difference: (f(x+h) − f(x)) / h. The calculus says the smaller h is, the better. The computer disagrees.
Proof sketch. A Taylor series gives f(x+h) = f(x) + h f′ + ½h²f″ + …, so the recipe is off by about ½h f″: truncation error, from throwing away the rest of the series. But the two function values agree in nearly every digit when h is tiny, so their difference keeps only a few meaningful digits and carries an absolute noise of about εf. Dividing by h magnifies that noise to ε/h. The real trick is to balance the two: ½h ≈ ε/h gives h ≈ √ε ≈ 10−8. Subtracting two nearly equal numbers is called cancellation, and it is the root of most numerical disasters: for example, (1 − cos x)/x² at x = 10−8 returns 0 on a computer, though the true answer is 0.5.
Slide left from large steps. The error falls along a straight line (truncation), reaches a floor, then climbs again as rounding takes over. The forward formula bottoms out near 10−8 with an error of about 3 × 10−9; the central difference, (f(x+h) − f(x−h))/2h, has a smaller truncation error (proportional to h²) and so reaches about 10−11 near h = 10−5. A better formula buys accuracy, but neither can beat the floor.
Newton's method: the number of correct digits doubles
To solve f(x) = 0, draw the tangent to the curve at your current guess, and jump to where the tangent crosses zero. Repeat.
Newton described the idea for polynomials in 1669 (published 1711), Joseph Raphson gave a simpler form in 1690, and Thomas Simpson extended it to general equations with calculus in 1740. For square roots the idea is far older: a Babylonian clay tablet from about 1800–1600 BC, YBC 7289, records √2 as 1;24,51,10 in base 60, about 1.41421296, accurate to six decimals. Historians infer that scribes used an averaging rule equivalent to Newton's, though the tablet shows the answer and not the method.
Proof sketch for √2. Newton's step for x² − 2 is xnext = (x + 2/x)/2, the average of a guess and 2 divided by the guess. A little algebra gives xnext − √2 = (x − √2)² / (2x). The real trick is the square: if the error is 10−3, the next error is about 10−6, and then 10−12. The number of correct digits roughly doubles each step, called quadratic convergence. Starting at 1 the errors are 0.41, 0.086, 0.0025, 0.0000021, 1.6 × 10−12, then rounding level. Bisection, which halves an interval known to bracket the root, gains only one binary digit per step: after the same five steps its error is still about 0.008.
Now try the third equation from a start of 0 or 1. Newton jumps from 0 to 1, back to 0, and never settles: it cycles forever. Try x² − 2 from 0: the slope is zero, so the tangent never crosses the axis. Newton is fast and fragile; bisection is slow and certain. Real solvers use both, bracketing with one and accelerating with the other. A numerical analyst always asks of a method: how fast, and when does it fail?
Some problems amplify errors, whatever method you use
Take two lines that are almost parallel. Their crossing point sits somewhere, but nudge one line by a hair and the crossing slides a long way. That is not the fault of the method used to find the point. It is a property of the problem, measured by its condition number.
Proof sketch. To solve Ax = b, suppose b is off by a small δb. The answer moves by δx = A−1δb. The size of A−1 bounds the output wobble, and the size of A bounds how small b can be relative to x. Putting both together: relative wobble in x ≤ ‖A‖ ‖A−1‖ × relative wobble in b. The real trick is to notice that the inverse of a nearly singular matrix is huge, so tiny input noise (a rounding error is noise of size ε) is blown up by cond(A). A condition number of 10k means you lose about k of your 16 digits, however clever the algorithm.
The famous test case is the Hilbert matrix, with entries 1/(i + j − 1). It looks harmless. Below, the page builds one, sets the true answer to a vector of ones, and solves it in ordinary double precision.
At 5 × 5 the condition number is about 106 and the answer is good to 12 digits. At 10 × 10 it is about 3.5 × 1013 and only three or four digits survive. At 13 × 13 (condition 1018) the "solution" is wrong in its first digit and the program raised no error. That silence is the lesson: a numerical answer does not announce its own accuracy, so you check the condition number, or a residual, or a second method.
The order of a method says what a smaller step buys
To find the area under a curve, cut it into n strips. The trapezoid rule replaces the curve by a straight chord in each strip; Simpson's rule uses a parabola through three points. Both converge as n grows, but at very different speeds.
The exponent p is the method's order. The trapezoid rule has order 2 (halve the strip, error ÷ 4); Simpson's has order 4 (error ÷ 16). On a log–log plot every method is a straight line whose slope is −p, which is how analysts read off the order from an experiment.
Proof sketch. A chord under a smooth curve misses by an amount proportional to the curvature f″ times the strip width squared, and there are 1/h strips, giving a total of order h². The real trick behind Simpson is that the leading h² error terms of two trapezoid estimates can be cancelled: Simpson's rule with 2n strips equals (4T2n − Tn)/3, where T is the trapezoid estimate. Combining two rough answers to cancel their main error is called Richardson extrapolation, and it is how a great deal of accuracy is bought cheaply.
With 8 strips the trapezoid rule is off by about 2 × 10−3 on ∫01 ex dx, and Simpson's by 2 × 10−6. At 1,024 strips Simpson's error hits roundoff and the slope bends flat: the V-shape of the second idea again. An engineer who cares about cost asks, "how many function evaluations does the order buy me?"
A method can be accurate in theory and explode in practice
The 34/03 piece showed Euler's method blowing up when the step is too large. The deeper fact is that the limit on the step often comes from a fast part of the problem you do not even care about. Take y′ = −50(y − cos t) − sin t, started at 1.2. The exact answer is cos t plus a fast transient 0.2e−50t that is gone within 0.1 seconds. After that the solution is a gentle cosine. A problem with a fast, uninteresting scale next to a slow, interesting one is called stiff.
Proof sketch. Explicit Euler multiplies the transient by (1 − 50h) each step, which is smaller than 1 in size only for h < 0.04, even after the transient is gone. Implicit (backward) Euler evaluates the rate at the new time, which gives a multiplier of 1/(1 + 50h), always below 1 for any h. The real trick is to look at what each method does to the pure decaying mode and ask for what steps its multiplier stays below 1. The price: each implicit step means solving an equation, which brings back Newton's method and linear systems from the earlier ideas.
Below h = 0.04 both methods follow the cosine. Push the slider past it and the explicit method flips sign each step, with growing swings: at h = 0.05 it ends at over two million, at 0.1 at over 1011, while the exact answer is about −0.416. The implicit method gets it right to three decimals at the same big steps. The general slogan of the field: a consistent method that is stable converges; without stability, accurate-looking formulas produce garbage.
Where analysts disagree is on how much accuracy is worth buying. Scientific computing usually insists on 64-bit doubles, while the machine-learning world deliberately trains with 16-bit formats that keep only two or three decimal digits, trading precision for speed. Which is right depends on whether the problem is well conditioned and stable, which is exactly what this course teaches you to check.
What a graduate just sees
- A computer keeps about 16 significant digits, and most decimal fractions, like 0.1, are not stored exactly.
- Subtracting two nearly equal numbers throws away digits: avoid cancellation, or rewrite the formula.
- Every method has a truncation error and a rounding error that pull in opposite directions; the best step size is a balance, not "as small as possible".
- Quadratic convergence means correct digits double; linear means one at a time. Always ask how fast and when it fails.
- Some problems are ill-conditioned: the answer is sensitive to the data no matter how good the algorithm is.
- Order p means halving the step divides the error by 2p; a log–log plot shows it as a straight line.
- Stability limits the step; stiff problems force you towards implicit methods.
- A number from a computer does not announce its own error. Check residuals, vary the step, try a second method.
Where it shows up next
- Partial Differential Equations (34/79): every weather, aircraft and reactor simulation is a numerical PDE solver, and the stability limit of the last idea becomes the CFL condition.
- Machine Learning (34/32): training a model is numerical optimisation, and low-precision arithmetic trades the rounding of idea one for speed.
- Linear algebra in practice: libraries such as LAPACK solve linear systems with pivoting precisely because of conditioning.
- Finance and engineering: option pricing, circuit simulators (SPICE) and finite-element structural analysis are all numerical methods with error control.
- Safety-critical software: the Patriot clock is the textbook case of an accumulated rounding error with real consequences.
If you want to go deeper
- 📘 Textbook: Richard Burden, J. Douglas Faires and Annette Burden, Numerical Analysis (10th ed., Cengage, 2015): a widely used first course, algorithm by algorithm.
- 📘 Textbook: Lloyd N. Trefethen and David Bau, Numerical Linear Algebra (SIAM, 1997): a short, elegant treatment of conditioning and stability.
- 📘 Gentler: David Goldberg, "What Every Computer Scientist Should Know About Floating-Point Arithmetic" (ACM Computing Surveys, 1991), free online, and Lloyd Trefethen's short essay "The Definition of Numerical Analysis".
- 🌐 Free course: MIT OpenCourseWare 18.330 Introduction to Numerical Analysis (Spring 2012, Laurent Demanet), with lecture notes.
Computers hold finite digits, so every answer is a rounded approximation. Good numerical work balances truncation and rounding, uses methods that converge fast, notices when the problem itself amplifies errors (conditioning), reads the order of a method from how its error shrinks, and respects stability. The result is an answer plus an honest idea of how wrong it is.
Try this: open a calculator or any programming prompt and compute 1 − (1 − 10−16), and then (1 + 10−16) − 1. Then try 0.1 × 3 and compare it with 0.3. You will see for yourself that arithmetic on a computer is not the arithmetic you learned at school.
Sources. Patriot failure at Dhahran, 25 February 1991: US General Accounting Office, report GAO/IMTEC-92-26 (1992), and the summary by Douglas Arnold (University of Minnesota); the roughly 0.34 s drift after about 100 hours and the figure for the displaced range gate come from secondary write-ups of that report, and are given as approximate. IEEE 754 first published 1985, revised 2008 and 2019. YBC 7289 (about 1800–1600 BC; √2 = 1;24,51,10 ≈ 1.41421296). Newton 1669 (published 1711), Raphson 1690, Simpson 1740: history of Newton's method. Computed in node: stored value of 0.1 = 0.1000000000000000055511…; 0.1 + 0.2 = 0.3000000000000000444; forward-difference and central-difference errors of sin at 1 against step size (minimum 3 × 10−9 near 10−8, and 1.1 × 10−11 near 10−5); (1 − cos x)/x² at x = 10−8 gives 0 against 0.5; Newton errors for √2 from 1; Hilbert matrix 1-norm condition numbers computed from the exact integer inverse formula (27 for n = 2, 9.4 × 105 for n = 5, 3.5 × 1013 for n = 10, 1.3 × 1018 for n = 13) and double-precision solution errors; trapezoid and Simpson errors for ∫01 ex dx; explicit and implicit Euler on the stiff example at several steps. Interactive datasets are illustrative problems chosen for the lesson.
- 🌐 MIT OpenCourseWare 18.330 Introduction to Numerical Analysis
- 📘 Burden, Faires and Burden, Numerical Analysis, 10th ed. (Cengage)
- 📘 Trefethen and Bau, Numerical Linear Algebra (SIAM)
- 🌐 The Patriot Missile Failure (Douglas Arnold)
- 🌐 GAO report IMTEC-92-26