Prev
Unit 1 · Computer Arithmetic and Error
Taylor Series and Truncation Error
60m100 XP
Next
🎯Computer Arithmetic and Error

Taylor Series and Truncation Error

Taylor's theorem with remainder, order of accuracy, the finite-difference sweet spot, and how libraries evaluate functions

Hewlett-Packard, 1972
In 1972 Hewlett-Packard releases the HP-35, the first scientific calculator that fits in a shirt pocket. It does sine, cosine, tangent, logarithms and exponentials — and makes the slide rule obsolete within a few years. Its processor cannot multiply in hardware. It works one decimal digit at a time, adding and shifting. So how does it compute sin⁡(0.7)\sin(0.7)? Not by looking it up: no table would fit. It uses an algorithm called CORDIC, invented by Jack Volder in 1959 for the navigation computer of the B-58 bomber, which turns a trigonometric function into a sequence of additions and shifts. Your laptop does something different — range reduction followed by a short polynomial — but both share one idea. Every elementary function a computer evaluates is an approximation with a known, deliberately accepted error. This lesson is about that deliberate error: where it comes from, how big it is, and how it fights rounding error for control of your answer.
Numerical methods replace something infinite or continuous — a series, a limit, an integral, a differential equation — by something finite. The error that this substitution makes, even in exact arithmetic, is called truncation error (or discretisation error). It is the error of the method. Rounding error, from the last two lessons, is the error of the machine.
The main tool for measuring truncation error is one theorem.
f(x+h)=∑k=0nf(k)(x)k! hk  +  f(n+1)(ξ)(n+1)! hn+1⏟Rn(h),ξ between x and x+hf(x+h) = \sum_{k=0}^{n} \frac{f^{(k)}(x)}{k!}\,h^k \;+\; \underbrace{\frac{f^{(n+1)}(\xi)}{(n+1)!}\,h^{n+1}}_{R_n(h)}, \qquad \xi \text{ between } x \text{ and } x+h

Taylor's theorem with the Lagrange remainder. If f has n+1 continuous derivatives, the error of the degree-n polynomial is exactly R_n for some unknown ξ — so it can be bounded by the largest |f^(n+1)| on the interval.

Derivation
Where Taylor's coefficients come from
Look for the polynomial in hh that behaves like ff near xx:
p(h)=c0+c1h+c2h2+c3h3+⋯+cnhn.p(h) = c_0 + c_1h + c_2h^2 + c_3h^3 + \dots + c_nh^n.
Ask it to match ff and its first nn derivatives at h=0h = 0. Differentiating pp kk times kills every term below hkh^k, turns ckhkc_kh^k into k! ckk!\,c_k, and leaves every higher term still multiplied by a power of hh, which vanishes at h=0h = 0. So
p(k)(0)=k! ck=f(k)(x)⟹ck=f(k)(x)k!.p^{(k)}(0) = k!\,c_k = f^{(k)}(x) \qquad\Longrightarrow\qquad c_k = \frac{f^{(k)}(x)}{k!}.
That is the whole origin of the 1/k!1/k! in every Taylor series.
The remainder needs one more idea. For n=0n = 0 it is the mean value theorem: somewhere between xx and x+hx + h the slope of ff equals the slope of the chord, so f(x+h)=f(x)+hf′(ξ)f(x + h) = f(x) + h f'(\xi). The Lagrange form for general nn is the same statement one level up, proved by applying Rolle's theorem repeatedly to the difference between ff and its Taylor polynomial — the same argument we will use for interpolation in lesson 7-1.
💡Order of accuracy and big-O
We write Rn(h)=O(hn+1)R_n(h) = O(h^{n+1}): as h→0h \to 0, the error shrinks like hn+1h^{n+1} times a constant. This is the single most useful fact about a numerical method. A method of order p has error ≈Chp\approx C h^p. Halve the step and the error falls by 2p2^p: by 2 for first order, by 16 for fourth order. Plot log⁡(error)\log(\text{error}) against log⁡h\log h and the slope is pp — that is how you check a method's order experimentally, and how you catch a bug that silently turns a second-order code into a first-order one.
📐How many terms does sin need?
For sin⁡x\sin x about 0 every derivative is bounded by 1, so the degree-(2n+1)(2n+1) Taylor polynomial has error at most ∣x∣2n+3/(2n+3)!|x|^{2n+3}/(2n+3)!. At x=0.5x = 0.5, eight terms (degree 15) give an error of about 2×10−202 \times 10^{-20} — far beyond double precision. At x=30x = 30 the same polynomial is wrong by about 101010^{10}, and at x=100x = 100 by almost 101810^{18}: the series stops while its terms are still growing, before the factorials catch up. Adding terms until they shrink does not rescue it either. At x=100x = 100 the largest term is about 104210^{42}, and cancelling terms that size down to an answer below 1 destroys every digit. Real libraries therefore reduce the argument first. Using sin⁡(x+kπ/2)\sin(x + k\pi/2) symmetries they map any xx to rr with ∣r∣≤π/4|r| \le \pi/4, then evaluate a short polynomial in rr. The lab builds this: range reduction plus an 8-term polynomial matches math.sin to about 10−1410^{-14} on [−100,100][-100, 100].
Far better an approximate answer to the right question, which is often vague, than an exact answer to the wrong question, which can always be made precise.
– John Tukey, The future of data analysis (1962)
The derivative and its sweet spot. Taylor's theorem gives the simplest numerical derivative directly. Rearranging f(x+h)=f(x)+hf′(x)+h22f′′(ξ)f(x+h) = f(x) + h f'(x) + \tfrac{h^2}{2} f''(\xi):
f′(x)≈f(x+h)−f(x)h    error≈h2∣f′′∣,f′(x)≈f(x+h)−f(x−h)2h    error≈h26∣f′′′∣f'(x) \approx \frac{f(x+h) - f(x)}{h} \;\; \text{error} \approx \frac{h}{2}|f''|, \qquad f'(x) \approx \frac{f(x+h) - f(x-h)}{2h} \;\; \text{error} \approx \frac{h^2}{6}|f'''|

Forward difference (first order) and central difference (second order): the odd terms of the two Taylor expansions cancel in the central formula.

Derivation
Both difference formulas from Taylor's theorem
Expand ff a step forward and a step back:
f(x+h)=f+hf′+h22f′′+h36f′′′+…f(x + h) = f + hf' + \tfrac{h^2}{2}f'' + \tfrac{h^3}{6}f''' + \dots
f(x−h)=f−hf′+h22f′′−h36f′′′+…f(x - h) = f - hf' + \tfrac{h^2}{2}f'' - \tfrac{h^3}{6}f''' + \dots
Forward difference. Rearrange the first line alone:
f(x+h)−f(x)h=f′+h2f′′+…\frac{f(x + h) - f(x)}{h} = f' + \tfrac{h}{2}f'' + \dots
The leading error is h2f′′\tfrac{h}{2}f'': proportional to hh, first order.
Central difference. Subtract the second line from the first. The ff terms cancel, the f′′f'' terms cancel, and the odd terms double:
f(x+h)−f(x−h)=2hf′+h33f′′′+…⟹f(x+h)−f(x−h)2h=f′+h26f′′′+…f(x + h) - f(x - h) = 2hf' + \tfrac{h^3}{3}f''' + \dots \qquad\Longrightarrow\qquad \frac{f(x + h) - f(x - h)}{2h} = f' + \tfrac{h^2}{6}f''' + \dots
The leading error is h26f′′′\tfrac{h^2}{6}f''': second order. Symmetry bought a whole order of accuracy for the same two function evaluations.
Truncation error says: make hh small. But the numerator subtracts two nearly equal numbers — the cancellation of the last lesson. Each value of ff carries a rounding error of about u∣f∣u|f|, and dividing their difference by hh amplifies it to about 2u∣f∣/h2u|f|/h. The total error of the forward difference is therefore roughly
E(h)≈h2∣f′′∣+2u∣f∣hE(h) \approx \dfrac{h}{2}|f''| + \dfrac{2u|f|}{h}
Truncation falls with hh; rounding grows as hh shrinks. The minimum is near h≈u≈10−8h \approx \sqrt{u} \approx 10^{-8}, with an error around 10−810^{-8} — half the digits of double precision are gone. The central difference does better, with an optimum near u1/3≈10−5u^{1/3} \approx 10^{-5} and an error near 10−1110^{-11}. Push hh to 10−1610^{-16} and the answer is garbage.
Derivation
The best step, found by calculus
Write M=∣f′′∣M = |f''|, T=∣f′′′∣T = |f'''| and F=∣f∣F = |f| near xx.
Forward difference. Truncation contributes h2M\tfrac{h}{2}M. Each of f(x+h)f(x + h) and f(x)f(x) carries a rounding error up to about uFuF, so their difference carries up to 2uF2uF, and dividing by hh gives 2uF/h2uF/h. Minimise the total:
E(h)=hM2+2uFh,E′(h)=M2−2uFh2=0  ⟹  h∗=2uFM.E(h) = \frac{hM}{2} + \frac{2uF}{h}, \qquad E'(h) = \frac{M}{2} - \frac{2uF}{h^2} = 0 \;\Longrightarrow\; h^* = 2\sqrt{\frac{uF}{M}}.
At h∗h^* the two parts are equal and E(h∗)=2uFME(h^*) = 2\sqrt{uFM}. For a function with FF and MM of order 1, h∗≈2u≈2×10−8h^* \approx 2\sqrt{u} \approx 2 \times 10^{-8} and E≈2×10−8E \approx 2 \times 10^{-8}: half of the sixteen digits.
Central difference. Truncation is h26T\tfrac{h^2}{6}T; the rounding error 2uF2uF is now divided by 2h2h:
E(h)=h2T6+uFh,E′(h)=hT3−uFh2=0  ⟹  h∗=(3uFT)1/3.E(h) = \frac{h^2T}{6} + \frac{uF}{h}, \qquad E'(h) = \frac{hT}{3} - \frac{uF}{h^2} = 0 \;\Longrightarrow\; h^* = \Big(\frac{3uF}{T}\Big)^{1/3}.
With F,T∼1F, T \sim 1: h∗≈(3.3×10−16)1/3≈7×10−6h^* \approx (3.3 \times 10^{-16})^{1/3} \approx 7 \times 10^{-6}, and the error there is 12h∗2T≈2×10−11\tfrac12 h^{*2}T \approx 2 \times 10^{-11} — about eleven digits. The general lesson: a method of order pp balanced against rounding has its best step near u1/(p+1)u^{1/(p+1)} and loses about 1/(p+1)1/(p+1) of the digits.
🕰A derivative with no cancellation
In 1967 Lyness and Moler found a trick that sidesteps the V entirely when ff can be evaluated for complex arguments. Expand f(x+ih)f(x + ih): the real part carries f(x)f(x), the imaginary part carries hf′(x)h f'(x). So f′(x)≈Im⁡f(x+ih)/hf'(x) \approx \operatorname{Im} f(x+ih)/h involves no subtraction, and h=10−100h = 10^{-100} gives a derivative accurate to the last bit. It is now used to verify aerodynamic and structural sensitivities in large engineering codes.
Derivation
Why the complex step works
Taylor's theorem holds for a step in any direction of the complex plane, provided ff is smooth (analytic) there. Step by ihih, an imaginary amount, and use i2=−1i^2 = -1, i3=−ii^3 = -i:
f(x+ih)=f(x)+ihf′(x)−h22f′′(x)−i h36f′′′(x)+…f(x + ih) = f(x) + ih f'(x) - \frac{h^2}{2}f''(x) - i\,\frac{h^3}{6}f'''(x) + \dots
For real xx and a function that is real on the real line, the real and imaginary parts separate cleanly:
Re⁡f(x+ih)=f−h22f′′+…,Im⁡f(x+ih)=hf′−h36f′′′+…\operatorname{Re} f(x + ih) = f - \tfrac{h^2}{2}f'' + \dots, \qquad \operatorname{Im} f(x + ih) = hf' - \tfrac{h^3}{6}f''' + \dots
Divide the imaginary part by hh:
Im⁡f(x+ih)h=f′(x)−h26f′′′(x)+…\frac{\operatorname{Im} f(x + ih)}{h} = f'(x) - \frac{h^2}{6}f'''(x) + \dots
Second order, like the central difference — but nothing was subtracted. The only rounding error is the one made in computing f(x+ih)f(x + ih) itself, so hh can be made as small as we like: at h=10−100h = 10^{-100} the truncation term is 10−20010^{-200} and the derivative is correct to the last bit.
💡Where you will use this
Every time step of an ODE solver in Orbital Mechanics, Attitude Dynamics and Rocket trajectories is a truncated Taylor series; its order decides your step size. Jacobians for Newton's method in Rocket Propulsion and Power Systems are often finite differences, and choosing hh is exactly this lesson. Finite-difference and finite-element meshes in Structures, Heat Transfer and Geotechnical Engineering are described by their order of accuracy.
📚Further reading
  • Michael Heath, Scientific Computing: An Introductory Survey (2nd ed.; SIAM Classics 2018), chapter 1 (truncation and rounding error) and the numerical-differentiation sections of chapter 8.
  • Jean-Michel Muller, Elementary Functions: Algorithms and Implementation (3rd ed., Birkhäuser, 2016) — range reduction, polynomial approximation and CORDIC, as libraries really do them.
  • William Squire and George Trapp, 'Using complex variables to estimate derivatives of real functions', SIAM Review 40(1), 110–112 (1998) — the three-page paper that revived the complex step.
↓But there is a problem
We now know where error comes from: rounding in every operation, truncation in every method. Good algorithms keep both small. But there is a third ingredient, and no algorithm can fix it. In 1961 a meteorologist at MIT re-ran a weather simulation from the middle, typing in a number from a printout: 0.506 instead of the 0.506127 the machine had held. A difference in the fourth decimal. Within two months of simulated weather the two runs had nothing in common. The machine had made no mistake, and neither had his algorithm. The problem itself was amplifying the error.