Computational Methods

Most equations you meet outside the textbook cannot be solved by a formula, and most integrals cannot be evaluated by anti-differentiation. Computational methods give up on exact answers and instead build a sequence of numbers that closes in on the true value, stopping when the answer is accurate enough for the job. This unit is the first place in your course where "accurate enough" is made precise. Note that the CDC syllabus offers Mechanics (Unit 21) as the alternative to this unit for the same 12 hours; you study one of the two, and Mechanics is covered in the next chapter.

20.1 Roots of algebraic and transcendental equations: bisection method

An algebraic equation is one built from powers of x, such as x^{3}-x-1=0. A transcendental equation involves trigonometric, exponential or logarithmic functions, such as \cos x = x or xe^{x} = 1. A number \alpha with f(\alpha)=0 is a root of f(x)=0. Beyond the quadratic (and the rarely used cubic) formula there is no algebraic recipe for the roots, so we locate them numerically.

Locating a root

Every method starts by trapping the root inside an interval.

To find such an interval, tabulate f at convenient integers and look for a change of sign. For f(x) = x^{3}-x-1:

x	0	1	2
f(x)	-1	-1	+5

The sign changes between x=1 and x=2, so a root lies in [1,2].

The method

Because the root always stays inside the current interval, the midpoint after n steps differs from the true root by at most half the current interval width:

This is the great virtue of bisection — it cannot fail, and it tells you in advance how many steps you need. To guarantee an error below \varepsilon, solve

20.2 Roots by Newton-Raphson method

Bisection is safe but slow: each step buys only one bit of accuracy. Newton-Raphson uses the gradient as well as the value of f, and converges dramatically faster.

Geometrically, x_{n+1} is the point where the tangent to y=f(x) at (x_{n}, f(x_{n})) meets the x-axis. Replacing the curve by its tangent is a good approximation near the root, so each new point lands much closer than the last.

Newton-Raphson for roots and reciprocals

Setting f(x) = x^{2}-N gives x_{n+1} = \frac{1}{2}\left(x_{n}+\frac{N}{x_{n}}\right), the classical square-root algorithm; f(x) = x^{3}-N gives x_{n+1} = \frac{1}{3}\left(2x_{n}+\frac{N}{x_{n}^{2}}\right). These are how calculators actually compute roots.

How fast, and when it fails

Newton-Raphson converges quadratically: near a simple root the error is roughly squared at each step, so the number of correct decimal places doubles every iteration. That is why three or four rows are normally enough, and why you should carry eight decimals in the working.

The price of the speed is that Newton-Raphson can fail. It breaks down if f'(x_{n}) = 0 (the tangent is horizontal and never meets the axis), it can be thrown far away if the starting value is poor, and it can fall into a cycle.

	Bisection	Newton-Raphson
Needs	an interval with a sign change	one starting value and f'(x)
Always converges?	yes	no
Speed	slow (linear)	fast (quadratic)
Error known in advance?	yes, \frac{b-a}{2^{n}}	no
Stop when	interval < tolerance	two successive values agree

20.3 Numerical integration: trapezoidal rule

When an integrand has no elementary anti-derivative — \sqrt{1+x^{3}}, e^{-x^{2}}, \frac{\sin x}{x} — or when f is known only as a table of measured values, we approximate the area directly.

Divide [a,b] into n equal strips of width h = \frac{b-a}{n} with ordinates

The trapezoidal rule replaces the curve over each strip by the chord joining its two endpoints, so each strip becomes a trapezium of area \frac{h}{2}(y_{r}+y_{r+1}). Adding all n of them, every interior ordinate is counted twice.

20.4 Numerical integration: Simpson's rule

Instead of joining the endpoints of each strip by a straight line, Simpson's rule fits a parabola through three consecutive points. A parabola bends, so it tracks a curve far better than a chord does.

Because three points are needed for one parabola, the strips are consumed in pairs, so n must be even. Over the pair [x_{0},x_{2}] the area under the fitted parabola works out to \frac{h}{3}(y_{0}+4y_{1}+y_{2}); adding all such pairs gives the alternating pattern of coefficients 4,2,4,2,\ldots,4.

	Trapezoidal rule	Simpson's one-third rule
Curve replaced by	straight chords	parabolic arcs
Restriction on n	none	must be even
Coefficients	1,2,2,\ldots,2,1	1,4,2,4,\ldots,4,1
Multiplier	\frac{h}{2}	\frac{h}{3}
Error	\propto h^{2}	\propto h^{4}
Exact for	linear f	cubic f and below

20.5 Extension: linear programming by the graphical method

Optimising a linear quantity subject to linear constraints is the other classical computational method of school mathematics, and it uses the same "compute, tabulate, compare" discipline. A linear programming problem asks you to maximise or minimise an objective function Z = ax+by subject to linear inequalities called constraints, together with the non-negativity conditions x \ge 0, y \ge 0.

Each inequality shades a half-plane; their intersection is the feasible region. Because Z is linear, its level curves ax+by = k form a family of parallel lines, and sliding that line across the feasible region shows that the optimum is always reached at a corner point (vertex) of the region. So the method is: draw the region, find its vertices, evaluate Z at each, and pick the best.

Chapter summary

• A root of f(x)=0 is bracketed whenever f(a) \,f(b) < 0 and f is continuous on [a,b].
• Bisection: repeatedly take c = \frac{a+b}{2} and keep the half in which the sign change lies. It always converges, and |x_{n}-\alpha| \le \frac{b-a}{2^{n}}, so n > \frac{\ln((b-a)/\varepsilon)}{\ln 2} iterations guarantee accuracy \varepsilon.
• Newton-Raphson: x_{n+1} = x_{n} - \frac{f(x_{n})}{f'(x_{n})}, obtained by replacing the curve with its tangent. Convergence is quadratic (correct digits double each step) but is not guaranteed; it fails when f'(x_{n}) is zero or the start is poor.
• Special cases: x_{n+1} = \frac12\left(x_{n}+\frac{N}{x_{n}}\right) for \sqrt{N}, and x_{n+1} = \frac{1}{3}\left(2x_{n}+\frac{N}{x_{n}^{2}}\right) for \sqrt[3]{N}.
• Trapezoidal rule with h = \frac{b-a}{n}: \int_{a}^{b}f\,dx \approx \frac{h}{2}\left[(y_{0}+y_{n})+2(y_{1}+\cdots+y_{n-1})\right], with |E_{T}| \le \frac{(b-a)h^{2}}{12}M_{2}. Any n is allowed.
• Simpson's one-third rule (n even): \int_{a}^{b}f\,dx \approx \frac{h}{3}\left[(y_{0}+y_{n})+4(\text{odd }y)+2(\text{even }y)\right], with |E_{S}| \le \frac{(b-a)h^{4}}{180}M_{4}. Simpson's three-eighths rule uses \frac{3h}{8} with coefficients 1,3,3,2,3,3,\ldots,1 and needs n a multiple of three.
• Simpson's rule is exact for any cubic; the trapezoidal rule is exact only for a straight line and always errs to the same side of a curve of fixed concavity.
• In a linear programming problem the optimum of Z = ax+by occurs at a corner of the feasible region, so evaluate Z at every vertex and compare.
• Round only at the end, keep six decimals throughout, and always state the accuracy claimed and the error estimate when the question asks for one.

Practice questions

Group A — Multiple choice (1 mark each)

1. The bisection method can be started on [a,b] only if [1]af(a) = fb (b) f(a)\,f(b) < 0cf(a)\,f(b) > 0df'(a)=0
2. After n bisections of [a,b] the error in the midpoint is at most [1]a\frac{b-a}{n}b\frac{b-a}{2n}c\frac{b-a}{2^{n}}d\frac{b-a}{n^{2}}
3. The Newton-Raphson iteration formula is [1]ax_{n+1}=x_{n}+\frac{f(x_{n})}{f'(x_{n})}bx_{n+1}=x_{n}-\frac{f(x_{n})}{f'(x_{n})}cx_{n+1}=x_{n}-\frac{f'(x_{n})}{f(x_{n})}dx_{n+1}=\frac{x_{n}}{f'(x_{n})}
4. Newton-Raphson fails at x_{n} when [1]af(x_{n})=0bf'(x_{n})=0cx_{n}>0df''(x_{n})=0
5. In Simpson's one-third rule the number of sub-intervals must be [1]aoddbevenca multiple of 3da prime
6. In the trapezoidal rule the multiplier outside the bracket is [1]ahb\frac{h}{2}c\frac{h}{3}d\frac{3h}{8}
7. Simpson's one-third rule integrates exactly every polynomial of degree at most [1]a1b2c3d4
8. Which equation is transcendental? [1]ax^{5}-3x+1=0b2x^{2}-7=0cxe^{x}-1=0dx^{3}=8
9. Using f(x)=x^{2}-5 with x_{0}=2, the first Newton-Raphson iterate is [1]a2.2500b2.2361c2.5000d2.0000
10. The error of the trapezoidal rule is proportional to [1]ahbh^{2}ch^{3}dh^{4}
11. In a linear programming problem the optimal value of the objective function occurs [1]aat the origin alwaysbat a corner point of the feasible regioncat the centre of the regiondoutside the region

Answers to Group A

1. (b) — opposite signs guarantee a crossing by the location theorem.

2. (c) — the interval halves every step, so its width is \frac{b-a}{2^{n}}.

3. (b) — subtract \frac{f}{f'}; this is the tangent meeting the x-axis.

4. (b) — a zero derivative means a horizontal tangent, which never meets the axis, and the formula divides by zero.

5. (b) — strips are used in pairs, one parabola per pair.

6. (b) — each trapezium has area \frac{h}{2}(y_r+y_{r+1}).

7. (c) — the cubic terms cancel in the error analysis, so Simpson's rule is exact for cubics.

8. (c) — it contains e^{x}; the others are polynomial equations.

9. (a) — x_1 = 2 - \frac{4-5}{4} = 2.25.

10. (b) — |E_T| \le \frac{(b-a)h^{2}}{12}M_2.

11. (b) — the level lines of a linear Z are parallel, so the extreme is pushed to a vertex.

Group B — Short answer (5 marks each)

1. Show that the equation x^{3}-x-1=0 has a root between 1 and 2, and find it correct to one decimal place by the bisection method. [5]
2. Find a root of xe^{x} = 1 correct to two decimal places using the bisection method on [0,1]. [5]
3. Derive the Newton-Raphson iteration formula from Taylor's theorem, and state two circumstances in which it fails. [5]
4. Use the Newton-Raphson method with x_{0}=2 to find the root of x^{3}-2x-5=0 correct to five decimal places. [5]
5. Evaluate \int_{0}^{1}e^{x}dx by the trapezoidal rule with four sub-intervals and compare with the exact value. [5]
6. Evaluate \int_{1}^{2}\ln x\,dx by Simpson's rule with four sub-intervals, giving your answer to five decimal places. [5]
7. Show that Simpson's rule gives the exact value of \int_{0}^{2}x^{3}dx with n=4, and explain why. [5]
8. Use x_{n+1} = \frac{1}{5}\left(4x_{n}+\frac{32}{x_{n}^{4}}\right), obtained from Newton-Raphson applied to x^{5}-32=0, with x_{0}=2.1 to compute \sqrt[5]{32} to five decimal places. [5]

Answers to Group B

1. f(1) = 1-1-1 = -1 < 0 and f(2) = 8-2-1 = 5 > 0, so by the location theorem a root lies in [1,2]. Bisecting:

n	a	b	c	f(c)
1	1.0000	2.0000	1.5000	+0.8750
2	1.0000	1.5000	1.2500	-0.2969
3	1.2500	1.5000	1.3750	+0.2246
4	1.2500	1.3750	1.3125	-0.0515
5	1.3125	1.3750	1.3438	+0.0826
6	1.3125	1.3438	1.3281	+0.0146
7	1.3125	1.3281	1.3203	-0.0187

The error bound is now \frac{1}{2^{7}} = 0.0078 < 0.05, so \alpha \approx 1.3 to one decimal place.

2. Let f(x) = xe^{x}-1. Then f(0) = -1 < 0 and f(1) = e-1 = 1.7183 > 0, so a root lies in [0,1].

n	a	b	c	f(c) = ce^{c}-1
1	0.0000	1.0000	0.5000	-0.1756
2	0.5000	1.0000	0.7500	+0.5878
3	0.5000	0.7500	0.6250	+0.1677
4	0.5000	0.6250	0.5625	-0.0128
5	0.5625	0.6250	0.5938	+0.0751
6	0.5625	0.5938	0.5781	+0.0306
7	0.5625	0.5781	0.5703	+0.0088
8	0.5625	0.5703	0.5664	-0.0020

The bound is \frac{1}{2^{8}} = 0.0039 < 0.005, so \alpha \approx 0.57 to two decimal places (true value 0.567143).

3. Let x_{n} approximate the root and let h satisfy f(x_{n}+h) = 0. By Taylor's theorem, f(x_{n}+h) = f(x_{n}) + hf'(x_{n}) + \ldots; neglecting the terms in h^{2} and higher because h is small,

f(x_{n}) + hf'(x_{n}) \approx 0 \;\Longrightarrow\; h \approx -\frac{f(x_{n})}{f'(x_{n})}

so x_{n+1} = x_{n} - \frac{f(x_{n})}{f'(x_{n})}. It fails (i) when f'(x_{n}) = 0, since the tangent is then parallel to the x-axis and the formula divides by zero, and (ii) when the starting value is far from the root, in which case the iterates may diverge or cycle (for example x^{3}-2x+2=0 with x_{0}=0 oscillates between 0 and 1).

4. f(x) = x^{3}-2x-5, f'(x) = 3x^{2}-2, so x_{n+1} = \frac{2x_{n}^{3}+5}{3x_{n}^{2}-2}.

n	x_{n}	f(x_{n})	f'(x_{n})	x_{n+1}
0	2.00000000	-1.00000000	10.00000000	2.10000000
1	2.10000000	+0.06100000	11.23000000	2.09456812
2	2.09456812	+0.00018572	11.16164684	2.09455148

Hence \alpha \approx 2.09455.

5. h = 0.25; ordinates y_r = e^{x_r}:

r	x_{r}	y_{r}
0	0.00	1.000000
1	0.25	1.284025
2	0.50	1.648721
3	0.75	2.117000
4	1.00	2.718282

Ends = 3.718282; interiors = 5.049746. So

T = \frac{0.25}{2}\left[3.718282 + 2(5.049746)\right] = 0.125 \times 13.817774 = 1.727222

The exact value is \left[e^{x}\right]_{0}^{1} = e-1 = 1.718282, so the trapezoidal rule overestimates by 0.008940 — as expected, since y=e^{x} is concave up.

6. h = 0.25; ordinates y_r = \ln x_r:

r	x_{r}	y_{r}	coefficient
0	1.00	0.000000	1
1	1.25	0.223144	4
2	1.50	0.405465	2
3	1.75	0.559616	4
4	2.00	0.693147	1

Ends = 0.693147; odd = 0.782760; even = 0.405465.

S = \frac{0.25}{3}\left[0.693147 + 4(0.782760) + 2(0.405465)\right] = 0.083333 \times 4.635117 = 0.386260

So \int_{1}^{2}\ln x\,dx \approx 0.38626 (exact value 2\ln 2 - 1 = 0.386294).

7. With a=0, b=2, n=4 we get h = 0.5 and ordinates y = 0,\ 0.125,\ 1,\ 3.375,\ 8. Ends = 8; odd = 0.125+3.375 = 3.5; even = 1.

S = \frac{0.5}{3}\left[8 + 4(3.5) + 2(1)\right] = \frac{0.5}{3}\times 24 = 4.000000

The exact value is \left[\frac{x^{4}}{4}\right]_{0}^{2} = 4. They agree because Simpson's error term involves the fourth derivative, and f^{(4)}(x) = 0 for a cubic, so the rule is exact for every polynomial of degree three or less.

8. With x_{0}=2.1:

n	x_{n}	x_{n+1} = \frac{1}{5}\left(4x_{n}+\frac{32}{x_{n}^{4}}\right)
0	2.10000000	2.00908099
1	2.00908099	2.00008172
2	2.00008172	2.00000001

So \sqrt[5]{32} \approx 2.00000, which is right, since 2^{5} = 32.

Group C — Long answer (8 marks each)

1. (a) State the location theorem and use it to show that \cos x = x has a root between 0 and 1 radian. [2](b) Find that root correct to two decimal places by the bisection method, showing the error bound at each stage. [6]
2. (a) Derive the trapezoidal rule for n sub-intervals. [3](b) Hence evaluate \int_{0}^{1}\sqrt{1+x^{3}}\,dx with four sub-intervals, and also apply Simpson's rule to the same ordinates. State which answer you trust more and why. [5]
3. (a) Obtain the Newton-Raphson formula for \sqrt{N} and use it to find \sqrt{13} correct to six decimal places starting from x_{0}=3.5. [5](b) Compare the bisection and Newton-Raphson methods under the headings reliability, speed and information required. [3]
4. A cottage industry in Bhaktapur makes two products. Maximise Z = 4x+3y subject to x+y \le 8, 2x+y \le 10, x \ge 0, y \ge 0, by the graphical method, listing all corner points and the value of Z at each. [8]

Answers to Group C

1. (a) Location theorem: if f is continuous on [a,b] and f(a)\,f(b) < 0, then f(x) = 0 has at least one root in (a,b). Put f(x) = \cos x - x, which is continuous. Then f(0) = 1 - 0 = +1 and f(1) = 0.540302 - 1 = -0.459698. The signs are opposite, so a root lies in (0,1).

(b)

n	a	b	c	f(c)	error bound \frac{1}{2^{n}}
1	0.000000	1.000000	0.500000	+0.377583	0.500000
2	0.500000	1.000000	0.750000	-0.018311	0.250000
3	0.500000	0.750000	0.625000	+0.185963	0.125000
4	0.625000	0.750000	0.687500	+0.085335	0.062500
5	0.687500	0.750000	0.718750	+0.033879	0.031250
6	0.718750	0.750000	0.734375	+0.007875	0.015625
7	0.734375	0.750000	0.742188	-0.005196	0.007813
8	0.734375	0.742188	0.738281	+0.001345	0.003906

The bound 0.003906 is below 0.005, so the root is 0.74 correct to two decimal places.

2. (a) Divide [a,b] into n strips of width h = \frac{b-a}{n} with ordinates y_{0},y_{1},\ldots,y_{n}. Replacing the curve over the r-th strip by the chord makes it a trapezium of parallel sides y_{r} and y_{r+1} and width h, so its area is \frac{h}{2}(y_{r}+y_{r+1}). Summing from r=0 to n-1, each interior ordinate appears in two trapezia and each end ordinate in one:

\int_{a}^{b}f(x)dx \approx \frac{h}{2}\left[(y_{0}+y_{n}) + 2(y_{1}+y_{2}+\cdots+y_{n-1})\right]

(b) h = 0.25 and y_r = \sqrt{1+x_r^{3}}:

r	x_{r}	y_{r}
0	0.00	1.000000
1	0.25	1.007782
2	0.50	1.060660
3	0.75	1.192424
4	1.00	1.414214

Trapezoidal: ends = 2.414214, interiors = 3.260866, so

T = \frac{0.25}{2}\left[2.414214 + 2(3.260866)\right] = 0.125 \times 8.935946 = 1.116993

Simpson: odd = 1.007782+1.192424 = 2.200206, even = 1.060660, so

S = \frac{0.25}{3}\left[2.414214 + 4(2.200206) + 2(1.060660)\right]

= 0.083333 \times 13.336158 = 1.111363

Simpson's value is the more reliable, because its error is of order h^{4} against h^{2} for the trapezoidal rule. (The true value is 1.111448, so Simpson is out by 0.000085 and the trapezoidal rule by 0.005545.)

3. (a) Put f(x) = x^{2}-N, so f'(x) = 2x. Then

x_{n+1} = x_{n} - \frac{x_{n}^{2}-N}{2x_{n}} = \frac{2x_{n}^{2}-x_{n}^{2}+N}{2x_{n}} = \frac{1}{2}\left(x_{n}+\frac{N}{x_{n}}\right)

With N=13 and x_{0}=3.5:

n	x_{n}	x_{n+1}
0	3.50000000	3.60714286
1	3.60714286	3.60555163
2	3.60555163	3.60555128

Hence \sqrt{13} \approx 3.605551.

(b) Reliability: bisection always converges provided the initial interval brackets a root; Newton-Raphson may diverge or cycle. Speed: bisection gains about one binary digit per iteration (linear convergence), whereas Newton-Raphson roughly doubles the number of correct digits each iteration (quadratic convergence), so it typically needs three or four steps against ten or more. Information required: bisection needs only values of f and a sign change; Newton-Raphson needs the derivative f' and a good starting value, but no bracketing interval.

4. Draw x+y=8 through (8,0) and (0,8), and 2x+y=10 through (5,0) and (0,10). Substituting the origin into each constraint gives 0 \le 8 and 0 \le 10, both true, so the feasible region lies on the origin side of both lines, in the first quadrant. Solving the two boundary equations simultaneously: subtracting x+y=8 from 2x+y=10 gives x=2, hence y=6. The corner points are therefore O(0,0), C(5,0), B(2,6) and A(0,8).

Corner point	Z = 4x+3y
O(0,0)	0
C(5,0)	20 + 0 = 20
B(2,6)	8 + 18 = 26
A(0,8)	0 + 24 = 24

The feasible region is a closed bounded polygon, so the maximum is attained at a corner. The greatest value is Z = 26, at x = 2 and y = 6.