← Back to CSIR NET

Numerical Analysis

Algorithms and error analysis for root finding, interpolation, quadrature, and the numerical solution of differential equations and linear systems.

Root Finding Interpolation Numerical Integration Numerical ODE Iterative Linear Solvers

0 of 5 units completed0%

Root Finding Methods

Iterative methods for solving \(f(x)=0\) and their convergence analysis.

🎯 Bisection & Fixed Point Iteration

Theorem — Bisection Method Convergence If \(f\) is continuous on \([a,b]\) with \(f(a)f(b) < 0\), then the bisection method converges to a root. After \(n\) iterations, the error satisfies \(|e_n| \le \frac{b-a}{2^{n+1}}\). The convergence is linear with rate \(1/2\).

Theorem — Fixed Point Iteration If \(g: [a,b] \to [a,b]\) is a contraction with \(|g'(x)| \le L < 1\) for all \(x \in [a,b]\), then the iteration \(x_{n+1} = g(x_n)\) converges to the unique fixed point \(x^*\) for any starting point in \([a,b]\). The convergence is linear: \(|x_{n+1} - x^*| \le L|x_n - x^*|\).

⚡ Newton-Raphson & Secant Methods

Definition — Newton-Raphson Method \[x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}\] This is a fixed-point iteration with \(g(x) = x - f(x)/f'(x)\).

Theorem — Quadratic Convergence of Newton's Method If \(f \in C^2\), \(f(x^*)=0\), \(f'(x^*)\neq 0\), and \(x_0\) is sufficiently close to \(x^*\), then Newton's method converges quadratically: \[|x_{n+1} - x^*| \le \frac{|f''(\xi)|}{2|f'(x^*)|}\,|x_n - x^*|^2\] For a root of multiplicity \(m > 1\), convergence degrades to linear with rate \(1 - 1/m\). The modified Newton method \(x_{n+1} = x_n - m\,f(x_n)/f'(x_n)\) restores quadratic convergence.

The secant method replaces \(f'(x_n)\) with a finite difference: \(x_{n+1} = x_n - f(x_n)\frac{x_n - x_{n-1}}{f(x_n)-f(x_{n-1})}\). Its order of convergence is \(\phi = \frac{1+\sqrt{5}}{2} \approx 1.618\) (the golden ratio).

Example

Apply one step of Newton's method to \(f(x) = x^3 - 2\) starting at \(x_0 = 1.5\).

\(f(1.5) = 3.375 - 2 = 1.375\), \(f'(x) = 3x^2\), \(f'(1.5) = 6.75\). \[x_1 = 1.5 - \frac{1.375}{6.75} \approx 1.2963\] The exact root is \(\sqrt[3]{2} \approx 1.2599\), so one iteration brings us much closer.

Interpolation

Polynomial interpolation, divided differences, Hermite interpolation, and splines.

📈 Lagrange & Newton Interpolation

Definition — Lagrange Interpolating Polynomial Given \(n+1\) distinct points \((x_0,y_0),\ldots,(x_n,y_n)\), the unique polynomial of degree \(\le n\) passing through all points is: \[P_n(x) = \sum_{k=0}^{n} y_k \prod_{\substack{j=0\\j\neq k}}^{n}\frac{x - x_j}{x_k - x_j}\]

Theorem — Interpolation Error If \(f \in C^{n+1}[a,b]\) and \(P_n\) interpolates \(f\) at \(x_0,\ldots,x_n \in [a,b]\), then for each \(x\): \[f(x) - P_n(x) = \frac{f^{(n+1)}(\xi)}{(n+1)!}\prod_{k=0}^{n}(x - x_k)\] for some \(\xi\) between the extreme values of \(x, x_0,\ldots,x_n\).

Newton's divided difference form: \(P_n(x) = f[x_0] + f[x_0,x_1](x-x_0) + \cdots + f[x_0,\ldots,x_n]\prod_{k=0}^{n-1}(x-x_k)\), where divided differences are defined recursively: \(f[x_i,\ldots,x_{i+k}] = \frac{f[x_{i+1},\ldots,x_{i+k}] - f[x_i,\ldots,x_{i+k-1}]}{x_{i+k} - x_i}\).

🔗 Hermite Interpolation & Splines

Definition — Hermite Interpolation Hermite interpolation matches both function values and derivatives at each node. Given \(\{(x_i, f(x_i), f'(x_i))\}_{i=0}^{n}\), the Hermite interpolant has degree \(\le 2n+1\).

Definition — Cubic Spline A cubic spline \(S\) on \([a,b]\) with knots \(a = x_0 < x_1 < \cdots < x_n = b\) is a piecewise cubic polynomial \(S_i\) on each \([x_i, x_{i+1}]\) satisfying \(S \in C^2[a,b]\) and \(S(x_i) = f(x_i)\). The natural spline has \(S''(a) = S''(b) = 0\); the clamped spline matches \(S'(a)=f'(a),\; S'(b) = f'(b)\).

Theorem — Spline Error Bound For a clamped cubic spline interpolating \(f \in C^4[a,b]\): \[\|f - S\|_\infty \le \frac{5}{384}\,h^4\|f^{(4)}\|_\infty\] where \(h = \max_i(x_{i+1}-x_i)\). This is \(O(h^4)\) — far better than high-degree polynomial interpolation for smooth functions.

Example

The Runge phenomenon: why does interpolating \(f(x) = \frac{1}{1+25x^2}\) on equispaced nodes fail?

On equispaced nodes in \([-1,1]\), the interpolation error grows exponentially with degree near the endpoints due to large oscillations of the node polynomial \(\prod(x-x_k)\). The remedy: use Chebyshev nodes \(x_k = \cos\!\left(\frac{2k+1}{2(n+1)}\pi\right)\), which minimize the maximum of \(|\prod(x-x_k)|\) and eliminate the Runge phenomenon.

Numerical Integration

Composite quadrature rules, Gaussian quadrature, and error analysis for numerical integration.

📊 Newton-Cotes & Composite Rules

Theorem — Error Bounds for Newton-Cotes Rules

Trapezoidal rule: \(\int_a^b f\,dx = \frac{h}{2}(f_0+f_n) + h\sum_{i=1}^{n-1}f_i - \frac{(b-a)h^2}{12}f''(\xi)\). Error: \(O(h^2)\).
Simpson's rule (composite): Error is \(-\frac{(b-a)h^4}{180}f^{(4)}(\xi)\). Order: \(O(h^4)\).
Simpson's 3/8 rule: Uses cubic interpolation; error is also \(O(h^4)\) but with a larger constant.

The degree of precision of a quadrature rule is the highest degree polynomial it integrates exactly. Trapezoidal has degree 1; Simpson's has degree 3 (one more than expected, due to symmetry).

⭐ Gaussian Quadrature

Theorem — Gaussian Quadrature Optimality An \(n\)-point Gaussian quadrature rule \[\int_{-1}^{1}f(x)\,dx \approx \sum_{i=1}^{n}w_i f(x_i)\] has degree of precision \(2n-1\) when the nodes \(x_i\) are the roots of the \(n\)-th Legendre polynomial \(P_n(x)\). No quadrature rule with \(n\) nodes can achieve higher precision.

Example

Compute \(\int_{-1}^{1}(3x^4 - x^2 + 2)\,dx\) using 2-point Gaussian quadrature. Is the result exact?

Two-point Gauss-Legendre: nodes \(x_{1,2} = \pm 1/\sqrt{3}\), weights \(w_1 = w_2 = 1\). \[I \approx f(-1/\sqrt{3}) + f(1/\sqrt{3}) = 2\left(\frac{3}{9} - \frac{1}{3} + 2\right) = 2\cdot 2 = 4\] Exact value: \(\left[\frac{3x^5}{5} - \frac{x^3}{3} + 2x\right]_{-1}^{1} = 2(\frac{3}{5}-\frac{1}{3}+2)=\frac{68}{15}\approx 4.533\). Not exact — the integrand is degree 4, but 2-point Gauss has precision 3.

Numerical Methods for ODEs

Single-step and multi-step methods for initial value problems, with convergence and stability analysis.

🔢 Euler & Runge-Kutta Methods

Definition — Euler Methods For \(y' = f(t,y)\), \(y(t_0) = y_0\) with step size \(h\):

Forward Euler: \(y_{n+1} = y_n + hf(t_n, y_n)\) — explicit, \(O(h)\) error
Backward Euler: \(y_{n+1} = y_n + hf(t_{n+1}, y_{n+1})\) — implicit, \(O(h)\), A-stable

Theorem — Classical RK4 The 4th-order Runge-Kutta method: \[y_{n+1} = y_n + \frac{h}{6}(k_1 + 2k_2 + 2k_3 + k_4)\] where \(k_1 = f(t_n, y_n)\), \(k_2 = f(t_n+h/2,\; y_n+hk_1/2)\), \(k_3 = f(t_n+h/2,\; y_n+hk_2/2)\), \(k_4 = f(t_n+h,\; y_n+hk_3)\). The local truncation error is \(O(h^5)\), giving a global error of \(O(h^4)\).

🔄 Multi-Step Methods & Stability

Definition — Linear Multi-Step Methods A \(k\)-step method has the form: \(\sum_{j=0}^{k}\alpha_j y_{n+j} = h\sum_{j=0}^{k}\beta_j f_{n+j}\). Adams-Bashforth methods have \(\beta_k = 0\) (explicit); Adams-Moulton have \(\beta_k \neq 0\) (implicit).

Theorem — Dahlquist Equivalence Theorem A linear multi-step method is convergent if and only if it is consistent (order \(\ge 1\)) and zero-stable (roots of the characteristic polynomial \(\rho(z) = \sum \alpha_j z^j\) lie in or on the unit circle, with those on the circle being simple).

Theorem — Dahlquist Barriers

The order of a stable \(k\)-step method is at most \(k+1\) (\(k\) odd) or \(k+2\) (\(k\) even).
No explicit linear multi-step method is A-stable.
Among A-stable implicit methods, the highest order is 2 (achieved by the trapezoidal rule).

Example

Determine the region of absolute stability for forward Euler applied to \(y' = \lambda y\).

Applying forward Euler: \(y_{n+1} = y_n + h\lambda y_n = (1 + h\lambda)y_n\). For stability, \(|1 + h\lambda| \le 1\). With \(z = h\lambda \in \mathbb{C}\), the stability region is the disk \(|1+z| \le 1\), centered at \((-1,0)\) with radius 1. For real \(\lambda < 0\), stability requires \(0 < h < 2/|\lambda|\).

Iterative Methods for Linear Systems

Jacobi, Gauss-Seidel, SOR, and their convergence criteria for solving \(Ax = b\).

🔧 Splitting Methods

Write \(A = D - L - U\) where \(D\) is diagonal, \(L\) is strict lower triangular, \(U\) is strict upper triangular.

Definition — Jacobi & Gauss-Seidel Iterations

Jacobi: \(x^{(k+1)} = D^{-1}(L+U)x^{(k)} + D^{-1}b\). Iteration matrix: \(T_J = D^{-1}(L+U)\).
Gauss-Seidel: \(x^{(k+1)} = (D-L)^{-1}Ux^{(k)} + (D-L)^{-1}b\). Iteration matrix: \(T_{GS} = (D-L)^{-1}U\).
SOR: \(x^{(k+1)} = (D-\omega L)^{-1}[(1-\omega)D + \omega U]x^{(k)} + \omega(D-\omega L)^{-1}b\).

Theorem — Convergence Criterion An iterative method \(x^{(k+1)} = Tx^{(k)} + c\) converges for any initial guess if and only if the spectral radius \(\rho(T) < 1\).

Theorem — Sufficient Conditions for Convergence

If \(A\) is strictly diagonally dominant (\(|a_{ii}| > \sum_{j\neq i}|a_{ij}|\) for all \(i\)), both Jacobi and Gauss-Seidel converge.
If \(A\) is symmetric positive definite, Gauss-Seidel converges.
For SOR: convergence requires \(0 < \omega < 2\). For SPD tridiagonal systems, optimal \(\omega = \frac{2}{1+\sqrt{1-\rho(T_J)^2}}\).

Example

For \(A = \begin{pmatrix}4&1\\1&3\end{pmatrix}\), compute the Jacobi iteration matrix and determine convergence.

\(D = \begin{pmatrix}4&0\\0&3\end{pmatrix}\), \(L+U = \begin{pmatrix}0&1\\1&0\end{pmatrix}\). So \(T_J = D^{-1}(L+U) = \begin{pmatrix}0&1/4\\1/3&0\end{pmatrix}\). Eigenvalues: \(\lambda^2 = 1/12\), so \(\lambda = \pm 1/\sqrt{12}\). Spectral radius \(\rho(T_J) = 1/\sqrt{12} < 1\). Jacobi converges. (Note: \(A\) is strictly diagonally dominant, confirming convergence.)

Key Takeaways

Newton's method converges quadratically for simple roots; the secant method achieves order ~1.618 without requiring derivatives.
High-degree polynomial interpolation on equispaced nodes suffers from the Runge phenomenon — use Chebyshev nodes or splines.
Gaussian quadrature with \(n\) nodes achieves precision \(2n-1\) — optimal among all \(n\)-point rules.
RK4 achieves \(O(h^4)\) global accuracy and is the workhorse single-step method.
Dahlquist's theorems set fundamental limits on accuracy and stability of multi-step methods.
Spectral radius of the iteration matrix determines convergence of iterative linear solvers.

Practice Problems

Problem 1

Newton's method applied to \(f(x) = x^2\) from \(x_0 = 1\). Find \(x_n\) in terms of \(n\) and determine the order of convergence.

Show Solution ▼

\(x_{n+1} = x_n - x_n^2/(2x_n) = x_n/2\). So \(x_n = 1/2^n \to 0\). The convergence is only linear (with rate 1/2) because \(x^*=0\) is a double root. This illustrates the degradation for multiple roots.

Problem 2

Find the Lagrange interpolating polynomial through \((0,1), (1,2), (2,5)\).

Show Solution ▼

\(L_0 = \frac{(x-1)(x-2)}{(0-1)(0-2)} = \frac{(x-1)(x-2)}{2}\), \(L_1 = \frac{x(x-2)}{(1)(1-2)} = -x(x-2)\), \(L_2 = \frac{x(x-1)}{(2)(1)} = \frac{x(x-1)}{2}\). \(P_2(x) = \frac{(x-1)(x-2)}{2} - 2x(x-2) + \frac{5x(x-1)}{2} = x^2 + 1\).

Problem 3

How many nodes does Gaussian quadrature need to integrate \(p(x) = 6x^5 - 2x^3 + x\) exactly on \([-1,1]\)?

Show Solution ▼

\(p\) has degree 5. An \(n\)-point Gauss rule has precision \(2n-1\). We need \(2n-1 \ge 5\), so \(n \ge 3\). 3 Gauss-Legendre nodes suffice.

Problem 4

Apply one step of forward Euler with \(h=0.1\) to \(y' = -2y,\; y(0) = 1\). Compare with the exact solution.

Show Solution ▼

\(y_1 = 1 + 0.1(-2)(1) = 0.8\). Exact: \(y(0.1) = e^{-0.2} \approx 0.8187\). Error \(\approx 0.0187\). The stability condition \(h < 2/|\lambda| = 1\) is satisfied.

Problem 5

Show that Gauss-Seidel diverges for \(A = \begin{pmatrix}1&2\\2&1\end{pmatrix}\).

Show Solution ▼

\(D = I\), \(L = \begin{pmatrix}0&0\\2&0\end{pmatrix}\), \(U = \begin{pmatrix}0&2\\0&0\end{pmatrix}\). \(T_{GS} = (D-L)^{-1}U = \begin{pmatrix}1&0\\-2&1\end{pmatrix}^{-1}\begin{pmatrix}0&2\\0&0\end{pmatrix} = \begin{pmatrix}0&2\\0&4\end{pmatrix}\). Eigenvalues: 0 and 4. Since \(\rho(T_{GS})=4 > 1\), Gauss-Seidel diverges. Note: \(A\) is not diagonally dominant.

Self-Assessment Quiz

1. The order of convergence of Newton's method for a simple root is:

A 1 (linear)

B 2 (quadratic)

C 3 (cubic)

D \(\phi \approx 1.618\)

2. The Runge phenomenon is resolved by:

A Increasing the polynomial degree

B Using Chebyshev nodes

C Reducing the number of nodes

D Extrapolation

3. The degree of precision of \(n\)-point Gaussian quadrature is:

A \(n-1\)

B \(2n-2\)

C \(2n-1\)

D \(2n\)

4. According to Dahlquist's second barrier, the highest-order A-stable linear multi-step method has order:

A 1

B 2

C 3

D 4

5. An iterative method \(x^{(k+1)} = Tx^{(k)} + c\) converges for all initial guesses iff:

A \(\|T\| < 1\) for some matrix norm

B \(\rho(T) < 1\)

C \(T\) is symmetric

D \(\det(T) < 1\)