Applied numerical computation
The computer does not calculate with real numbers but with approximations, and that difference has consequences. Numerical methods solve what has no analytical solution, always with an error that you must know how to bound.
01The computer does not calculate with real numbers
A double stores 64 bits: one sign bit, 11 of exponent and 52 of mantissa. With that it represents
about 16 significant decimal digits and nothing more. Everything else is rounded.
- 0.1 + 0.2 does not give 0.3. It gives 0.30000000000000004, because 0.1 does not exist exactly in binary.
Floats are never compared with
==: you compare \( |a-b| \lt \varepsilon \). - Catastrophic cancellation. Subtracting two similar numbers destroys the significant digits: if each has 16 good digits and they agree in 15, the subtraction leaves just one.
- Order matters. Adding a million small numbers to a large one loses the small ones; adding from smallest to largest preserves more.
The machine epsilon is the smallest number that, added to 1, gives something different from 1: about \( 2.2 \times 10^{-16} \) in double precision. It sets the floor of the relative error of any calculation.
02Error, convergence and stopping criterion
| Type of error | Where it comes from |
|---|---|
| Round-off | From the machine’s finite arithmetic |
| Truncation | From cutting off an infinite process: a series, a limit |
| In the data | From the measurement that went into the calculation |
A method converges linearly if each step multiplies the error by a constant, and quadratically if it squares it: in that case the correct digits double with each iteration. Since the exact answer is almost never known, the practical stopping criterion combines three things: the change between iterations is small, the residual \( |f(x)| \) is small, and an iteration limit so as not to keep spinning if it does not converge.
Choose the function and the method, and advance iteration by iteration watching the error fall.
03Roots of an equation
| Method | Iteration | Convergence | Requires |
|---|---|---|---|
| Bisection | \( x_{k+1} = \dfrac{a_k + b_k}{2} \) | Linear, one bit per step | A sign change |
| Newton-Raphson | \( x_{k+1} = x_k - \dfrac{f(x_k)}{f'(x_k)} \) | Quadratic | The derivative |
| Secant | \( x_{k+1} = x_k - f(x_k)\dfrac{x_k - x_{k-1}}{f(x_k) - f(x_{k-1})} \) | Order 1.618 | Two starting points |
| Fixed point | \( x_{k+1} = g(x_k) \) | Linear, if \( |g'| \lt 1 \) | Rewriting the equation |
Bisection never fails if there is a sign change, but it is slow: each step gains one bit, so
about 50 iterations are needed to exhaust a double. Newton flies —it doubles the correct
digits at each step— but it can go anywhere if the derivative is small or the starting point is
far away. In practice they are combined: bisection to get close, Newton to finish.
04The other three calculations that always show up
Numerical integration
Least-squares fitting
For a line \( y = mx + b \) that passes as close as possible to \( n \) measured points:
It is the calculation behind every sensor calibration: pairs are measured and the line is fitted that converts converter counts into physical units.
Systems of equations
For a few nodes, Gaussian elimination with partial pivoting is enough —choosing as pivot the largest element of the column— which avoids dividing by small numbers. For large, sparse systems, such as those that come from a resistor network, the iterative Jacobi and Gauss-Seidel methods are the better choice.
05Where it appears in electronics
- Calibrating a sensor. Least squares on the measured pairs.
- Operating point of a nonlinear circuit. Solving for the intersection of the load line with the diode curve is exactly the problem of this page.
- Simulation. SPICE solves a nonlinear system at each step with Newton-Raphson, and advances in time by integrating numerically.
- Digital filtering. Moving averages, median filters and accumulators, taking care not to lose precision through the order of the sums.
- Microcontrollers without a floating-point unit. There you work in fixed point: faster, but the scale has to be managed by hand.
06In the lab
Print 0.1 + 0.2 == 0.3, and then the value with twenty decimals. Calculate the machine
epsilon with a loop that divides by two until 1 + e == 1.
Solve \( x^2 = 2 \) by bisection and by Newton, tabulating the error at each iteration. Plot the logarithm of the error: bisection gives a straight line and Newton, a curve that plunges.
With a temperature sensor and a reference thermometer, measure ten pairs of values, fit the line by least squares and program that conversion into the microcontroller. Compare the error before and after calibrating.
07Common mistakes
- Comparing floats with
==. - Subtracting similar quantities without realizing that digits are lost.
- Iterating without a limit and ending up in an infinite loop when the method does not converge.
- Using Newton with an arbitrary starting point and ending up at another root, or diverging.
- Confusing precision with accuracy: printing fifteen decimals of a quantity measured with two significant figures.
- Not checking the result by substituting it into the original equation.
08Self-assessment
Why does 0.1 + 0.2 not give exactly 0.3?
Because 0.1 and 0.2 have no finite binary representation: they are stored rounded and the error shows up when adding.
How many bisection iterations are needed to gain three decimals?
Each iteration divides the interval by two, so it gains one bit. For a factor of 1000 about 10 iterations are needed, because \( 2^{10} = 1024 \).
What does it mean that Newton converges quadratically?
That the new error is proportional to the square of the previous one: the correct digits double at each step. From 2 to 4, from 4 to 8, from 8 to 16.
When does Newton fail?
When the derivative is small or zero near the point, when the starting point is far away, or when the function oscillates: it can diverge or jump to another root.
What error does Simpson’s rule have compared with the trapezoidal rule?
\( O(h^4) \) versus \( O(h^2) \): when the step is halved, Simpson’s error falls by a factor of 16 and the trapezoidal rule’s by 4.
What is partial pivoting for in Gaussian elimination?
To avoid dividing by a very small number, which would amplify the round-off error. The largest element of the column is chosen as the pivot.
09Further reading
- Richard L. Burden and J. Douglas Faires. Numerical Analysis. 10th ed. (Spanish edition, “Análisis numérico”), Cengage Learning, 2016. The standard text for the subject, with algorithms and error analysis.
- Steven C. Chapra and Raymond P. Canale. Numerical Methods for Engineers. 7th ed. (Spanish edition, “Métodos numéricos para ingenieros”), McGraw-Hill, 2015. Strongly application-oriented, with code and engineering cases.
- William H. Press et al. Numerical Recipes. 3rd ed., Cambridge University Press, 2007. Real implementations, with honest comments on when each method fails.
- David Goldberg. What Every Computer Scientist Should Know About Floating-Point Arithmetic. ACM Computing Surveys, 1991. The article to read once in your life on floating point.