Numerical Methods - Computational Algorithms

Table of Contents

1. Root Finding
2. Numerical Integration
3. Numerical Differentiation
4. Interpolation
5. Numerical Linear Algebra
6. Initial Value Problems for ODEs
7. Finite Difference Methods
8. Finite Element Method

---

Root Finding

Bisection Method

Assumption: $f$ continuous on $[a,b]$ with $f(a)f(b) < 0$

Algorithm:

1. Set $c = (a+b)/2$
2. If $f(c) = 0$ (or $|f(c)| < \epsilon$): done
3. If $f(a)f(c) < 0$: set $b = c$, else set $a = c$
4. Repeat

Convergence: Guaranteed, linear rate (error reduces by factor 1/2 each step)

Error bound after $n$ steps: $|c_n - r| \leq (b-a)/2^{n+1}$

Fixed Point Iteration

Find fixed point of g: $x = g(x)$

Algorithm: $x_{n+1} = g(x_n)$

Convergence: If $|g'(x)| < 1$ near fixed point $r$

Rate: If $g'(r) = 0$ (higher order), convergence is quadratic.

Newton’s Method

For $f(x) = 0$:

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

Convergence: Quadratic if $f'(r) \neq 0$

Order: $p = 2$

Geometric: Tangent line intersection with x-axis

Issue: May diverge or oscillate for poor initial guess

Secant Method

Approximates derivative:

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

Convergence: Superlinear, order $\approx 1.618$

Advantage: No derivative needed

Disadvantage: Two initial guesses needed

Convergence Criteria

Absolute error: $|x_n - x_{n-1}| < \epsilon$

Relative error: $|x_n - x_{n-1}|/|x_n| < \epsilon$

Function value: $|f(x_n)| < \epsilon$

Multi-dimensional Newton

For $f(x) = 0$ (system):

$x_{n+1} = x_n - J^{-1}(x_n) f(x_n)$

where $J$ is Jacobian matrix.

Requires: Solving linear system at each iteration

---

Numerical Integration

Newton-Cotes Formulas

Closed formulas: Use endpoints

Trapezoidal Rule

$\int_a^b f(x)dx \approx \frac{b-a}{2}[f(a) + f(b)]$

Composite:

$\int_a^b f(x)dx \approx h\left[\frac{f(a)}{2} + \sum_{i=1}^{n-1} f(x_i) + \frac{f(b)}{2}\right]$

where $h = (b-a)/n, x_i = a + ih$

Error: $O(h^2)$

Simpson’s Rule

$\int_a^b f(x)dx \approx \frac{b-a}{6}[f(a) + 4f(\frac{a+b}{2}) + f(b)]$

Composite ($n$ even):

$\int_a^b f(x)dx \approx \frac{h}{3}[f(x_0) + 4\sum_{i=1}^{n/2} f(x_{2i-1}) + 2\sum_{i=1}^{n/2-1} f(x_{2i}) + f(x_n)]$

Error: $O(h^4)$

Requires: Even number of subintervals

Adaptive Integration

Adapt subinterval size based on local error estimates.

Romberg Integration: Richardson extrapolation on trapezoidal rule

Gauss Quadrature: Best $n$-node formula is orthogonal polynomial method

Gauss-Legendre: Nodes are roots of Legendre polynomials

---

Numerical Differentiation

Forward Difference

$f'(x) \approx \frac{f(x+h) - f(x)}{h}$

Error: $O(h)$

Backward Difference

$f'(x) \approx \frac{f(x) - f(x-h)}{h}$

Error: $O(h)$

Central Difference

$f'(x) \approx \frac{f(x+h) - f(x-h)}{2h}$

Error: $O(h^2)$

Second Derivative

$f''(x) \approx \frac{f(x+h) - 2f(x) + f(x-h)}{h^2}$

Error: $O(h^2)$

Sensitivity to Roundoff

Numerical derivatives amplify rounding errors:

Relative error: $\epsilon/h$

Optimal $h$: Balance truncation and roundoff errors.

---

Interpolation

Polynomial Interpolation

Given $(x_i, y_i)$: Find polynomial $p(x)$ with $p(x_i) = y_i$

Existence: Unique polynomial of degree $\leq n$ for $n+1$ points

Lagrange Interpolation

$p(x) = \sum_{i=0}^n y_i L_i(x)$

where $L_i(x) = \prod_{j \neq i} \frac{x-x_j}{x_i-x_j}$

Properties:

• $L_i(x_i) = 1$
• $L_i(x_j) = 0$ for $j \neq i$
• Basis for polynomials

Newton Divided Differences

Divided differences:

$f[x_0] = f(x_0)$ $f[x_0,x_1] = \frac{f(x_1) - f(x_0)}{x_1 - x_0}$ $f[x_0,...,x_k] = \frac{f[x_1,...,x_k] - f[x_0,...,x_{k-1}]}{x_k - x_0}$

Newton form:

$p_n(x) = f[x_0] + \sum_{k=1}^n f[x_0,...,x_k] \prod_{j=0}^{k-1}(x-x_j)$

Advantage: Easy to add points (compute recursively)

Hermite Interpolation

Interpolate function and derivatives

Given $(x_i, y_i, y'_i)$, find $p$ with $p(x_i) = y_i$ and $p'(x_i) = y'_i$

Splines

Cubic spline: Piecewise cubic with continuous 1st and 2nd derivatives

Natural spline: $S'' = 0$ at endpoints

Clamped spline: $S'$ prescribed at endpoints

Setup: System of equations for cubic coefficients

---

Numerical Linear Algebra

Gaussian Elimination with Pivoting

Partial pivoting: Swap row with largest $|a_{ij}|$ in column

Complete pivoting: Swap both rows and columns

Gauss-Jordan: Continue elimination to reduce form

LU Decomposition

$PA = LU$ where $P$ permutation, $L$ lower, $U$ upper triangular

Solve: $Ly = Pb$, then $Ux = y$

Advantage: $L, U$ computed once, solves multiple systems

Cholesky Decomposition

For positive definite A: $A = LL^T$ ($L$ lower triangular)

More efficient than LU

QR Decomposition

$A = QR$ where $Q$ orthogonal, $R$ upper triangular

Gram-Schmidt: Compute column by column

Householder reflections: More numerically stable

$Q^T Q = I$

Eigenvalue Algorithms

Power method: For dominant eigenvalue

$x_{k+1} = \frac{A x_k}{||A x_k||}$

Converges to eigenvector, eigenvalue from Rayleigh quotient

Inverse iteration: For eigenvalue near $\mu$

Solve $(A - \mu I)x_{k+1} = x_k$

QR algorithm: Iterative method for all eigenvalues

Shifts accelerate convergence

Condition Number

$\kappa(A) = ||A|| ||A^{-1}||$

Indicates sensitivity to errors in linear systems

Large $\kappa \implies$ nearly singular

---

Initial Value Problems for ODEs

Problem

$y' = f(t,y), \quad y(t_0) = y_0$

Euler’s Method

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

where $h$ = step size

Error: $O(h)$

Local error: $O(h^2)$

Global error: $O(h)$ over entire interval

Runge-Kutta Methods

RK4 (Classical):

$k_1 = h f(t_n, y_n)$ $k_2 = h f(t_n + h/2, y_n + k_1/2)$ $k_3 = h f(t_n + h/2, y_n + k_2/2)$ $k_4 = h f(t_n + h, y_n + k_3)$ $y_{n+1} = y_n + (k_1 + 2k_2 + 2k_3 + k_4)/6$

Error: $O(h^4)$

Multistep Methods

Use previous values:

$y_{n+1} = \sum_{j=0}^{k-1} \alpha_j y_{n-j} + h \sum_{j=-1}^{k-1} \beta_j f_{n-j}$

Adams-Bashforth: Explicit ($\beta_{-1} = 0$)

Adams-Moulton: Implicit ($\beta_{-1} \neq 0$, requires solve)

Advantage: More efficient per step

Disadvantage: Need starting values (e.g., RK4)

Stability

Stability region: {$z \in \mathbb{C}$ : method stable}

A-stable: Includes left half-plane

Implicit methods often more stable (larger stability regions)

Adaptive Step Size

Estimate local error: Compare two methods

Adjust $h$: Larger $h$ if error small, smaller if large

---

Finite Difference Methods

Approximating Derivatives

Second derivative:

$u''(x_i) \approx \frac{u(x_{i+1}) - 2u(x_i) + u(x_{i-1})}{h^2}$

Boundary Value Problems

Example: $u'' = f(x), u(0) = \alpha, u(1) = \beta$

Discretize: $x_i = ih, i = 0,...,n+1$ ($h = 1/(n+1)$)

Discrete equations:

$\frac{u_{i+1} - 2u_i + u_{i-1}}{h^2} = f(x_i)$

Linear system: $Au = b$

Tridiagonal matrix (banded)

Two-Point BVP

Eigenvalue problem: $-u'' = \lambda u, u(0) = u(1) = 0$

Eigenvalues: $\lambda_k = (k\pi)^2, k = 1, 2, \dots$

Use QZ algorithm for discrete problem

2D Problems

Laplace equation: $u_{xx} + u_{yy} = 0$

Discretize: $u_{i,j} \approx u(x_i, y_j)$

Five-point stencil:

$\frac{u_{i+1,j} - 2u_{i,j} + u_{i-1,j}}{h_x^2} + \frac{u_{i,j+1} - 2u_{i,j} + u_{i,j-1}}{h_y^2} = 0$

Block tridiagonal system

---

Finite Element Method

Weak Formulation

Instead of solving strong form, solve weak form:

Find $u$ such that $\int_{\Omega} \nabla u \cdot \nabla v \, d\Omega = \int_{\Omega} fv \, d\Omega$

for all test functions $v$

Galerkin Method

Approximate $u$ as: $u \approx \sum u_i \phi_i$

Choose test functions $\psi_j$

Discretize: $\sum u_i \int \nabla\phi_i \cdot \nabla\psi_j = \int f \psi_j$

Galerkin: Choose $\psi_j = \phi_j$

Linear system: $Ku = F$

where $K_{ij} = \int \nabla\phi_i \cdot \nabla\phi_j$

Basis Functions

Piecewise linear: $\phi_i(x)$ hat function at node $i$

Piecewise polynomial: Higher order elements

Advantages

Adaptive meshing Handles complex geometries Irregular meshes

Mesh Generation

Triangular elements (2D) Tetrahedral elements (3D)

---

Supplementary Numerical Methods Reference (from mathematics_GPT)

Error Estimates for Numerical Differentiation

Forward Difference:

$f'(x) \approx \frac{f(x+h) - f(x)}{h}$

Error: $O(h)$

Backward Difference:

$f'(x) \approx \frac{f(x) - f(x-h)}{h}$

Error: $O(h)$

Central Difference:

$f'(x) \approx \frac{f(x+h) - f(x-h)}{2h}$

Error: $O(h^2)$

Error Estimates for Numerical Integration

Trapezoidal Rule:

$\int_a^b f(x)dx \approx \frac{b-a}{2}[f(a) + f(b)]$

Composite:

$\int_a^b f(x)dx \approx h\left[\frac{f(a)}{2} + \sum_{i=1}^{n-1} f(x_i) + \frac{f(b)}{2}\right]$

Error: $O(h^2)$

Simpson’s Rule:

$\int_a^b f(x)dx \approx \frac{b-a}{6}[f(a) + 4f(\frac{a+b}{2}) + f(b)]$

Composite ($n$ even):

$\int_a^b f(x)dx \approx \frac{h}{3}[f(x_0) + 4\sum_{i=1}^{n/2} f(x_{2i-1}) + 2\sum_{i=1}^{n/2-1} f(x_{2i}) + f(x_n)]$

Error: $O(h^4)$

Error Estimates for ODE Solvers

Euler’s Method:

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

Error: $O(h)$

Local error: $O(h^2)$

Global error: $O(h)$ over entire interval

Runge-Kutta Methods:

RK4 (Classical):

$k_1 = h f(t_n, y_n)$ $k_2 = h f(t_n + h/2, y_n + k_1/2)$ $k_3 = h f(t_n + h/2, y_n + k_2/2)$ $k_4 = h f(t_n + h, y_n + k_3)$ $y_{n+1} = y_n + (k_1 + 2k_2 + 2k_3 + k_4)/6$

Error: $O(h^4)$

Error Estimates for Linear Systems

Gaussian Elimination with Pivoting:

Partial pivoting: Swap row with largest $|a_{ij}|$ in column

Complete pivoting: Swap both rows and columns

Error: $O(n^3)$

LU Decomposition:

$PA = LU$ where $P$ permutation, $L$ lower, $U$ upper triangular

Error: $O(n^3)$

Cholesky Decomposition:

For positive definite A: $A = LL^T$ ($L$ lower triangular)

Error: $O(n^3/6)$

Error Estimates for Eigenvalue Problems

Power method: For dominant eigenvalue

$x_{k+1} = \frac{A x_k}{||A x_k||}$

Error: $O(\epsilon$)

Inverse iteration: For eigenvalue near $\mu$

Error: $O(\epsilon$)

QR algorithm: Iterative method for all eigenvalues

Error: $O(n^3)$

Error Estimates for Condition Number

$\kappa(A) = ||A|| ||A^{-1}||$

Indicates sensitivity to errors in linear systems

Large $\kappa \implies$ nearly singular

Error: $O(n^2)$

---

Next: Optimization or Topology

---

Last updated: Comprehensive numerical methods reference covering root finding, integration, and differential equation solvers.