Fundamental Numerical Methods

Why This Matters

Most real-world equations don't have neat, closed-form solutions. Numerical methods transform those impossible analytical problems into solvable computational ones. The core skill here is choosing the right method for a given problem, understanding why certain approaches converge faster, and recognizing the trade-offs between accuracy and computational cost.

These methods underpin everything from climate simulations to machine learning optimization. The concepts you need to master include convergence rates, stability conditions, error propagation, and computational complexity. Don't just memorize formulas. Know when each method shines, when it fails, and what mathematical principles make it work.

---

Finding Where Functions Equal Zero

Root-finding algorithms locate values where $f(x) = 0$. This task appears everywhere, from solving nonlinear equations to optimization subproblems. The key distinction is between bracketing methods (guaranteed but slow) and open methods (fast but potentially unstable).

Bisection Method

• Guaranteed convergence if $f(a)$ and $f(b)$ have opposite signs (by the Intermediate Value Theorem). The method repeatedly halves the interval containing the root.
• Linear convergence rate means each iteration adds roughly one binary digit of accuracy, requiring about 3.3 iterations per decimal digit.
• Robust but slow. Use it when reliability matters more than speed, or as a fallback when faster methods fail.

Newton-Raphson Method

• Quadratic convergence near a simple root doubles the number of correct digits each iteration, using the update formula $x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}$
• Requires derivative computation and a sufficiently close initial guess. It fails when $f'(x) \approx 0$ (division by near-zero) or when the initial guess lands near an inflection point.
• Best for smooth, well-behaved functions where you can compute derivatives analytically or via automatic differentiation.

Secant Method

Worth mentioning as a middle ground: the secant method approximates the derivative using two previous iterates, so you don't need $f'(x)$. Its convergence order is approximately 1.618 (the golden ratio), placing it between bisection and Newton-Raphson in speed.

---

Estimating Values Between Known Points

Interpolation constructs functions that pass exactly through given data points. The challenge is balancing polynomial degree against numerical stability and computational efficiency.

Lagrange Interpolation

• Explicit polynomial formula constructs the unique degree-$n$ polynomial through $n+1$ points using basis polynomials $L_i(x) = \prod_{j \neq i} \frac{x - x_j}{x_i - x_j}$
• Conceptually simple but computationally expensive. Adding one new point requires recalculating all basis polynomials from scratch.
• Prone to Runge's phenomenon: high-degree polynomials on equally spaced nodes can oscillate wildly near interval endpoints. Using Chebyshev nodes mitigates this.

Newton's Divided Differences

• Incremental construction allows adding new data points without complete recalculation, using the form $P(x) = \sum_{i=0}^{n} [y_0, \ldots, y_i] \prod_{j=0}^{i-1}(x - x_j)$
• The divided difference table stores intermediate calculations, making the method more efficient for growing datasets.
• Mathematically equivalent to Lagrange (both produce the same unique polynomial), but superior for practical computation and error estimation.

---

Approximating Definite Integrals

Numerical integration (quadrature) approximates $\int_a^b f(x)\,dx$ when antiderivatives are unavailable or impractical. Higher-order methods achieve better accuracy by using more sophisticated function approximations within each subinterval.

Trapezoidal Rule

• Approximates the integrand as linear segments on each subinterval. The composite formula is $\int_a^b f(x)\,dx \approx \frac{h}{2}[f(a) + 2f(x_1) + \cdots + 2f(x_{n-1}) + f(b)]$
• The error is $O(h^2)$, meaning halving the step size cuts the error by roughly a factor of 4.
• Works surprisingly well for periodic functions integrated over complete periods, where endpoint errors cancel out due to symmetry.

Simpson's Rule

• Uses parabolic approximations on pairs of subintervals, achieving error of $O(h^4)$, a dramatic improvement for smooth functions.
• Requires an even number of subintervals. The composite formula is $\int_a^b f(x)\,dx \approx \frac{h}{3}[f(a) + 4f(x_1) + 2f(x_2) + 4f(x_3) + \cdots + f(b)]$
• Exact for polynomials up to degree 3. This is the sweet spot between simplicity and accuracy for most applications.

---

Computing Derivatives from Discrete Data

Numerical differentiation estimates derivatives when only discrete function values are available. Unlike integration, differentiation amplifies noise: small data errors become large derivative errors.

Finite Difference Approximations

Three standard formulas for $f'(x)$:

• Forward difference: $f'(x) \approx \frac{f(x+h) - f(x)}{h}$ with $O(h)$ error
• Backward difference: $f'(x) \approx \frac{f(x) - f(x-h)}{h}$ with $O(h)$ error
• Central difference: $f'(x) \approx \frac{f(x+h) - f(x-h)}{2h}$ with $O(h^2)$ error

The central difference is more accurate because the leading error terms cancel by symmetry.

The step size dilemma is critical to understand. Smaller $h$ reduces truncation error but amplifies round-off error from floating-point arithmetic. There's an optimal $h$ that balances these two effects (roughly $h \sim \sqrt{\epsilon_{\text{mach}}}$ for first derivatives with central differences, where $\epsilon_{\text{mach}}$ is machine epsilon).

Higher-order stencils use more points for better accuracy. For example, a five-point stencil achieves fourth-order accuracy for first derivatives.

---

Solving Linear Systems

Linear algebra operations form the backbone of scientific computing. Direct methods give exact answers (up to round-off) in a fixed number of operations, while iterative methods approximate solutions with controllable accuracy.

Gaussian Elimination

1. Form the augmented matrix $[A | b]$.
2. Use row operations to reduce to upper triangular (row-echelon) form. This is the forward elimination phase with $O(n^3)$ cost.
3. Apply partial pivoting at each step: swap rows to place the largest-magnitude element on the diagonal. This prevents division by small numbers and improves numerical stability.
4. Solve via back substitution in $O(n^2)$ operations.

LU Decomposition

• Factors the matrix as $A = LU$ (or $PA = LU$ with pivoting), where $L$ is lower triangular and $U$ is upper triangular.
• Efficient for multiple right-hand sides. Once you've computed the factorization ($O(n^3)$), each new system $Ax = b$ requires only $O(n^2)$ forward and back substitution.
• Useful for determinants: $\det(A) = \pm \prod_i U_{ii}$ (the sign depends on the number of row swaps from pivoting).

---

Analyzing Matrix Properties

Eigenvalue problems reveal intrinsic matrix properties essential for understanding linear transformations, stability, and system dynamics. The eigenvalue equation $Ax = \lambda x$ identifies special directions (eigenvectors) that are only scaled, not rotated, under the transformation $A$.

Eigenvalue Computation

• Power iteration finds the dominant eigenvalue by repeatedly multiplying: $x_{k+1} = \frac{Ax_k}{\|Ax_k\|}$. It converges when $|\lambda_1| > |\lambda_2|$, and the convergence rate depends on the ratio $|\lambda_2 / \lambda_1|$.
• QR algorithm iteratively factors $A_k = Q_k R_k$ and forms $A_{k+1} = R_k Q_k$. The matrices $A_k$ converge to a form that reveals all eigenvalues. This is the standard method used in practice (e.g., by MATLAB's eig).
• Applications span disciplines: stability analysis (eigenvalues inside the unit circle for discrete systems, negative real parts for continuous), PCA (largest eigenvalues capture the most variance), vibration analysis (natural frequencies).

---

Fitting Models to Data

Curve fitting finds mathematical relationships in noisy data. Least squares minimizes the sum of squared residuals, providing optimal estimates under the assumption of Gaussian noise.

Least Squares Approximation

• The normal equations $A^T A \mathbf{c} = A^T \mathbf{b}$ solve for coefficients $\mathbf{c}$ that minimize $\|\mathbf{b} - A\mathbf{c}\|_2^2$. Here $A$ is the design matrix (e.g., columns of $1, x, x^2, \ldots$ for polynomial fitting).
• Overfitting vs. underfitting: complex models fit training data closely but generalize poorly; simple models may miss real patterns. The model degree is a key choice.
• QR factorization approach is numerically superior to forming $A^T A$ directly, because forming the normal equations squares the condition number ($\kappa(A^T A) = \kappa(A)^2$).

---

Solving Differential Equations Numerically

ODE solvers advance solutions through time when analytical solutions don't exist. The fundamental trade-off is between accuracy (higher-order methods) and stability (implicit methods for stiff problems).

Euler's Method

• First-order explicit method that updates via $y_{n+1} = y_n + h f(t_n, y_n)$. Simple to implement but requires small step sizes for accuracy.
• Local truncation error is $O(h^2)$, and global error is $O(h)$. Doubling the number of steps only halves the error.
• Valuable as an educational foundation for understanding more sophisticated methods, but rarely used in production code.

Runge-Kutta Methods

The classic fourth-order RK4 evaluates $f$ at four carefully chosen points per step:

1. $k_1 = f(t_n, y_n)$
2. $k_2 = f(t_n + h/2,\; y_n + h k_1/2)$
3. $k_3 = f(t_n + h/2,\; y_n + h k_2/2)$
4. $k_4 = f(t_n + h,\; y_n + h k_3)$
5. $y_{n+1} = y_n + \frac{h}{6}(k_1 + 2k_2 + 2k_3 + k_4)$

This achieves $O(h^5)$ local error ($O(h^4)$ global). Adaptive step-size control (e.g., RK45 embedded pairs) compares solutions at different orders to estimate error and adjust $h$ automatically.

RK4 is the workhorse of scientific computing for non-stiff problems, balancing accuracy, stability, and implementation complexity.

---

Discretizing Continuous Problems

Finite difference methods convert differential equations into algebraic systems by replacing derivatives with discrete approximations on a grid. Grid spacing, boundary conditions, and stability constraints determine solution quality.

Finite Difference Schemes

• Spatial discretization replaces continuous derivatives with algebraic expressions. For example, the second derivative becomes $\frac{\partial^2 u}{\partial x^2} \approx \frac{u_{i+1} - 2u_i + u_{i-1}}{(\Delta x)^2}$
• The CFL condition $\frac{c \,\Delta t}{\Delta x} \leq 1$ constrains the time step for explicit schemes solving hyperbolic PDEs (like the wave equation). Violating it causes the numerical solution to blow up.
• Implicit vs. explicit schemes: explicit methods are simple but only conditionally stable (must satisfy CFL-type constraints). Implicit methods require solving a linear system at each time step but allow much larger time steps, which is critical for stiff problems or diffusion equations.

---

Understanding Errors and Stability

Error analysis and stability assessment determine whether numerical results are trustworthy. Every computation accumulates two kinds of error: truncation error (from mathematical approximations) and round-off error (from finite-precision arithmetic).

Error Propagation and Stability

• Condition number $\kappa(A)$ measures how much input perturbations get amplified into output errors. For a linear system $Ax = b$, the relative error in $x$ can be up to $\kappa(A)$ times the relative error in $b$. A large condition number means the problem itself is sensitive, regardless of what algorithm you use.
• Stability is a property of the algorithm, not the problem. A stable method keeps rounding errors bounded as the computation proceeds. An unstable method can produce garbage regardless of step size.
• The Lax equivalence theorem (for linear PDEs): a consistent finite difference scheme converges to the true solution if and only if it is stable. Convergence = consistency + stability.

---

Quick Reference Table

---

Self-Check Questions

1. Convergence comparison: Why does Newton-Raphson converge faster than bisection near a simple root, and under what conditions might you prefer bisection anyway?

2. Method selection: You need to integrate a smooth function with high accuracy using minimal function evaluations. Would you choose the trapezoidal rule or Simpson's rule, and why does the error order matter?

3. Stability analysis: Explain why numerical differentiation is inherently less stable than numerical integration, and describe one strategy to mitigate this problem.

4. Compare and contrast: Both interpolation and least squares fitting produce functions from data points. When would you choose each approach, and what happens if you use interpolation on noisy experimental data?

5. Synthesis: A colleague proposes using Euler's method with very small time steps to solve a stiff ODE system. Explain why this approach is problematic and recommend an alternative, justifying your choice in terms of stability properties.