Appearance
10.P — Worked Problems: Numerical Methods & Optimisation
Ten problems, each run iteration by iteration so you can watch the digits arrive — and, in two cases, watch them leave. Every formula used is derived in the chapter that introduces it.
Problem 1 — Floating point, and a subtraction that destroys an answer
Compute \sqrt{10^{10}+1}-\sqrt{10^{10}} two ways and explain the difference.
Solution
The direct route.
\sqrt{10^{10}+1} = \sqrt{10000000001} = 100000.000005
\sqrt{10^{10}} = 100000
A double holds about 16 significant digits. 100000.000005 needs 12 digits to be written, so it is stored fine, but only the last few digits carry the information we want. Subtracting:
100000.000005 - 100000 = 0.000005
We started with two numbers known to 16 digits and ended with an answer known to about 5 digits. The leading 11 digits cancelled and took their precision with them.
The stable route. Multiply top and bottom by the conjugate:
\sqrt{A+1}-\sqrt A = \frac{\left(\sqrt{A+1}-\sqrt A\right)\left(\sqrt{A+1}+\sqrt A\right)}{\sqrt{A+1}+\sqrt A} = \frac{(A+1)-A}{\sqrt{A+1}+\sqrt A} = \frac{1}{\sqrt{A+1}+\sqrt A}
= \frac{1}{100000.000005+100000} = \frac{1}{200000.000005} = 4.99999999988\times10^{-6}
Now there is no subtraction of nearly equal quantities, so all 16 digits survive.
Answers: 5.0×10⁻⁶ from the naive route with 5 good digits; 4.99999999988×10⁻⁶ from the stable route with 16.
The general lesson. The problem itself is perfectly well conditioned — a small change in A produces a small change in the answer. It was the algorithm that was unstable. Whenever a formula subtracts two nearly equal things, look for an algebraic rewrite that removes the subtraction. The conjugate trick, the half-angle identity, and \ln(1+x) in place of \ln(1+x)-\ln(1) are the three that come up most.
Problem 2 — Bisection, step by step
Find the root of f(x) = x^3-x-2 in [1,2] to within 0.05 by bisection.
Solution
Step 0 — check the bracket.
f(1) = 1-1-2 = -2 \qquad f(2) = 8-2-2 = +4
Opposite signs, so a root lies between. The method cannot fail from here.
Step 1 — how many iterations will be needed?
\frac{2-1}{2^n}<0.05 \quad \Rightarrow \quad 2^n>20 \quad \Rightarrow \quad n\ge5
Knowing the count in advance is a property only bisection has.
Step 2 — iterate. At each step, take the midpoint, evaluate, and keep the half that still straddles zero.
| n | a | b | m | f(m) | keep |
|---|---|---|---|---|---|
| 1 | 1 | 2 | 1.5 | -0.125 | right half |
| 2 | 1.5 | 2 | 1.75 | +1.609 | left half |
| 3 | 1.5 | 1.75 | 1.625 | +0.666 | left half |
| 4 | 1.5 | 1.625 | 1.5625 | +0.252 | left half |
| 5 | 1.5 | 1.5625 | 1.53125 | +0.059 | left half |
Working for the first row: f(1.5) = 3.375-1.5-2 = -0.125, which is negative like f(1), so the root is between 1.5 and 2.
Step 3 — report. After 5 steps the bracket is [1.5, 1.53125], width 0.03125 < 0.05 ✓
Answer: x \approx 1.516, with the root guaranteed to lie in [1.5, 1.53125].
The true root is 1.521380. Our midpoint estimate 1.5156 is out by 0.0058, comfortably inside the promised bound.
What makes bisection worth keeping. It is slow — one bit per step, so 16 digits would need 53 iterations — but it cannot diverge, cannot cycle, and gives a guaranteed interval rather than an estimate. Every serious root-finder uses it as a fallback when a faster method misbehaves.
Problem 3 — Newton's method, and watching digits double
Find \sqrt{2} using Newton's method on f(x)=x^2-2, starting from x_0 = 1.
Solution
Step 1 — derive the iteration. With f(x) = x^2-2 and f'(x) = 2x:
x_{n+1} = x_n - \frac{x_n^2-2}{2x_n}
Simplify by putting it over a common denominator:
= \frac{2x_n^2 - x_n^2+2}{2x_n} = \frac{x_n^2+2}{2x_n} = \frac12\left(x_n+\frac{2}{x_n}\right)
That final form is worth recognising: it is the average of your guess and 2 divided by your guess. If the guess is too small, \frac2x is too big, and the average lands between — this is the Babylonian method, in use around 1700 BC, and Newton's method rediscovers it exactly.
Step 2 — iterate.
x_1 = \tfrac12\left(1+2\right) = 1.5
x_2 = \tfrac12\left(1.5+\tfrac{2}{1.5}\right) = \tfrac12(1.5+1.333333) = 1.4166667
x_3 = \tfrac12\left(1.4166667+1.4117647\right) = 1.4142157
x_4 = \tfrac12\left(1.4142157+1.4142114\right) = 1.4142136
Step 3 — count the correct digits at each step. The true value is 1.41421356237.
| n | x_n | error | correct digits |
|---|---|---|---|
| 0 | 1 | 4\times10^{-1} | 0 |
| 1 | 1.5 | 8.6\times10^{-2} | 1 |
| 2 | 1.4166667 | 2.5\times10^{-3} | 2 |
| 3 | 1.4142157 | 2.1\times10^{-6} | 5 |
| 4 | 1.4142136 | 1.6\times10^{-12} | 11 |
Answer: after four iterations, \sqrt2 = 1.41421356, correct to every digit shown.
The doubling is visible in the error column: 10^{-3}\to10^{-6}\to10^{-12}. That is the quadratic convergence law e_{n+1}\approx Ce_n^2 with C = \frac{f''}{2f'} = \frac{2}{2(1.414)} = 0.354.
Check the prediction. From e_3 = 2.1\times10^{-6}, the law predicts
e_4 \approx 0.354\times(2.1\times10^{-6})^2 = 1.56\times10^{-12}
and the actual error was 1.6\times10^{-12} ✓
One more iteration would exceed double precision entirely, which is why library square-root routines run about four or five Newton steps from a good table lookup and stop.
Problem 4 — Where Newton fails
Apply Newton's method to f(x) = x^3-2x+2 starting from x_0 = 0, and explain what happens.
Solution
Step 1 — the iteration. f'(x) = 3x^2-2, so
x_{n+1} = x_n - \frac{x_n^3-2x_n+2}{3x_n^2-2}
Step 2 — iterate from x_0 = 0.
x_1 = 0 - \frac{0-0+2}{0-2} = 0 - \frac{2}{-2} = 0+1 = 1
x_2 = 1 - \frac{1-2+2}{3-2} = 1 - \frac{1}{1} = 0
x_3 = 0 \to 1, \qquad x_4 = 1\to 0, \qquad \ldots
The iteration is trapped in a two-cycle: 0\to1\to0\to1 forever. It never converges and never diverges, and no amount of patience helps.
Step 3 — where is the actual root? f(-2) = -8+4+2 = -2 and f(-1) = -1+2+2 = 3, so there is a root between -2 and -1; it is at x = -1.76929. Newton starting from 0 never goes anywhere near it.
Step 4 — what a fix looks like. Three options, in order of how much they cost you.
Start elsewhere. From x_0 = -2: x_1 = -2-\frac{-2}{10} = -1.8, then -1.76951, then -1.76929. Three steps, converged.
Bracket first. Use a few bisection steps to get into the neighbourhood of the root, then switch to Newton. This is what production solvers do.
Damp the step. Take x_{n+1} = x_n - \lambda\frac{f}{f'} with \lambda<1 chosen so that |f| actually decreases. This is a line search, and it converts most cycles into convergence.
Answer: the iteration cycles between 0 and 1 and never finds the root at -1.769.
The lesson. Newton's quadratic convergence is a local promise, not a global one. It says that once you are close enough, the digits double. It says nothing about getting close, and a plausible-looking starting point can produce a cycle, a divergence, or convergence to the wrong root. Every real solver hedges against this.
Problem 5 — Numerical differentiation, and the optimal step size
Estimate f'(1) for f(x)=e^x using the central difference with h = 0.1, 0.01, 10^{-6} and 10^{-12}. The true value is e = 2.718281828.
Solution
Step 1 — the formula.
f'(x)\approx\frac{f(x+h)-f(x-h)}{2h}
Step 2 — compute each case.
h = 0.1:
\frac{e^{1.1}-e^{0.9}}{0.2} = \frac{3.0041660-2.4596031}{0.2} = \frac{0.5445629}{0.2} = 2.7228147
Error: 4.53\times10^{-3}
h = 0.01:
\frac{e^{1.01}-e^{0.99}}{0.02} = \frac{2.7456011-2.6912345}{0.02} = 2.7183271
Error: 4.53\times10^{-5} — a hundredfold improvement for a tenfold smaller h, exactly as the O(h^2) law predicts.
h = 10^{-6}: the truncation error should now be about 4.5\times10^{-13}, but rounding takes over. The two values of e^x agree to about 6 digits, so their difference has only 10 good digits left, and dividing by 2\times10^{-6} gives roughly
Error: \approx 3\times10^{-11} — still excellent, and about the best achievable.
h = 10^{-12}: now e^{1+h} and e^{1-h} agree in their first 12 digits. Their difference retains only about 4 significant digits, and dividing by 2\times10^{-12} amplifies the rounding noise enormously.
Error: \approx 10^{-4} or worse, and the result varies unpredictably with the exact bits involved.
Step 3 — the summary.
| h | error |
|---|---|
| 0.1 | 4.5\times10^{-3} |
| 0.01 | 4.5\times10^{-5} |
| 10^{-6} | \approx3\times10^{-11} |
| 10^{-12} | \approx10^{-4} |
Answer: accuracy improves down to about h = 10^{-6} and then gets dramatically worse.
Step 4 — where the optimum comes from. Total error is roughly
E(h) \approx \frac{h^2}{6}\left|f'''\right| + \frac{\varepsilon|f|}{h}
Differentiate and set to zero:
\frac{h}{3}|f'''| = \frac{\varepsilon|f|}{h^2} \quad \Rightarrow \quad h^3 = \frac{3\varepsilon|f|}{|f'''|} \quad \Rightarrow \quad h \approx \varepsilon^{1/3}
With \varepsilon = 2.2\times10^{-16} that is h\approx 6\times10^{-6} ✓ matching the table.
Why this matters beyond the exercise. Anyone computing a gradient by finite differences — in optimisation, in sensitivity analysis, in a physics simulation — faces this. The instinct "use the smallest h the computer allows" is exactly wrong, and it produces gradients that look like noise.
Problem 6 — Trapezoid against Simpson
Estimate \int_0^1 e^{-x^2}dx with 4 intervals, by both rules, and compare with the true value 0.7468241.
Solution
Step 1 — tabulate. With h = 0.25:
| i | x_i | f(x_i) = e^{-x_i^2} |
|---|---|---|
| 0 | 0.00 | 1.0000000 |
| 1 | 0.25 | 0.9394131 |
| 2 | 0.50 | 0.7788008 |
| 3 | 0.75 | 0.5697829 |
| 4 | 1.00 | 0.3678794 |
Step 2 — trapezoid. Ends counted once, interiors twice:
T = \frac{0.25}{2}\left[1.0000000 + 2(0.9394131+0.7788008+0.5697829)+0.3678794\right]
The interior sum is 2.2879968, doubled to 4.5759936.
T = 0.125\left[1.0000000+4.5759936+0.3678794\right] = 0.125\times5.9438730 = 0.7429841
Error: 0.7468241-0.7429841 = 3.84\times10^{-3}
Step 3 — Simpson. Pattern 1, 4, 2, 4, 1:
S = \frac{0.25}{3}\left[1.0000000+4(0.9394131)+2(0.7788008)+4(0.5697829)+0.3678794\right]
4(0.9394131) = 3.7576524, \quad 2(0.7788008)=1.5576016, \quad 4(0.5697829) = 2.2791316
Sum: 1.0000000+3.7576524+1.5576016+2.2791316+0.3678794 = 8.9622650
S = 0.0833333\times8.9622650 = 0.7468554
Error: 3.13\times10^{-5}
Answers: trapezoid 0.742984 (error 3.8\times10^{-3}); Simpson 0.746855 (error 3.1\times10^{-5}).
Simpson is 123 times more accurate for the same five function evaluations. That is the h^4 error law against the h^2 one.
Step 4 — confirm the scaling by halving h. With 8 intervals, theory predicts the trapezoid error falls by 4 and Simpson's by 16:
- trapezoid: 3.84\times10^{-3}/4 = 9.6\times10^{-4}
- Simpson: 3.13\times10^{-5}/16 = 2.0\times10^{-6}
Both match the computed values to within a few percent ✓
Step 5 — Richardson extrapolation for free. Combine the two trapezoid estimates:
\frac{4T_{h/2}-T_h}{3} = \frac{4(0.7458656)-0.7429841}{3} = 0.7468261
which is better than Simpson at the same cost, because it is Simpson — the extrapolation of two trapezoid rules is algebraically identical to Simpson's rule, which is a pleasant thing to discover rather than be told.
Why this integral matters. \int e^{-x^2}dx has no elementary antiderivative — this is the Gaussian of 7.5 — the distribution zoo, and every normal-distribution table ever printed was produced by numerical integration exactly like this.
Problem 7 — Euler against Runge–Kutta
Solve y' = y with y(0)=1 from t=0 to t=1, using one step of h=1, then four steps of h=0.25, by Euler and by RK4. The true answer is e = 2.718282.
Solution
Euler, one step of h=1.
y_1 = y_0+hf(t_0,y_0) = 1+1(1) = 2
Error: 0.718. Badly wrong — the method assumed the slope stayed at 1 for the whole unit interval, when in fact it grew to 2.7.
Euler, four steps of h=0.25. Since f = y, each step multiplies by (1+h) = 1.25:
y_4 = (1.25)^4 = 2.441406
Error: 0.277. Better, and the error fell by roughly a factor of 2.6 for a 4-fold reduction in h — consistent with O(h) once the accumulation over more steps is accounted for.
RK4, one step of h=1. Compute the four slopes:
k_1 = f(y_0) = 1
k_2 = f\left(y_0+\tfrac h2k_1\right) = 1+0.5(1) = 1.5
k_3 = f\left(y_0+\tfrac h2k_2\right) = 1+0.5(1.5) = 1.75
k_4 = f\left(y_0+hk_3\right) = 1+1(1.75) = 2.75
y_1 = 1+\frac16\left[1+2(1.5)+2(1.75)+2.75\right] = 1+\frac{1+3+3.5+2.75}{6} = 1+\frac{10.25}{6} = 2.708333
Error: 0.0099 — in a single step, and 28 times more accurate than Euler with four.
RK4, four steps of h=0.25. Each step multiplies by the same factor, which for f=y works out to
1+h+\frac{h^2}{2}+\frac{h^3}{6}+\frac{h^4}{24} = 1.2840254
That expression is the first five terms of the series for e^{h} — RK4 is reproducing the Taylor series of the true solution up to h^4, which is exactly what "fourth-order" means.
y_4 = (1.2840254)^4 = 2.718209
Error: 7.3\times10^{-5}.
Summary.
| Method | steps | result | error |
|---|---|---|---|
| Euler | 1 | 2.000000 | 7.2\times10^{-1} |
| Euler | 4 | 2.441406 | 2.8\times10^{-1} |
| RK4 | 1 | 2.708333 | 9.9\times10^{-3} |
| RK4 | 4 | 2.718209 | 7.3\times10^{-5} |
Answer: RK4 with four steps is about 3,800 times more accurate than Euler with four steps.
The cost comparison that matters. RK4 evaluates f four times per step, so four RK4 steps cost 16 evaluations against Euler's 4. To reach RK4's accuracy, Euler would need roughly 2.8\times10^{-1}/7.3\times10^{-5} = 3800 times more steps — about 15,000 evaluations. Paying four times more per step to need four thousand times fewer steps is why nobody uses Euler for real work.
Problem 8 — Gradient descent by hand
Minimise f(x,y) = x^2+10y^2 starting from (1,1), with learning rate \alpha = 0.1. Do four steps and explain what you see.
Solution
Step 1 — the gradient.
\nabla f = (2x,\ 20y)
Step 2 — the update rule.
x_{n+1} = x_n - 0.1(2x_n) = 0.8x_n
y_{n+1} = y_n - 0.1(20y_n) = y_n - 2y_n = -y_n
Step 3 — iterate.
| n | x | y | f |
|---|---|---|---|
| 0 | 1.000 | 1.000 | 11.000 |
| 1 | 0.800 | -1.000 | 10.640 |
| 2 | 0.640 | 1.000 | 10.410 |
| 3 | 0.512 | -1.000 | 10.262 |
| 4 | 0.410 | 1.000 | 10.168 |
What is happening. The x coordinate is converging nicely, shrinking by 20% each step. The y coordinate is flipping sign and never shrinking at all — it oscillates between +1 and -1 forever. The function value is barely moving, and it is stuck near 10 because the 10y^2 term never improves.
Step 4 — diagnose it. Convergence in a coordinate needs |1-\alpha\lambda|<1, where \lambda is that direction's curvature.
- For x: curvature 2, so 1-0.1(2) = 0.8 ✓ converges
- For y: curvature 20, so 1-0.1(20) = -1.0 ✗ exactly on the boundary, giving a permanent oscillation
Step 5 — the fix, and its cost. Convergence needs
\alpha < \frac{2}{\lambda_{\max}} = \frac{2}{20} = 0.1
so any \alpha below 0.1 works. Take \alpha = 0.09:
x: \times0.82 \text{ per step}, \qquad y: \times(-0.8) \text{ per step}
Both converge, but the x direction now shrinks by only 18% per step — it is being held back by the step size that the y direction requires.
Answer: with \alpha=0.1 the y coordinate oscillates forever; any \alpha<0.1 converges, but slowly.
The general principle, and why it has a name. The problem's condition number is the ratio of the largest to the smallest curvature:
\kappa = \frac{20}{2} = 10
Gradient descent needs roughly \kappa iterations per digit of accuracy. Here \kappa=10 is mild. In a real machine-learning problem \kappa can be 10^4 or worse, which is why the raw method is unusable and why momentum, adaptive step sizes and normalisation exist — every one of them is an attempt to reduce the effective condition number rather than to compute the gradient better.
Problem 9 — Linear programming, solved by inspection of corners
A workshop makes chairs (£40 profit) and tables (£70 profit). Each chair needs 2 hours of carpentry and 1 hour of finishing; each table needs 5 hours of carpentry and 1.5 hours of finishing. There are 80 carpentry hours and 30 finishing hours available. What should they make?
Solution
Step 1 — write the problem. Let c = chairs and t = tables.
\text{maximise } P = 40c+70t
\text{subject to } 2c+5t\le80 \quad\text{(carpentry)}
c+1.5t\le30 \quad\text{(finishing)}
c\ge0,\ t\ge0
Step 2 — find the corners of the feasible region. By the fundamental theorem, the optimum is at one of them, so listing them is enough.
Corner A: (0,0). Make nothing. P = 0.
Corner B: t=0, finishing binding. c = 30. Check carpentry: 2(30) = 60\le80 ✓ feasible. P = 1200.
Corner C: both constraints tight. Solve the pair:
2c+5t = 80
c+1.5t = 30 \quad \Rightarrow \quad c = 30-1.5t
Substitute:
2(30-1.5t)+5t = 80 \quad \Rightarrow \quad 60-3t+5t = 80 \quad \Rightarrow \quad 2t = 20 \quad \Rightarrow \quad t = 10
Then c = 30-15 = 15. P = 40(15)+70(10) = 600+700 = 1300.
Corner D: c=0, carpentry binding. t = 16. Check finishing: 1.5(16) = 24\le30 ✓ feasible. P = 1120.
Step 3 — compare.
| Corner | (c,t) | Profit |
|---|---|---|
| A | (0,0) | £0 |
| B | (30,0) | £1200 |
| C | (15,10) | £1300 |
| D | (0,16) | £1120 |
Answer: make 15 chairs and 10 tables, for £1,300 profit.
Step 4 — check the resources are used properly.
Carpentry: 2(15)+5(10) = 80 — fully used. Finishing: 15+1.5(10) = 30 — fully used.
Both constraints are tight, which is why this corner won: any slack resource is money left on the table.
Step 5 — the shadow prices, which are the useful output. Solve the dual, or equivalently ask what one extra hour of each would be worth. Adding one carpentry hour changes the optimum to c = 14.5, t=10.5, giving P = 580+735 = 1315 — an increase of £15. Adding one finishing hour instead gives c=17.5, t=9, so P = 700+630 = 1330 — an increase of £30.
y_{\text{carpentry}} = £15/\text{hour}, \qquad y_{\text{finishing}} = £30/\text{hour}
This is the answer a manager actually wants. Not "make 15 and 10", which will change next week, but "overtime in finishing is worth up to £30 an hour and overtime in carpentry only £15" — which tells them where to hire, where to invest, and what to refuse.
Problem 10 — Everything at once: a satellite orbit, computed
A satellite's altitude above a target follows \frac{dh}{dt} = -0.5h+3\sin(t), with h(0)=10. (a) Estimate h(2) with RK4 using h_{\text{step}} = 1. (b) Find the time in [0,4] where the altitude is greatest, by finding a root of the derivative. (c) Estimate the average altitude over [0,4] by Simpson's rule.
Solution
Write f(t,h) = -0.5h+3\sin t.
(a) Two RK4 steps of size 1.
Step 1, from t=0, h=10:
k_1 = f(0,10) = -5+0 = -5
k_2 = f(0.5,\ 10+0.5(-5)) = f(0.5, 7.5) = -3.75+3(0.479426) = -3.75+1.438 = -2.312
k_3 = f(0.5,\ 10+0.5(-2.312)) = f(0.5, 8.844) = -4.422+1.438 = -2.984
k_4 = f(1,\ 10+1(-2.984)) = f(1, 7.016) = -3.508+3(0.841471) = -3.508+2.524 = -0.984
h_1 = 10+\frac16\left[-5+2(-2.312)+2(-2.984)+(-0.984)\right] = 10+\frac{-16.576}{6} = 10-2.763 = 7.237
Step 2, from t=1, h=7.237:
k_1 = f(1, 7.237) = -3.619+2.524 = -1.095
k_2 = f(1.5,\ 7.237-0.548) = f(1.5, 6.689) = -3.345+3(0.997495) = -3.345+2.992 = -0.353
k_3 = f(1.5,\ 7.237-0.177) = f(1.5, 7.060) = -3.530+2.992 = -0.538
k_4 = f(2,\ 7.237-0.538) = f(2, 6.699) = -3.350+3(0.909297) = -3.350+2.728 = -0.622
h_2 = 7.237+\frac16\left[-1.095+2(-0.353)+2(-0.538)+(-0.622)\right] = 7.237+\frac{-3.499}{6} = 7.237-0.583 = 6.654
Answer (a): h(2)\approx 6.654.
The exact solution, obtainable by the integrating factor of 6.1 — first order, is h(t) = 11.4e^{-0.5t}+1.2\sin t - 2.4\cos t, giving h(2) = 4.194+1.091+0.999 = 6.284. Our two coarse steps are out by about 6% — respectable for a step size as large as the timescale of the problem, and it would fall to about 0.04% at h_{\text{step}} = 0.25.
(b) The maximum altitude. The altitude peaks where \frac{dh}{dt} = 0:
-0.5h(t)+3\sin t = 0
Using the exact solution, define
g(t) = 3\sin t - 0.5\left(11.4e^{-0.5t}+1.2\sin t-2.4\cos t\right)
Evaluate:
g(2) = 2.728 - 0.5(6.284) = 2.728-3.142 = -0.414
g(1.5) = 2.992 - 0.5(5.385+1.197+(-0.170)) = 2.992 - 3.206 = -0.214
g(1) = 2.524-0.5(6.915+1.010-1.297) = 2.524-3.314 = -0.790
All negative — the altitude is falling throughout, so within [0,2] the maximum is at the start. Check further out:
g(3) = 0.423-0.5(2.544+0.169+2.376) = 0.423-2.545 = -2.121
Still falling. The altitude decreases monotonically on [0,4], so the maximum on the interval is h(0)=10 at t=0.
Answer (b): the greatest altitude on [0,4] is 10 at t=0; the derivative never returns to zero because the decaying term 11.4e^{-0.5t} dominates the driving term over this range.
(c) The average altitude, which is \frac{1}{4}\int_0^4h\,dt. Use Simpson with four intervals, h_{\text{step}} = 1, from the exact solution:
| t | h(t) |
|---|---|
| 0 | 10.000 |
| 1 | 6.628 |
| 2 | 6.284 |
| 3 | 5.089 |
| 4 | 3.673 |
\int_0^4 h\,dt \approx \frac13\left[10.000+4(6.628)+2(6.284)+4(5.089)+3.673\right]
= \frac13\left[10.000+26.512+12.568+20.356+3.673\right] = \frac{73.109}{3} = 24.370
\text{average} = \frac{24.370}{4} = 6.09
Answer (c): the average altitude over the four units of time is about 6.09.
What this problem demonstrates. Three different numerical tools were needed for three different questions about the same system — an integrator for the trajectory, a root-finder for the turning point, and a quadrature rule for the average — and each carried its own error behaviour. Knowing which tool answers which question, and how much to trust each answer, is what numerical analysis is for. The alternative, in a real mission, is a satellite in the wrong place.
Next: Part 11 — The Great Stories, where the mathematics of the previous ten Parts is put back into the hands of the people who found it.