6.12 Numerical Methods: Root Finding, Integration & ODE Solvers

Why Numerical Methods Appear on a PSU Paper

Closed-form solutions run out quickly in engineering practice: the Colebrook equation for pipe friction, the transcendental frequency equation for a continuous beam, and most fin and radiation problems admit no algebraic solution. The syllabus therefore lists three numerical topics, and objective items on them are highly formulaic — recall the iteration formula, do two rounds of arithmetic.

Root Finding: Bisection Method

If $f(a)$ and $f(b)$ have opposite signs, a root lies between them. Repeatedly halve the interval, keeping the half in which the sign change persists:

$x_{\text{mid}} = \frac{a + b}{2}$

Bisection is always convergent but only linearly so: each iteration gains exactly one bit of accuracy, and after $n$ steps the error bound is $(b-a)/2^{n}$. It cannot find roots of even multiplicity, where the function touches but does not cross the axis.

Root Finding: Newton-Raphson Method

The workhorse method, derived from a first-order Taylor expansion:

$x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}$

Convergence is quadratic near a simple root — the number of correct digits roughly doubles each iteration. The failure modes are equally examinable:

• $f'(x_n) \approx 0$ sends the next iterate far away (a near-horizontal tangent).
• A poor initial guess can diverge or cycle between two values.
• At a root of multiplicity greater than one, convergence degrades to linear.

Worked example. Find $\sqrt{10}$ by solving $f(x) = x^2 - 10 = 0$ with $x_0 = 3$.

$x_1 = 3 - \frac{9 - 10}{6} = 3 + 0.16667 = 3.16667$

$x_2 = 3.16667 - \frac{10.0278 - 10}{6.33333} = 3.16667 - 0.00439 = 3.16228$

Two iterations already give six-figure accuracy against the true value 3.162278. Note the special case: applying Newton-Raphson to $x^2 - N$ gives the classical square-root iteration $x_{n+1} = \tfrac{1}{2}\left(x_n + N/x_n\right)$.

Root Finding: Secant and Regula Falsi

The secant method replaces the derivative by a finite difference:

$x_{n+1} = x_n - f(x_n)\frac{x_n - x_{n-1}}{f(x_n) - f(x_{n-1})}$

It needs two starting points and no derivative, converging at order about 1.618. Regula falsi (false position) uses the same formula but always retains a bracketing pair, trading some speed for guaranteed convergence.

Numerical Integration

Both rules divide $[a,b]$ into $n$ strips of width $h = (b-a)/n$ with ordinates $y_0, y_1, \ldots, y_n$.

Trapezoidal rule

$\int_a^b y\,dx \approx \frac{h}{2}\Big[(y_0 + y_n) + 2(y_1 + y_2 + \cdots + y_{n-1})\Big]$

Any number of intervals is allowed. The rule fits straight lines between points, so it is exact for linear integrands and its error is proportional to $h^2$.

Simpson's one-third rule

$\int_a^b y\,dx \approx \frac{h}{3}\Big[(y_0 + y_n) + 4(y_1 + y_3 + \cdots) + 2(y_2 + y_4 + \cdots)\Big]$

The weighting pattern is $1, 4, 2, 4, 2, \ldots, 4, 1$ — odd-indexed ordinates get 4, even-indexed interior ordinates get 2. It requires an even number of intervals (an odd number of ordinates) because it fits a parabola through each group of three points. It is exact up to cubic integrands and its error goes as $h^4$.

Simpson's three-eighths rule

$\int_a^b y\,dx \approx \frac{3h}{8}\Big[(y_0 + y_n) + 3(y_1 + y_2 + y_4 + y_5 + \cdots) + 2(y_3 + y_6 + \cdots)\Big]$

This requires the number of intervals to be a multiple of 3.

Worked example. Evaluate $\int_0^2 x^3 dx$ with $n = 2$, so $h = 1$ and ordinates $y_0 = 0$, $y_1 = 1$, $y_2 = 8$.

Simpson: $\dfrac{1}{3}\left[(0 + 8) + 4(1)\right] = \dfrac{12}{3} = 4$ — the exact answer, since Simpson is exact for cubics.

Trapezoidal: $\dfrac{1}{2}\left[(0+8) + 2(1)\right] = 5$, an error of 25%. This contrast is a favourite examination point.

Numerical Solution of Differential Equations

For $\dfrac{dy}{dx} = f(x,y)$ with $y(x_0) = y_0$:

Euler's method (single-step, first order):

$y_{n+1} = y_n + h\,f(x_n, y_n)$

Simple but crude; the error per step is proportional to $h^2$ and accumulates.

Modified Euler / Heun (predictor-corrector, second order):

$y_{n+1} = y_n + \frac{h}{2}\left[f(x_n,y_n) + f(x_{n+1}, y^{*}_{n+1})\right]$

Fourth-order Runge-Kutta, the standard multi-stage method:

$k_1 = hf(x_n,y_n), \quad k_2 = hf\!\left(x_n + \tfrac{h}{2}, y_n + \tfrac{k_1}{2}\right)$

$k_3 = hf\!\left(x_n + \tfrac{h}{2}, y_n + \tfrac{k_2}{2}\right), \quad k_4 = hf(x_n + h, y_n + k_3)$

$y_{n+1} = y_n + \frac{1}{6}\left(k_1 + 2k_2 + 2k_3 + k_4\right)$

Four function evaluations per step buy fourth-order accuracy, which is why RK4 is the default in almost every engineering solver. Euler's method is the special case obtained by keeping only $k_1$.