Numerical Analysis: Roots, Quadrature, and ODEs

Numerical analysis is the study of algorithms that obtain approximate solutions to mathematical problems where exact analytical solutions are either impossible to find or computationally expensive. This lesson covers the rigorous foundations of root-finding, integration, and the numerical solution of differential equations.

1. Root-Finding Algorithms

Root-finding involves determining a value $x^{\ast }$ such that $f(x^{\ast })=0$.

The Bisection Method

Based on the Intermediate Value Theorem (IVT), if $f(x)$ is continuous on $[a,b]$ and $f(a)f(b)<0$, then there exists at least one $c\in (a,b)$ such that $f(c)=0$.

The algorithm iteratively halves the interval $[a,b]$. The absolute error after $n$ iterations is bounded by: $|e_n|=|x^{\ast }-x_n|\le \frac{b-a}{2^n}$ This method is robust but converges linearly.

Newton-Raphson Method

Derived from the linear Taylor expansion of $f(x)$ around $x_n$: $f(x)\approx f(x_n)+f^{\prime }(x_n)(x-x_n)$ Setting $f(x)=0$ and solving for $x_{n+1}$ yields: $x_{n+1}=x_n-\frac{f(x_n)}{f^{\prime }(x_n)}$

Convergence Analysis: If $f\in C^2[a,b]$ and $x^{\ast }$ is a simple root ($f^{\prime }(x^{\ast })\ne 0$), Newton’s method exhibits quadratic convergence: ${\lim }_{n\to \infty }\frac{|e_{n+1}|}{|e_n{|}^2}=\frac{|f^{\prime \prime }(x^{\ast })|}{2|f^{\prime }(x^{\ast })|}$ However, it requires a “sufficiently close” initial guess $x_0$.

Fixed-Point Iteration (FPI)

An equation $f(x)=0$ can often be rewritten as $x=g(x)$. The iteration $x_{n+1}=g(x_n)$ converges to a fixed point $x^{\ast }$ if $g$ satisfies the conditions of the Banach Fixed-Point Theorem.

Contraction Mapping Principle: If $g$ is a contraction mapping on a closed interval $I$ (i.e., $|g^{\prime }(x)|\le L<1$ for all $x\in I$), then:

1. $x^{\ast }$ is unique in $I$.
2. The sequence $x_{n+1}=g(x_n)$ converges to $x^{\ast }$ for any $x_0\in I$.
3. Convergence is linear with rate $L$.

2. Order of Convergence

We formally define the order of convergence $p$ and the asymptotic error constant $C$: ${\lim }_{n\to \infty }\frac{|e_{n+1}|}{|e_n{|}^p}=C$

• $p=1$ (Linear): Bisection, FPI.
• $1<p<2$ (Superlinear): Secant Method ($p\approx 1.618$).
• $p=2$ (Quadratic): Newton’s Method.

3. Numerical Integration (Quadrature)

Numerical integration approximates $I={\int }_a^{b}f(x)dx$ via a weighted sum ${\sum }_{i=0}^n{w}_{i}f(x_i)$.

Newton-Cotes Formulas

These use equally spaced nodes.

• Trapezoidal Rule: Approximates $f$ with a first-degree polynomial. ${\int }_a^{b}f(x)dx\approx \frac{h}{2}[f(a)+f(b)],E=-\frac{h^3}{12}{f}^{\prime \prime }(\xi )$
• Simpson’s 1/3 Rule: Fits a parabola through three points ($h=\frac{b-a}{2}$). ${\int }_a^{b}f(x)dx\approx \frac{h}{3}[f(x_0)+4f(x_1)+f(x_2)],E=-\frac{h^5}{90}{f}^{(4)}(\xi )$

Gaussian Quadrature

Unlike Newton-Cotes, Gaussian quadrature chooses both weights $w_i$ and nodes $x_i$ to maximize the degree of polynomials integrated exactly. For an $n$-point rule, nodes are the roots of the $n$-th degree Legendre Polynomial $P_n(x)$ on $[-1,1]$. An $n$-point Gaussian quadrature rule is exact for all polynomials of degree $2n-1$ or less.

4. Numerical Ordinary Differential Equations (ODEs)

Consider the Initial Value Problem (IVP): $\frac{dy}{dt}=f(t,y),y(t_0)=y_0$

Euler’s Method

The simplest approach using the first term of the Taylor expansion: $y_{n+1}=y_n+hf(t_n,y_n)$

• Local Truncation Error (LTE): $O(h^2)$ (error in one step).
• Global Truncation Error (GTE): $O(h)$ (accumulated error).

Runge-Kutta Methods (RK4)

The “Classic” RK4 method achieves $O(h^4)$ accuracy by taking a weighted average of four incremental slopes: $k_1=f(t_n,y_n)$ $k_2=f(t_n+\frac{h}{2},y_n+\frac{h}{2}{k}_1)$ $k_3=f(t_n+\frac{h}{2},y_n+\frac{h}{2}{k}_2)$ $k_4=f(t_n+h,y_n+h{k}_3)$ $y_{n+1}=y_n+\frac{h}{6}(k_1+2k_2+2k_3+k_4)$

5. Stability and Stiffness

Stiffness

A differential equation is stiff if certain terms cause rapid variations in the solution, requiring extremely small step sizes for explicit methods like Euler or RK4 to remain stable.

Implicit Methods

To solve stiff equations, we use Implicit Euler (Backward Euler): $y_{n+1}=y_n+hf(t_{n+1},y_{n+1})$ This requires solving an algebraic equation at each step (often via Newton’s method) but allows for much larger step sizes because it is A-stable.

Python Implementation: RK4 vs. Euler

Below we compare the accuracy of Euler’s method against the 4th-order Runge-Kutta method for the ODE $y^{\prime }=-2t{y}^2$ with $y(0)=1$. The exact solution is $y(t)=\frac{1}{1+t^2}$.

import numpy as np

def f(t, y):
    return -2 * t * y**2

def exact_sol(t):
    return 1 / (1 + t**2)

def euler_step(t, y, h):
    return y + h * f(t, y)

def rk4_step(t, y, h):
    k1 = f(t, y)
    k2 = f(t + h/2, y + h/2 * k1)
    k3 = f(t + h/2, y + h/2 * k2)
    k4 = f(t + h, y + h * k3)
    return y + (h/6) * (k1 + 2*k2 + 2*k3 + k4)

# Simulation parameters
t0, y0 = 0, 1
t_end = 2
h = 0.2
steps = int((t_end - t0) / h)

t_vals = np.linspace(t0, t_end, steps + 1)
y_euler = np.zeros(steps + 1)
y_rk4 = np.zeros(steps + 1)
y_euler[0] = y0
y_rk4[0] = y0

for i in range(steps):
    y_euler[i+1] = euler_step(t_vals[i], y_euler[i], h)
    y_rk4[i+1] = rk4_step(t_vals[i], y_rk4[i], h)

# Error analysis at t = 2
exact = exact_sol(t_end)
print(f"Exact solution: {exact:.6f}")
print(f"Euler result:   {y_euler[-1]:.6f} (Error: {abs(exact - y_euler[-1]):.6e})")
print(f"RK4 result:     {y_rk4[-1]:.6f} (Error: {abs(exact - y_rk4[-1]):.6e})")

Knowledge Check