1.5 Numerical Methods, Root-Finding, and Numerical Integration

1.5 Numerical Methods, Root-Finding, and Numerical Integration

Engineering systems frequently involve non-linear equations, complicated experimental datasets, or differential equations lacking analytical closed-form solutions. Numerical methods provide approximate algorithms to solve such mathematical models within specified tolerances.

---

1. Numerical Root-Finding Algorithms

Root-finding algorithms locate values of $x$ for which $f(x) = 0$.

Bisection Method (Bracketing Method)

Requires an initial interval $[a_0, b_0]$ where $f(a_0)$ and $f(b_0)$ have opposite signs ($f(a_0) \cdot f(b_0) < 0$). By the Intermediate Value Theorem, a root exists in $[a_0, b_0]$.

1. Midpoint calculation: $c_k = \frac{a_k + b_k}{2}$.
2. Evaluate $f(c_k)$:
   • If $f(a_k) \cdot f(c_k) < 0$, the root lies in $[a_k, c_k] \implies b_{k+1} = c_k, a_{k+1} = a_k$.
   • If $f(c_k) \cdot f(b_k) < 0$, the root lies in $[c_k, b_k] \implies a_{k+1} = c_k, b_{k+1} = b_k$.
3. Maximum absolute error bound after $n$ iterations: $E_n \le \frac{b_0 - a_0}{2^n}$

Newton-Raphson Method (Open Method)

Uses the tangent line at current guess $x_n$ to find the next approximation $x_{n+1}$: $x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}$

Convergence & Limitations

• Convergence Rate: Quadratic convergence ($E_{n+1} \propto E_n^2$) near simple roots.
• Failure Modes: Fails if $f'(x_n) = 0$ (horizontal tangent), oscillates near local extrema, or diverges if initial guess $x_0$ is too far from the root.

Secant Method

Replaces derivative $f'(x_n)$ in Newton-Raphson with a backward finite difference approximation using two points $x_n$ and $x_{n-1}$: $x_{n+1} = x_n - f(x_n) \frac{x_n - x_{n-1}}{f(x_n) - f(x_{n-1})}$

---

2. Numerical Integration (Quadrature)

Numerical integration estimates definite integrals $I = \int_a^b f(x) dx$ over $n$ subintervals of equal width $h = \frac{b - a}{n}$ with grid points $x_i = a + i h$.

Trapezoidal Rule

Approximates $f(x)$ as linear segments between grid points: $I_{trap} = \frac{h}{2} \left[ f(x_0) + 2 f(x_1) + 2 f(x_2) + \cdots + 2 f(x_{n-1}) + f(x_n) \right]$ $I_{trap} = \frac{h}{2} \left[ f(x_0) + 2 \sum_{i=1}^{n-1} f(x_i) + f(x_n) \right]$

• Global Truncation Error: $E_T = -\frac{(b - a) h^2}{12} f''(\xi)$, where $a \le \xi \le b$. Error is $O(h^2)$.

Simpson's 1/3 Rule

Approximates $f(x)$ using parabolic arcs across pairs of subintervals. Strict Requirement: The number of subintervals $n$ must be an even integer.

$I_{simp} = \frac{h}{3} \left[ f(x_0) + 4 f(x_1) + 2 f(x_2) + 4 f(x_3) + 2 f(x_4) + \cdots + 4 f(x_{n-1}) + f(x_n) \right]$ $I_{simp} = \frac{h}{3} \left[ f(x_0) + 4 \sum_{i \text{ odd}} f(x_i) + 2 \sum_{i \text{ even}} f(x_i) + f(x_n) \right]$

• Global Truncation Error: $E_S = -\frac{(b - a) h^4}{180} f^{(4)}(\xi)$. Error is $O(h^4)$.
• Exact for polynomials up to degree 3!

---

3. Numerical Solution of Initial Value ODEs

For first-order initial value problem $\frac{dy}{dx} = f(x,y)$ with $y(x_0) = y_0$ and step size $h$:

Euler's Method (Explicit First-Order)

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

• Local Error: $O(h^2)$, Global Error: $O(h)$ (First-order accurate).

Improved Euler Method (Heun's Method - Predictor-Corrector)

1. Predictor: $y_{n+1}^* = y_n + h \cdot f(x_n, y_n)$
2. Corrector: $y_{n+1} = y_n + \frac{h}{2} \left[ f(x_n, y_n) + f(x_{n+1}, y_{n+1}^*) \right]$

• Global Error: $O(h^2)$ (Second-order accurate).

---

4. Linear Interpolation Formula

For a tabular dataset where $y_1 = f(x_1)$ and $y_2 = f(x_2)$, the interpolated value $y$ at $x$ ($x_1 \le x \le x_2$) is: $y = y_1 + \left( \frac{y_2 - y_1}{x_2 - x_1} \right) (x - x_1)$

---

5. Worked Engineering Problems

Worked Example 1: Newton-Raphson Iteration

Problem: Solve $f(x) = x^3 - 2x - 5 = 0$ using two iterations of Newton-Raphson starting with $x_0 = 2.0$.

Solution:

Step 1: Compute derivative $f'(x)$. $f'(x) = 3x^2 - 2$

Step 2: First Iteration ($n=0$). $f(x_0) = f(2.0) = 2^3 - 2(2) - 5 = 8 - 4 - 5 = -1.0$ $f'(x_0) = f'(2.0) = 3(2.0)^2 - 2 = 12 - 2 = 10.0$ $x_1 = x_0 - \frac{f(x_0)}{f'(x_0)} = 2.0 - \frac{-1.0}{10.0} = 2.0 + 0.1 = 2.10$

Step 3: Second Iteration ($n=1$). $f(x_1) = f(2.10) = (2.10)^3 - 2(2.10) - 5 = 9.261 - 4.20 - 5 = 0.061$ $f'(x_1) = f'(2.10) = 3(2.10)^2 - 2 = 3(4.41) - 2 = 13.23 - 2 = 11.23$ $x_2 = 2.10 - \frac{0.061}{11.23} = 2.10 - 0.00543 = 2.09457$

---

Worked Example 2: Numerical Integration Comparison

Problem: Estimate $I = \int_0^4 (x^3 + 2) dx$ using (a) Trapezoidal Rule with $n=4$ subintervals ($h=1$) and (b) Simpson's 1/3 Rule with $n=4$ ($h=1$). Compare both against the exact analytical solution.

Solution:

**Grid points & Function evaluations ($h = 1$):

• $x_0 = 0 \implies f(0) = 0^3 + 2 = 2$
• $x_1 = 1 \implies f(1) = 1^3 + 2 = 3$
• $x_2 = 2 \implies f(2) = 2^3 + 2 = 10$
• $x_3 = 3 \implies f(3) = 3^3 + 2 = 29$
• $x_4 = 4 \implies f(4) = 4^3 + 2 = 66$

(a) Trapezoidal Rule: $I_{trap} = \frac{1}{2} \left[ f(0) + 2(f(1) + f(2) + f(3)) + f(4) \right]$ $I_{trap} = \frac{1}{2} \left[ 2 + 2(3 + 10 + 29) + 66 \right] = \frac{1}{2} \left[ 2 + 2(42) + 66 \right] = \frac{1}{2} [2 + 84 + 66] = \frac{152}{2} = 76.0$

(b) Simpson's 1/3 Rule: $I_{simp} = \frac{1}{3} \left[ f(0) + 4(f(1) + f(3)) + 2(f(2)) + f(4) \right]$ $I_{simp} = \frac{1}{3} \left[ 2 + 4(3 + 29) + 2(10) + 66 \right] = \frac{1}{3} \left[ 2 + 4(32) + 20 + 66 \right] = \frac{1}{3} [2 + 128 + 20 + 66] = \frac{216}{3} = 72.0$

(c) Exact Analytical Integral: $I_{exact} = \left[ \frac{x^4}{4} + 2x \right]_0^4 = \left(\frac{256}{4} + 8\right) - 0 = 64 + 8 = 72.0$

Notice: Simpson's 1/3 Rule yields the exact value (72.0) because its truncation error is zero for cubic polynomials!