Numerical methods

In short

Numerical methods replace exact algebra by finite, repeatable arithmetic: a root is reached by iteration, an integral by a weighted sum of ordinates, a differential equation by stepping forward in small increments. They exist because most equations met in practice — transcendental equations, large linear systems, integrals with no closed form, non-linear ODEs — have no analytic solution. In a B.Sc (CBCS) course the topic carries a full paper plus a practical lab in C/Python/Scilab; it is examined in GATE MA (and in GATE CS/ME as "numerical methods"), CSIR-UGC NET, IIT JAM-based M.Sc entrances and in state PSC technical papers. In jobs it is the mathematics underneath every simulation and analytics product — FEA/CFD at Ansys, Altair, Siemens and L&T, trajectory and orbit computation at ISRO and DRDO, weather models at IMD/NCMRWF, pricing and risk engines at banks and actuarial firms, and SciPy/NumPy/MATLAB work in any data-science role.

Explained

Errors, significant digits and conditioning

• Absolute error Ea = |true − approx|; relative error Er = Ea/|true|; percentage error = 100 Er. Relative error is what matters, because it is scale-free.
• Round-off error comes from storing a number in finitely many digits (IEEE 754 double ≈ 16 significant decimal digits); truncation error comes from cutting an infinite process short (stopping a Taylor series, using finite h).
• An approximation is correct to n significant digits if Er ≤ 0.5 × 10¹⁻ⁿ. Rounding to n digits gives Er ≤ 0.5 × 10¹⁻ⁿ.
• Error propagation: in a sum, absolute errors add; in a product or quotient, relative errors add; for y = f(x), Δy ≈ |f′(x)| Δx.
• Catastrophic cancellation: subtracting two nearly equal numbers destroys significant digits. Compute √(x+1) − √x as 1/(√(x+1)+√x); solve a quadratic with the stable root formula.
• Ill-conditioned problem: a small change in data causes a large change in the answer (condition number large). Unstable algorithm: errors grow as the computation proceeds. A well-conditioned problem can still be ruined by an unstable method.

Roots of algebraic and transcendental equations

• All methods need a bracket or a starting guess. By the intermediate value theorem, if f(a)·f(b) < 0 and f is continuous, a root lies in (a, b).
• Bisection: repeatedly halve the bracket, keeping the half where the sign changes. Always converges, linearly (error halves each step); after n steps error ≤ (b−a)/2ⁿ, so n ≥ log₂((b−a)/ε).
• Regula falsi (false position): replace f by the chord joining (a, f(a)) and (b, f(b)); the bracket is retained, so convergence is guaranteed but can be slow and one-sided.
• Newton–Raphson: x_(n+1) = x_n − f(x_n)/f′(x_n) — the tangent method. Quadratic convergence (digits roughly double) near a simple root. Fails or diverges if f′(x_n) ≈ 0, at a point of inflexion, if the guess is poor, or oscillates between two points; at a root of multiplicity m convergence degrades to linear (cure: x − m f/f′).
• Secant method: Newton with the derivative replaced by a difference quotient; needs no f′, order ≈ 1.618 (golden ratio), but no bracket, so it can diverge.
• Fixed-point (successive substitution): write f(x)=0 as x = g(x) and iterate x_(n+1)=g(x_n). Converges in an interval containing the root if |g′(x)| < 1 there; smaller |g′| means faster convergence.
• Order of convergence p: |e_(n+1)| ≈ C|e_n|^p. Bisection p=1 (C=½), regula falsi p=1, secant p≈1.618, Newton p=2.
• Special purpose: Horner's scheme to evaluate a polynomial cheaply; synthetic division and Birge–Vieta for polynomial roots; Graeffe's root-squaring and Bairstow (complex pairs) appear in some syllabi.

Finite differences, interpolation and curve fitting

• Operators: forward Δf(x)=f(x+h)−f(x); backward ∇f(x)=f(x)−f(x−h); central δ; shift E f(x)=f(x+h), with Δ = E − 1 and the relation (1+Δ)(1−∇)=1.
• Newton's forward difference formula suits interpolation near the start of an equally spaced table; Newton's backward near the end; central difference formulas (Gauss, Stirling, Bessel) near the middle.
• Newton's divided difference and Lagrange's formula handle unequally spaced data. Divided differences are symmetric in the arguments and allow a new point to be added without redoing the work; Lagrange does not.
• Through n+1 distinct points there is exactly one polynomial of degree ≤ n — all these formulas are rearrangements of the same polynomial.
• Error: f(x) − P_n(x) = fⁿ⁺¹(ξ)·Π(x−x_i)/(n+1)!. High-degree interpolation over a wide range oscillates badly (Runge phenomenon); prefer low-degree piecewise or cubic spline interpolation, which matches value, slope and curvature at each knot.
• Inverse interpolation: find x for a given y, e.g. by Lagrange with x and y interchanged.
• Extrapolation (outside the table) is far less reliable than interpolation — error grows rapidly.
• Least squares curve fitting is different in kind: the curve need not pass through the points, it minimises Σ(residual)². Used for experimental data with scatter; interpolation is used for exact tabulated values.

Numerical differentiation and integration

• Differentiation from a difference table: forward/backward differences give O(h) accuracy, the central difference (f(x+h)−f(x−h))/2h gives O(h²). Differentiation is noise-amplifying: too small an h makes round-off dominate, so an optimum h exists.
• Newton–Cotes rules integrate the interpolating polynomial over equal subintervals:
  • Trapezoidal rule — straight lines, error O(h²), exact for degree 1, needs any number of intervals.
  • Simpson's 1/3 rule — parabolas, error O(h⁴), exact for degree 3, needs an even number of intervals.
  • Simpson's 3/8 rule — cubics, error O(h⁴), number of intervals a multiple of 3.
  • Boole's and Weddle's rules — higher order, intervals in multiples of 4 and 6 respectively.
• Romberg integration: systematically combine trapezoidal values at h and h/2 to cancel error terms, R = (4 T(h/2) − T(h))/3 (which is Simpson's value).
• Gauss–Legendre quadrature: choose both weights and nodes; an n-point rule is exact for polynomials of degree 2n−1 — far more accurate per evaluation than Newton–Cotes, but the nodes are irrational and not equally spaced, so it cannot be used on a fixed data table.
• Double integrals are handled by applying a 1-D rule in each direction (trapezoidal or Simpson, repeatedly).

Systems of linear equations and eigenvalues

• Direct methods finish in a fixed number of operations:
  • Gauss elimination — forward elimination to upper triangular form, then back substitution; about 2n³/3 operations.
  • Partial pivoting — swap in the largest available pivot; essential for stability, and mandatory when a pivot is zero.
  • Gauss–Jordan — reduce to the identity; gives the inverse, costs more than elimination.
  • LU (Doolittle/Crout) decomposition — A = LU once, then solve cheaply for many right-hand sides; Cholesky A = LLᵀ for symmetric positive-definite A.
  • Thomas algorithm for tridiagonal systems, O(n) — used by spline and finite-difference codes.
• Iterative methods suit large sparse systems:
  • Jacobi uses only previous-iteration values; Gauss–Seidel uses each new value immediately and is usually about twice as fast; SOR adds a relaxation factor ω (1<ω<2).
  • Sufficient condition for convergence: strict diagonal dominance, |a_ii| > Σ_{j≠i}|a_ij| for every row (also convergent for symmetric positive-definite A). Rows may be rearranged to achieve it.
• Eigenvalues: the power method converges to the largest-magnitude eigenvalue and its vector; the inverse power method gives the smallest; Jacobi's rotation method gives all eigenvalues of a symmetric matrix; Gerschgorin's circles bound them.
• An ill-conditioned system (near-singular A, large condition number ‖A‖‖A⁻¹‖) gives wildly different answers for tiny data changes, however good the method.

Ordinary differential equations and initial-value problems

• Problem: dy/dx = f(x, y), y(x₀) = y₀; march forward in steps of h.
• Taylor series method: use the expansion of y about x₀; accurate but needs successive derivatives of f.
• Euler's method y_(n+1) = y_n + h f(x_n, y_n): first order, O(h) global error, cheap but inaccurate.
• Modified Euler / Heun (improved Euler): predict with Euler, correct with the average slope — second order.
• Runge–Kutta methods get high order without derivatives: RK2 is Heun's; RK4 needs four function evaluations per step and is fourth order — the workhorse, self-starting, easy to code.
• Predictor–corrector (multistep) methods reuse past values: Milne–Simpson and Adams–Bashforth/Adams–Moulton. One function evaluation per step (cheaper than RK4) but they are not self-starting — RK4 supplies the first few values — and they need equal steps.
• Single-step vs multistep: single-step methods allow easy step-size change; multistep methods are cheaper per step and give a built-in error estimate (predictor − corrector difference).
• Stability and stiffness: explicit methods need h small enough to stay in the stability region; stiff equations (very different time scales) force absurdly small h, so implicit methods (backward Euler, trapezoidal/Crank–Nicolson) are used.
• Boundary-value problems: finite-difference replacement of derivatives giving a tridiagonal system, or the shooting method (guess the initial slope, iterate on it). Picard's method gives a successive-approximation analytic solution.

Key points

• Truncation error – from approximating the process (reduce by smaller h or higher order); round-off error – from finite precision (grows as h shrinks).
• Bisection – always converges, linear, error ≤ (b−a)/2ⁿ; slow but safe.
• Newton–Raphson – quadratic convergence, needs f′ and a good guess; fails when f′ ≈ 0.
• Secant – order 1.618, no derivative needed, no bracket guarantee.
• Fixed-point convergence – requires |g′(x)| < 1 near the root.
• Interpolating polynomial – unique for n+1 distinct points; Newton and Lagrange give the same polynomial.
• Divided differences – symmetric and incremental; new data points cost almost nothing.
• Runge phenomenon – high-degree interpolation oscillates near the ends; use splines instead.
• Trapezoidal rule – error O(h²), any n; Simpson's 1/3 – O(h⁴), n even; Simpson's 3/8 – n a multiple of 3.
• Gauss quadrature – n points, exact up to degree 2n−1; needs freely chosen nodes.
• Partial pivoting – required for numerical stability in Gauss elimination, not just for zero pivots.
• LU decomposition – factor once, solve many right-hand sides cheaply.
• Diagonal dominance – sufficient condition for Jacobi and Gauss–Seidel to converge.
• Gauss–Seidel vs Jacobi – Seidel uses updated values at once and typically converges in about half the iterations.
• Euler O(h), RK4 O(h⁴) – RK4 with four evaluations beats Euler with four small steps by a huge margin.
• Predictor–corrector – cheaper per step than RK4 but not self-starting and fixed-step.
• Ill-conditioning – a property of the problem; instability – a property of the algorithm.
• Order of accuracy p – halving h reduces error by about 2^p; verify a code by checking this ratio.

Formulas and examples

Symbols: h = step size or tabular interval (same units as x); n = number of subintervals or data points; p = interpolation parameter (dimensionless); ξ = some point in the interval; errors carry the units of the quantity computed.

Example 1 — Newton–Raphson for x³ − 2x − 5 = 0, x₀ = 2.

• f′(x) = 3x² − 2. f(2) = 8 − 4 − 5 = −1, f′(2) = 10 → x₁ = 2 − (−1)/10 = 2.1.
• f(2.1) = 9.261 − 4.2 − 5 = 0.061, f′(2.1) = 11.23 → x₂ = 2.1 − 0.061/11.23 = 2.094568.
• f(2.094568) = 9.189320 − 4.189136 − 5 = 0.000184, f′ = 11.161645 → x₃ = 2.0945515.
• Correct digits: 2 → 5 → 8 — the signature of quadratic convergence. Bisection on [2, 3] would need about 24 steps for the same accuracy.

Example 2 — ∫₀¹ dx/(1+x²) with h = 0.25 (true value π/4 = 0.78539816). Ordinates: y₀ = 1, y₁ = 0.9411765, y₂ = 0.8, y₃ = 0.64, y₄ = 0.5.

• Trapezoidal: (0.25/2)[1 + 0.5 + 2(0.9411765 + 0.8 + 0.64)] = 0.125 × 6.262353 = 0.7827941 (error 2.6 × 10⁻³).
• Simpson's 1/3: (0.25/3)[1 + 0.5 + 4(0.9411765 + 0.64) + 2(0.8)] = 0.0833333 × 9.4247059 = 0.7853922 (error 6.0 × 10⁻⁶).
• Hence π ≈ 4 × 0.7853922 = 3.141569 from only five ordinates; Simpson is about 430 times more accurate here for the same data.

Example 3 — Newton's forward difference interpolation. Table x = 1, 2, 3, 4 with y = 1, 8, 27, 64. Differences: Δy = 7, 19, 37; Δ²y = 12, 18; Δ³y = 6. Estimate y(1.5): p = (1.5 − 1)/1 = 0.5. y = 1 + 0.5(7) + [0.5(−0.5)/2](12) + [0.5(−0.5)(−1.5)/6](6) = 1 + 3.5 − 1.5 + 0.375 = 3.375. Exact 1.5³ = 3.375 — a cubic is reproduced exactly because the third differences are constant and the fourth vanish.

Example 4 — Gauss–Seidel, 20x + y − 2z = 17, 3x + 20y − z = −18, 2x − 3y + 20z = 25. Each row is strictly diagonally dominant (20 > 1+2, 20 > 3+1, 20 > 2+3), so convergence is assured. Starting from (0, 0, 0):

Round 1: x = (17 − 0 + 0)/20 = 0.85; y = (−18 − 3(0.85) + 0)/20 = −1.0275; z = (25 − 2(0.85) + 3(−1.0275))/20 = 1.010875. Three sweeps give x = 1, y = −1, z = 1 to four decimals.

Example 5 — RK4 for dy/dx = x + y, y(0) = 1, h = 0.1.

• k₁ = 0.1(0 + 1) = 0.1; k₂ = 0.1(0.05 + 1.05) = 0.11; k₃ = 0.1(0.05 + 1.055) = 0.1105; k₄ = 0.1(0.1 + 1.1105) = 0.12105.
• Δy = (0.1 + 0.22 + 0.221 + 0.12105)/6 = 0.66205/6 = 0.1103417 → y(0.1) = 1.1103417.
• Exact y = 2eˣ − x − 1 = 1.1103418, so the error is ≈ 1 × 10⁻⁷. Euler gives y(0.1) = 1 + 0.1(1) = 1.1, error 1.03 × 10⁻² — about 10⁵ times worse.

# RK4 in four lines (NumPy/SciPy style used in labs and industry)
def rk4(f, x, y, h, n):
    for _ in range(n):
        k1 = h*f(x, y);        k2 = h*f(x+h/2, y+k1/2)
        k3 = h*f(x+h/2, y+k2/2); k4 = h*f(x+h, y+k3)
        y += (k1 + 2*k2 + 2*k3 + k4)/6; x += h
    return y

Where it is used

• Simulation and CAE: finite-element and CFD solvers (Ansys, Altair, Siemens, Dassault centres in Bengaluru, Pune and Hyderabad; L&T and Tata Technologies services teams) are built on sparse linear solvers, quadrature and time-stepping schemes.
• Space and defence: ISRO/VSSC orbit propagation, re-entry trajectories and launch-vehicle dynamics integrate ODEs numerically; DRDO labs use the same methods for flight simulation.
• Weather and climate: IMD and NCMRWF run finite-difference and spectral models on HPC systems; numerical stability limits the permissible time step.
• Power, water and civil systems: load-flow analysis in POSOCO (Grid-India)/state SLDC practice solves large non-linear systems by Newton–Raphson; pipe-network and structural analysis use iterative sparse solvers.
• Finance and actuarial work: IRR and bond-yield calculations are root-finding problems; option pricing uses finite-difference and quadrature methods at bank quant desks and GIFT City firms.
• Data science and engineering roles: gradient descent, numerical optimisation, spline smoothing and regression fitting in NumPy/SciPy, MATLAB/Octave, R and Scilab — the starting point for analyst and scientific-computing jobs, and the basis of placement coding questions on bisection, Simpson and RK4.

Interview questions

• Q: Why would you use bisection when Newton–Raphson is much faster? A: Bisection cannot fail once a sign change is bracketed, needs no derivative, and gives a guaranteed error bound; Newton can diverge from a bad guess. In practice a hybrid (bisection to localise, then Newton) is used.
• Q: What is the order of convergence of Newton–Raphson and when is it lost? A: Quadratic for a simple root; it degrades to linear at a multiple root and the method fails where f′ ≈ 0 or the initial guess is far from the root.
• Q: Difference between interpolation and curve fitting? A: Interpolation forces the curve through every data point and is for exact tabulated values; least-squares fitting minimises squared residuals and is for scattered experimental data.
• Q: Why must Simpson's 1/3 rule have an even number of intervals? A: It fits one parabola over each pair of subintervals, so the subintervals must pair up exactly.
• Q: Simpson's rule is built from quadratics — why is it exact for cubics? A: The cubic part of the error integrates to zero over the symmetric pair of intervals, so the leading error term involves the fourth derivative.
• Q: Why is pivoting needed in Gauss elimination? A: A zero pivot breaks the algorithm and a small pivot amplifies round-off error; partial pivoting picks the largest available pivot and keeps the solution stable.
• Q: When would you choose an iterative solver over Gauss elimination? A: For large sparse systems — Gauss–Seidel keeps only the non-zeros, costs O(non-zeros) per sweep and avoids the 2n³/3 cost and fill-in of elimination.
• Q: What condition guarantees Gauss–Seidel converges? A: Strict diagonal dominance in every row (or a symmetric positive-definite matrix); rows may be reordered to obtain it.
• Q: Why does reducing h not always reduce error? A: Truncation error falls with h, but round-off error rises (especially in differentiation, where you divide by h), so total error passes through a minimum.
• Q: Difference between an ill-conditioned problem and an unstable method? A: Ill-conditioning is inherent to the problem — small data changes cause large answer changes; instability is the method's fault, errors amplifying as it proceeds.
• Q: RK4 uses four evaluations per step. Why not take four Euler steps instead? A: RK4's global error is O(h⁴) against Euler's O(h); in the worked example RK4 was about 10⁵ times more accurate for the same number of evaluations.
• Q: What is a stiff ODE and how is it handled? A: One with widely separated time scales, where an explicit method needs an impractically small h for stability; implicit methods such as backward Euler or Crank–Nicolson are used.
• Q: Why is high-degree polynomial interpolation discouraged? A: The Runge phenomenon — the polynomial oscillates violently near the ends of the range; piecewise cubic splines are better behaved.

Common mistakes

• Using Simpson's 1/3 rule with an odd number of intervals, or 3/8 with a number not divisible by 3, and still quoting the error order.
• Taking h as the number of intervals instead of (b−a)/n, or forgetting to convert degrees to radians before evaluating trigonometric ordinates.
• Applying Newton–Raphson without checking f′ at the guess, then reporting a divergent or oscillating sequence as the root.
• Stopping iteration on |x_(n+1) − x_n| alone when the function is flat — check both the step and |f(x)|, and prefer a relative tolerance.
• Running Gauss–Seidel on a non-diagonally-dominant system without rearranging the rows, then concluding the method "does not work".
• Mixing up Jacobi and Gauss–Seidel by substituting newly computed values inside a Jacobi sweep (or old ones inside Seidel).
• Rounding intermediate results to two or three decimals: in exam work keep 5–6 decimals throughout, since iterative error accumulates.
• Extrapolating with an interpolation formula, or using a forward-difference formula near the end of a table where the required differences do not exist.

Related topics

• Calculus: limits, derivatives and Taylor series · Linear algebra: matrices, rank and eigenvalues · Ordinary differential equations · Linear programming · Probability and statistics: regression and least squares · Programming in C/Python for scientific computing