Numerical Methods

Numerical methods are the board's answer to problems that have no clean closed form: a root that will not factor, an integral with no elementary antiderivative, a table of measured values you must read between. The CELE tests a small, dependable kit, root finding by bisection and Newton-Raphson, linear interpolation between tabulated points, the trapezoidal and Simpson rules for area under a curve, and the finite-difference formulas that turn derivatives into arithmetic. Every one of them is a short, mechanical recipe, and this lesson carries each recipe all the way to a finished number so you can check your own habits against it.

Root finding by bisection

Bisection is the slow but sure method. It needs a starting interval [a, b] on which the function changes sign, that is, f(a) and f(b) have opposite signs, which by the intermediate value theorem guarantees at least one root between them. Each step evaluates the midpoint c = (a + b) / 2, checks the sign of f(c), and keeps whichever half still brackets the sign change. The interval is cut in half every step, so after n steps the bracket width is (b - a) / 2^n, and the midpoint estimate is within (b - a) / 2^(n+1) of the true root.

Worked example: Find sqrt(2) as the positive root of f(x) = x^2 - 2 on [1, 2]. Check the ends: f(1) = 1 - 2 = -1 and f(2) = 4 - 2 = 2, opposite signs, so a root is bracketed. Step 1: c = 1.5, f(1.5) = 2.25 - 2 = 0.25, positive, so the root is in [1, 1.5]. Step 2: c = 1.25, f(1.25) = 1.5625 - 2 = -0.4375, negative, so the root is in [1.25, 1.5]. Step 3: c = 1.375, f(1.375) = 1.890625 - 2 = -0.109375, negative, so the root is in [1.375, 1.5]. Step 4: c = 1.4375, f(1.4375) = 2.066 - 2 = 0.066, positive, so the root is in [1.375, 1.4375]. After 4 steps the bracket width is 1 / 2^4 = 0.0625 and the best estimate, the midpoint 1.40625, is within 0.03125 of the true value sqrt(2) = 1.41421. Slow, but it never fails once you have a valid bracket.

Newton-Raphson method

Newton-Raphson trades the guarantee of bisection for speed. Starting from a single guess x_n, it replaces the curve with its tangent line at that point and takes the next estimate to be where that tangent crosses the x-axis:

x_(n+1) = x_n - f(x_n) / f'(x_n)

Near a simple root the error roughly squares each step (quadratic convergence), so the number of correct digits doubles per iteration. The price is that you must supply the derivative f'(x), and the method can misbehave: if f'(x_n) is near zero the tangent is nearly flat and the step flies off, and a poor starting guess can send the sequence to a different root or diverge entirely.

Worked example: Solve x^2 - 2 = 0 again, this time by Newton-Raphson, so you can compare. Here f(x) = x^2 - 2 and f'(x) = 2x, so the iteration is x_(n+1) = x_n - (x_n^2 - 2) / (2 x_n), which simplifies to the classic average x_(n+1) = (x_n + 2 / x_n) / 2. Start at x_0 = 1.5. Then x_1 = (1.5 + 2 / 1.5) / 2 = (1.5 + 1.333333) / 2 = 1.416667. Next x_2 = (1.416667 + 2 / 1.416667) / 2 = (1.416667 + 1.411765) / 2 = 1.414216. One more step reproduces 1.414214, which matches sqrt(2) to six digits. Newton reached in two steps what bisection was still crawling toward after four.

Step	Bisection estimate	Newton estimate
start	bracket [1, 2]	x_0 = 1.5
1	1.5	1.416667
2	1.25	1.414216
3	1.375	1.414214
4	1.4375	1.414214

The contrast is the whole point: bisection adds about one bit of accuracy per step, while Newton doubles the correct digits per step once it is close.

Linear interpolation

When you only have a table of values, linear interpolation estimates a value between two known points by drawing a straight line between them. Given points (x0, y0) and (x1, y1), the value at an intermediate x is:

y = y0 + (y1 - y0) times (x - x0) / (x1 - x0)

The bracketed factor (x - x0) / (x1 - x0) is just the fraction of the way from x0 to x1, and you apply that same fraction to the change in y. Reading a value the other way, finding the x that gives a target y, is called inverse interpolation and is exactly how the false-position (regula falsi) root-finding method locates a root: interpolate the x where the straight line between (a, f(a)) and (b, f(b)) crosses zero.

Worked example (reading a table): A rating table lists a weir discharge coefficient of C = 0.611 at a head of 0.20 m and C = 0.593 at a head of 0.40 m. Estimate C at a head of 0.30 m. Here x0 = 0.20, y0 = 0.611, x1 = 0.40, y1 = 0.593, and x = 0.30. The fraction of the way across is (0.30 - 0.20) / (0.40 - 0.20) = 0.10 / 0.20 = 0.5, so C = 0.611 + (0.593 - 0.611)(0.5) = 0.611 + (-0.018)(0.5) = 0.611 - 0.009 = 0.602. The estimate is C = 0.602.

Worked example (inverse interpolation for a root): Take f(x) = x^2 - 2 with the bracket a = 1, b = 2, so f(a) = -1 and f(b) = 2. The straight line between (1, -1) and (2, 2) crosses zero at x = a - f(a) times (b - a) / (f(b) - f(a)) = 1 - (-1)(2 - 1) / (2 - (-1)) = 1 + 1 / 3 = 1.3333. That single false-position estimate, 1.3333, already beats the first bisection midpoint 1.5 as an approximation to sqrt(2) = 1.41421.

Numerical integration: the trapezoidal rule

To integrate a function you cannot integrate by hand, or one known only at sample points, slice the interval [a, b] into n equal strips of width h = (b - a) / n and approximate the area of each strip. The trapezoidal rule caps each strip with a straight chord, making it a trapezoid, and sums them:

integral from a to b of f(x) dx = (h / 2)(y0 + 2 y1 + 2 y2 + ... + 2 y_(n-1) + yn)

The two end ordinates carry a coefficient of 1 and every interior ordinate carries a coefficient of 2. Forgetting to double the interior ordinates is the most common slip.

Worked example: Estimate the integral of x^2 from 0 to 2 with n = 4, so h = 0.5. The nodes are x = 0, 0.5, 1.0, 1.5, 2.0 with ordinates y = 0, 0.25, 1.0, 2.25, 4.0. Apply the rule: T = (0.5 / 2)[0 + 2(0.25 + 1.0 + 2.25) + 4.0] = 0.25[0 + 2(3.5) + 4.0] = 0.25(11.0) = 2.75. The exact value is 2^3 / 3 = 8 / 3 = 2.6667, so the trapezoidal estimate 2.75 overshoots by 0.0833, exactly the concave-up behavior the figure shows.

Simpson's one-third rule

Simpson's one-third rule fits a parabola through each pair of strips instead of a straight chord, so it captures curvature and is far more accurate for the same effort. It requires an even number of strips n (an odd number of ordinates). The pattern of coefficients is 1 at the ends, then 4, 2, 4, 2, ..., 4, 1 through the interior, odd-indexed ordinates weighted 4 and even-indexed ordinates weighted 2:

integral from a to b of f(x) dx = (h / 3)(y0 + 4 y1 + 2 y2 + 4 y3 + ... + 4 y_(n-1) + yn)

A useful fact for the exam: Simpson's rule integrates any cubic or lower polynomial exactly, so it returns the exact answer for parabolas and cubics no matter how coarse the spacing.

Worked example: Redo the integral of x^2 from 0 to 2 with n = 4, h = 0.5, and the same ordinates 0, 0.25, 1.0, 2.25, 4.0. Simpson gives S = (0.5 / 3)[0 + 4(0.25) + 2(1.0) + 4(2.25) + 4.0] = (0.16667)[0 + 1.0 + 2.0 + 9.0 + 4.0] = (0.16667)(16.0) = 2.6667. That equals the exact value 8 / 3 to the last digit, because x^2 is a polynomial of degree below four.

Worked example (a curve neither rule nails): Estimate the integral of 1 / x from 1 to 2 with n = 4, h = 0.25. The ordinates are 1 / 1 = 1, 1 / 1.25 = 0.8, 1 / 1.5 = 0.666667, 1 / 1.75 = 0.571429, and 1 / 2 = 0.5. Trapezoidal: T = (0.25 / 2)[1 + 2(0.8 + 0.666667 + 0.571429) + 0.5] = 0.125(5.576192) = 0.697024. Simpson: S = (0.25 / 3)[1 + 4(0.8) + 2(0.666667) + 4(0.571429) + 0.5] = (0.083333)(8.319048) = 0.693254. The exact value is the natural logarithm ln 2 = 0.693147, so Simpson is off by only 0.000107 while the trapezoidal rule is off by 0.003877, about thirty-six times worse for the same five ordinates.

Method	Interior coefficients	Estimate of integral 1/x from 1 to 2, n = 4	Error against ln 2 = 0.693147
Trapezoidal	1, 2, 2, 2, 1 (times h/2)	0.697024	0.003877
Simpson's 1/3	1, 4, 2, 4, 1 (times h/3)	0.693254	0.000107

Finite differences

Finite differences turn derivatives, defined by a limit, into arithmetic on nearby function values by simply using a small step h instead of the limit. Three first-derivative formulas matter. The forward and backward differences use one neighbor each and carry an error proportional to h. The central difference straddles the point with both neighbors and carries an error proportional to h^2, so halving h cuts its error by about a factor of four rather than two. There is also a central formula for the second derivative.

Approximation	Formula	Error order
Forward difference (f')	[f(x + h) - f(x)] / h	h
Backward difference (f')	[f(x) - f(x - h)] / h	h
Central difference (f')	[f(x + h) - f(x - h)] / (2h)	h^2
Central second difference (f'')	[f(x + h) - 2 f(x) + f(x - h)] / h^2	h^2

Worked example: Estimate the derivatives of f(x) = x^2 at x = 2 using h = 0.1, so f(1.9) = 3.61, f(2.0) = 4.00, f(2.1) = 4.41. The forward difference gives (4.41 - 4.00) / 0.1 = 0.41 / 0.1 = 4.1, and the backward difference gives (4.00 - 3.61) / 0.1 = 0.39 / 0.1 = 3.9; each misses the true value f'(2) = 4 by 0.1. The central difference gives (4.41 - 3.61) / (2 times 0.1) = 0.80 / 0.2 = 4.0, exact here because the two one-sided errors cancel. The central second difference gives (4.41 - 2 times 4.00 + 3.61) / (0.1)^2 = (4.41 - 8.00 + 3.61) / 0.01 = 0.02 / 0.01 = 2.0, matching f''(2) = 2.

The same forward-difference idea drives Euler's method for a first-order differential equation dy/dx = g(x, y): replacing the derivative by [y_(n+1) - y_n] / h and rearranging gives the marching formula y_(n+1) = y_n + h times g(x_n, y_n). Worked example: For dy/dx = x + y with y(0) = 1 and step h = 0.1, the first Euler step is y_1 = 1 + 0.1 times (0 + 1) = 1.1 at x = 0.1, and the second is y_2 = 1.1 + 0.1 times (0.1 + 1.1) = 1.1 + 0.1(1.2) = 1.22 at x = 0.2.

Root method	Needs	Convergence	Watch out for
Bisection	a bracket [a, b] with a sign change	slow but guaranteed (linear)	requires f(a) and f(b) of opposite sign
Newton-Raphson	one start x0 and the derivative f'(x)	very fast (quadratic) near the root	fails if f'(x) is near 0; can diverge from a poor start
False position	a bracket [a, b] with a sign change	faster than bisection	one endpoint can stall for strongly curved functions

Exam-day strategy

• Before starting bisection or false position, confirm the bracket is valid by checking that f(a) and f(b) have opposite signs; without a sign change these methods have nothing to converge to.
• On a Newton-Raphson item, write the iteration x_(n+1) = x_n - f(x_n) / f'(x_n) first, then compute f(x_n) and f'(x_n) separately before combining; a near-zero derivative is the warning sign that the method may blow up.
• For linear interpolation, form the fraction (x - x0) / (x1 - x0) first and confirm it is between 0 and 1; a value outside that range means you are extrapolating, not interpolating, and the estimate is far less reliable.
• Match the integration rule to the ordinate count: the trapezoidal rule works for any number of strips, but Simpson's one-third rule needs an even number of strips (an odd number of ordinates); if the count is odd, you cannot apply Simpson to the whole interval.
• Lock in the coefficient patterns: trapezoidal is 1, 2, 2, ..., 2, 1 times h/2, and Simpson is 1, 4, 2, 4, ..., 4, 1 times h/3; the single most common mistake is dropping the factor of 2 on the interior trapezoidal ordinates or swapping the 4 and 2 weights in Simpson.
• Prefer the central difference over a one-sided difference when the point has neighbors on both sides; its error shrinks like h^2, so it is markedly more accurate for the same step size.