Numerical Analysis

Syllabus:

Interpolation: Lagrangian Interpolating Formula, Divided Differences and Their Properties, Newton's Fundamental Interpolation Formula, Finite Differences, Forward Differences, Backward Differences, relation between Forward and Backward Differences, Newton's Forward Difference Interpolation Formula.

Solution of Equations: Method of Bisections, Regula-Falsi Method, Newton-Raphson's Method.

Linear Equation: Gauss Elimination Method, Cholesky's Method, Gauss-Seidel Method, Method of Relaxation.

Ordinary Differential Equations: Euler's Method, Single Step Method, Runge-Kutta Method, Milne-Simpson Method.

In this notes, we'll cover -
Part 1. Interpolation

• Lagrange Interpolating Formula
• Divided Differences and Their Properties
• Newton's Divided Difference Formula
• Finite Differences (Forward, Backward)
• Relation Between Forward and Backward Differences
• Newton's Forward Difference Formula

Part 2. Solution of Equations

• Bisection Method
• Regula-Falsi Method
• Newton-Raphson Method

Part 3. Systems of Linear Equations

• Gauss Elimination
• Cholesky's Method
• Gauss-Seidel Method
• Method of Relaxation

Part 4. Ordinary Differential Equations

• Euler's Method
• Single-Step Methods
• Runge-Kutta Method
• Milne-Simpson Method

Part 1. Interpolation

The Interpolation Problem

Definition 1. Interpolation:
Let $y = f(x)$ be a function, and suppose the values of $f$ are known at distinct points
$x_0,\ x_1,\ \ldots,\ x_n$
as $y_0 = f(x_0),\ y_1 = f(x_1),\ \ldots,\ y_n = f(x_n)$.

Interpolation means finding an approximate value of $f(x)$ at some intermediate point $x$ lying between $x_0$ and $x_n$, by constructing a polynomial (or other simple function) that passes through the given data points.

The polynomial $P(x)$ of degree at most $n$ satisfying
$P(x_i) = y_i \quad (i = 0,1,\ldots,n)$
is called an interpolating polynomial.

Lagrange Interpolating Formula

Definition 2. Lagrange Fundamental Polynomials:
Given distinct nodes $x_0, x_1, \ldots, x_n$, the $i$-th Lagrange basis polynomial is

Clearly $\ell_i(x_k) = 0$ if $k \neq i$, and $\ell_i(x_i) = 1$.

Theorem 1 (Lagrange Interpolation Formula)

Let $y_i = f(x_i)$ be given at $(n+1)$ distinct points $x_0, x_1, \ldots, x_n$. Then there exists a unique polynomial $P_n(x)$ of degree at most $n$ such that $P_n(x_i) = y_i$ for each $i$, and

Proof:

Let the data $(x_i, y_i)$, $i = 0,1,\ldots,n$, be given, with all $x_i$ distinct.
Define
$P_n(x) = \sum_{i=0}^{n} y_i\,\ell_i(x).$

Step 1. Degree of $P_n$.
Each $\ell_i(x)$ is a product of $n$ linear factors, so $\deg(\ell_i) = n$.
Hence $P_n$ is a polynomial of degree at most $n$.

Step 2. Interpolation conditions.
For a fixed index $k$,
$P_n(x_k) = \sum_{i=0}^{n} y_i\,\ell_i(x_k) = y_k\,\ell_k(x_k) = y_k \cdot 1 = y_k,$
since $\ell_i(x_k) = 0$ for $i \neq k$.
Thus $P_n(x_i) = y_i$ for all $i$.

Step 3. Uniqueness.
Suppose $Q_n(x)$ is another polynomial of degree at most $n$ with $Q_n(x_i) = y_i$.
Then $R(x) = P_n(x) - Q_n(x)$ has degree at most $n$ and has $(n+1)$ roots $x_0,\ldots,x_n$.
A nonzero polynomial of degree $\leq n$ cannot have more than $n$ roots.
Hence $R \equiv 0$, so $P_n = Q_n$.

Therefore the Lagrange formula gives the unique interpolating polynomial of degree $\leq n$.

Error term (statement): If $f$ is $(n+1)$ times continuously differentiable, then for some $\xi$ between the smallest and largest of $x, x_0, \ldots, x_n$,

E.g., Find $P(x)$ through $(0,1)$, $(1,3)$, $(2,7)$.

Divided Differences

Definition 3. Divided Differences:
Let $y_i = f(x_i)$. The zeroth divided difference is
$f[x_i] = y_i.$

The first divided difference is
$f[x_i, x_{i+1}] = \frac{f[x_{i+1}] - f[x_i]}{x_{i+1} - x_i} = \frac{y_{i+1} - y_i}{x_{i+1} - x_i}.$

The second divided difference is

In general,

Theorem 2 (Properties of Divided Differences)

(i) Divided differences are symmetric in their arguments: the value of $f[x_0,x_1,\ldots,x_k]$ does not depend on the order of $x_0,\ldots,x_k$.
(ii) If $f$ is a polynomial of degree $m$, then the $(m+1)$-th and higher divided differences vanish.
(iii) (Mean value form) If $f^{(k)}$ exists, then
$f[x_0,x_1,\ldots,x_k] = \frac{f^{(k)}(\xi)}{k!}$
for some $\xi$ in the interval containing $x_0,\ldots,x_k$.

Proof of (i) (sketch):
One shows by induction that

which is visibly symmetric in the $x$'s (and $y$'s attached to them).

Proof of (ii):
If $\deg(f) = m$, then $f^{(m+1)} \equiv 0$. By (iii), the $(m+1)$-th divided difference is zero.

Newton's Divided Difference Formula

Theorem 3 (Newton's Fundamental / Divided Difference Formula)

The unique interpolating polynomial of degree at most $n$ may be written as

Proof:

Let $P_n(x)$ be the interpolating polynomial. Write it in the Newton form

Putting $x = x_0$: $P_n(x_0) = a_0 = y_0 = f[x_0]$.

Putting $x = x_1$:
$y_1 = a_0 + a_1(x_1-x_0) \implies a_1 = \frac{y_1-y_0}{x_1-x_0} = f[x_0,x_1].$

Continuing inductively, one finds $a_k = f[x_0,x_1,\ldots,x_k]$.

Hence Newton's divided difference formula holds.
(It agrees with Lagrange at the nodes, so by uniqueness it is the interpolating polynomial.)

Finite Differences

Assumption: Nodes are equally spaced: $x_i = x_0 + ih$, $i = 0,1,2,\ldots$, with step $h > 0$. Write $y_i = f(x_i)$.

Definition 4. Forward Difference Operator $\Delta$:
$\Delta y_i = y_{i+1} - y_i.$
Higher differences:
$\Delta^2 y_i = \Delta(\Delta y_i) = \Delta y_{i+1} - \Delta y_i = y_{i+2} - 2y_{i+1} + y_i,$
and in general
$\Delta^{k} y_i = \Delta^{k-1} y_{i+1} - \Delta^{k-1} y_i.$

Definition 5. Backward Difference Operator $\nabla$:
$\nabla y_i = y_i - y_{i-1}.$
Higher: $\nabla^2 y_i = \nabla y_i - \nabla y_{i-1}$, and so on.

Definition 6. Shift Operator $E$:
$E y_i = y_{i+1}, \qquad E^k y_i = y_{i+k}.$

Theorem 4 (Operator Relations)

The operators $E$, $\Delta$, and $\nabla$ are related by
$E = 1 + \Delta, \qquad \nabla = 1 - E^{-1}, \qquad \Delta = E\nabla = \nabla E.$
Also,
$\Delta^n y_i = \nabla^n y_{i+n}.$

Proof:

$E y_i = y_{i+1} = y_i + (y_{i+1}-y_i) = y_i + \Delta y_i = (1+\Delta)y_i$.
Hence $E = 1 + \Delta$.

Further, $\nabla y_{i+1} = y_{i+1} - y_i = \Delta y_i$, so $\Delta y_i = \nabla y_{i+1} = \nabla E y_i$, giving $\Delta = \nabla E$.
Similarly $\Delta = E\nabla$.

Finally, $\Delta y_i = y_{i+1}-y_i = \nabla y_{i+1}$, and by induction $\Delta^n y_i = \nabla^n y_{i+n}$.

Theorem 5 (Relation Between Forward and Backward Differences)

At a fixed tabular point,
$\Delta^n y_0 = \nabla^n y_n,$
and more generally the forward differences starting at $y_i$ equal the backward differences ending at $y_{i+n}$.

Proof:
This is the special case $i = 0$ of $\Delta^n y_i = \nabla^n y_{i+n}$ from Theorem 4.

Newton's Forward Difference Formula

Theorem 6 (Newton-Gregory Forward Difference Formula)

Let $x = x_0 + uh$, where $u = \dfrac{x-x_0}{h}$. Then the interpolating polynomial is

Proof:

Start from Newton's divided difference formula with equal spacing $x_i = x_0 + ih$.

One shows that
$f[x_0,x_1,\ldots,x_k] = \frac{\Delta^k y_0}{k!\, h^k}.$

Also
$(x-x_0)(x-x_1)\cdots(x-x_{k-1}) = h^k\, u(u-1)\cdots(u-k+1).$

Substituting these into the divided difference formula yields Newton's forward formula.

When to use: Prefer this formula when $x$ lies near the beginning of the table (small $u$).

Newton's backward formula (for reference): With $u = \dfrac{x-x_n}{h}$ and backward differences at $y_n$,

Prefer when $x$ is near the end of the table.

Part 2. Solution of Equations

We seek a root of $f(x) = 0$, where $f$ is continuous (and, for Newton-Raphson, differentiable).

Bisection Method

Idea: If $f(a)$ and $f(b)$ have opposite signs, then by the Intermediate Value Theorem there is at least one root in $(a,b)$. Repeatedly bisect the interval and keep the half where the sign change persists.

Theorem 7 (Bisection Method — Algorithm and Convergence)

Let $f$ be continuous on $[a,b]$ with $f(a)f(b) < 0$. Define

where $[x_n^{(L)}, x_n^{(R)}]$ is the current bracket containing a sign change. Then the sequence of midpoints converges to a root $\xi \in (a,b)$, and after $n$ steps
$|x_n - \xi| \leq \frac{b-a}{2^n}.$

Proof / Derivation:

Let $f$ be continuous on $[a,b]$ with $f(a)f(b) < 0$.
By IVT, there exists $\xi \in (a,b)$ with $f(\xi) = 0$.

Algorithm (one step):
Set $c = \dfrac{a+b}{2}$.

• If $f(c) = 0$, stop; $c$ is a root.
• If $f(a)f(c) < 0$, replace $b$ by $c$.
• If $f(c)f(b) < 0$, replace $a$ by $c$.

The new interval still contains a root and has length half of the previous one.

After $n$ bisections, the interval length is $\dfrac{b-a}{2^n}$.
Any point $x_n$ of the current interval (in particular the midpoint) satisfies
$|x_n - \xi| \leq \frac{b-a}{2^n}.$
Hence $x_n \to \xi$ as $n \to \infty$.

Rate: Linear convergence (error roughly halves each time). Guaranteed, but slow.

E.g., Solve $x^3 - x - 1 = 0$ on $[1,2]$.
$f(1) = -1 < 0$, $f(2) = 5 > 0$.
$c_1 = 1.5$, $f(1.5) = 0.875 > 0$ $\to$ new interval $[1, 1.5]$.
Continue until the required accuracy.

Regula-Falsi Method (Method of False Position)

Idea: Instead of the midpoint, take the point where the chord joining $(a,f(a))$ and $(b,f(b))$ meets the $x$-axis.

Theorem 8 (Regula-Falsi Formula)

Derivation:

Let $f$ be continuous on $[a,b]$ with $f(a)f(b) < 0$.
The equation of the chord through $(a,f(a))$ and $(b,f(b))$ is

Set $y = 0$ to find the intercept $x_1$:

If $f(a)f(x_1) < 0$, set $b := x_1$; else set $a := x_1$. Repeat.

Note: Convergence is usually faster than bisection, but one endpoint may stay fixed for many iterations (linear convergence still).

Newton-Raphson Method

Idea: Replace $f$ near $x_n$ by its tangent line; take the next iterate as the root of the tangent.

Theorem 9 (Newton-Raphson Formula)

Let $f$ be continuously differentiable and $f'(x) \neq 0$ near a root. The Newton-Raphson iteration is

If $f(\xi)=0$, $f'(\xi)\neq 0$, and $f''$ is continuous near $\xi$, then for $x_0$ sufficiently close to $\xi$ the method converges quadratically:

Derivation:

Let $x_n$ be an approximation to the root.
The tangent to $y = f(x)$ at $x_n$ is
$y - f(x_n) = f'(x_n)\,(x - x_n).$
It meets the $x$-axis where $y = 0$:

Geometric meaning: Each step follows the tangent down (or up) to the axis.

Caution: Fails or diverges if $f'(x_n) \approx 0$, or if $x_0$ is poorly chosen. For multiple roots, convergence drops to linear unless a modified formula is used.

E.g., Find $\sqrt{2}$ by solving $f(x) = x^2 - 2 = 0$.
$x_{n+1} = x_n - \frac{x_n^2-2}{2x_n} = \frac{1}{2}\left(x_n + \frac{2}{x_n}\right).$
With $x_0 = 1.5$: $x_1 = 1.4167$, $x_2 = 1.4142$, … (rapid convergence).

Part 3. Systems of Linear Equations

We solve $A\mathbf{x} = \mathbf{b}$, where $A$ is an $n\times n$ nonsingular matrix.

Gauss Elimination Method

Definition 7. Gauss Elimination:
A direct method that reduces $A$ to upper triangular form by elementary row operations, then solves by back substitution.

Algorithm:

Step 1 (Forward elimination). For $k = 1,2,\ldots,n-1$:
Assume $a_{kk} \neq 0$ (pivot). For each row $i = k+1,\ldots,n$, set
$\ell_{ik} = \frac{a_{ik}}{a_{kk}},$
and replace row $i$ by
$\text{row}_i \leftarrow \text{row}_i - \ell_{ik}\cdot\text{row}_k$
(including the right-hand side $b_i$).

Step 2 (Back substitution). The system is now $U\mathbf{x} = \mathbf{c}$ with $U$ upper triangular.

Partial pivoting: If $|a_{kk}|$ is small, swap row $k$ with the row below it that has the largest $|a_{ik}|$ in column $k$. This improves numerical stability.

E.g., Solve

Eliminate $x$ from rows 2,3 using multipliers $\ell_{21} = -3/2$, $\ell_{31} = -1$; continue to triangular form; back-substitute.

Cholesky's Method

Definition 8. Cholesky Decomposition:
If $A$ is symmetric positive definite, then there exists a unique lower triangular matrix $L$ with positive diagonal entries such that
$A = L L^{T}.$
(Sometimes written $A = L D L^{T}$ with $D$ diagonal.)

Algorithm (compute $L = (\ell_{ij})$):
For $i = 1,\ldots,n$:

Then solve $L\mathbf{y} = \mathbf{b}$ (forward substitution) and $L^{T}\mathbf{x} = \mathbf{y}$ (back substitution).

Advantage: Roughly half the work of Gauss elimination when $A$ is SPD; numerically stable.

Gauss-Seidel Method

Definition 9. Iterative Methods:
Write $A = L_* + D + U_*$ (strictly lower + diagonal + strictly upper).
Gauss-Seidel updates unknowns in place, using the newest values as soon as they are available:

In matrix form:
$(D+L_*)\mathbf{x}^{(k+1)} = \mathbf{b} - U_*\mathbf{x}^{(k)}.$

Theorem 10 (Sufficient Condition)

If $A$ is strictly diagonally dominant, or symmetric positive definite, then the Gauss-Seidel iteration converges for any initial vector $\mathbf{x}^{(0)}$.

Remark: Gauss-Seidel usually converges faster than Jacobi (which uses only old values in a full sweep).

Method of Relaxation

Definition 10. Successive Over-Relaxation (SOR):
Accelerate Gauss-Seidel by extrapolating with a relaxation factor $\omega$:

where $\hat{x}_i^{(k+1)}$ is the Gauss-Seidel update for component $i$.

• $\omega = 1$: ordinary Gauss-Seidel.
• $0 < \omega < 1$: under-relaxation (may stabilize divergent processes).
• $1 < \omega < 2$: over-relaxation (often faster when it converges).

Choosing $\omega$: For many structured SPD systems there is an optimal $\omega \in (1,2)$; in exams one often takes a given $\omega$ (e.g. $1.2$) and iterates until consecutive approximations agree to the required tolerance.

Part 4. Ordinary Differential Equations

We consider the IVP

A numerical method produces values $y_n \approx y(x_n)$ at $x_n = x_0 + nh$.

Euler's Method

Theorem 11 (Euler's Method)

The Euler step is

Geometrically: advance along the tangent to the solution curve at $(x_n,y_n)$.
Local truncation error is $O(h^2)$; global error is $O(h)$ (first-order method).

Derivation:

By Taylor's theorem,

Dropping the $O(h^2)$ term and replacing $y(x_n)$ by $y_n$ gives Euler's formula.

E.g., $y' = x + y$, $y(0) = 1$, $h = 0.1$.
$y_1 = 1 + 0.1(0+1) = 1.1$,
$y_2 = 1.1 + 0.1(0.1+1.1) = 1.22$, and so on.

Single-Step Methods

Definition 11. Single-Step Method:
A method that computes $y_{n+1}$ using only the data at the current point $(x_n, y_n)$ (and evaluations of $f$ between $x_n$ and $x_{n+1}$).

Euler and classical Runge-Kutta are single-step.
Multistep methods (e.g. Milne-Simpson) use several previous points $y_n, y_{n-1}, \ldots$.

Modified Euler / Heun (improved Euler):

This is a second-order single-step method (average of slopes at the ends).

Runge-Kutta Method

Theorem 12 (Classical 4th-Order Runge-Kutta)

The RK4 formulae are

Local truncation error is $O(h^5)$; the method is of order $4$.

Idea: Weighted average of four slope estimates in $[x_n, x_n+h]$ — one at the left, two at the midpoint, one at the right — chosen so that the Taylor expansions match through order $h^4$.

Why popular in exams: One formula, no need for higher derivatives of $f$, good accuracy for moderate $h$.

Milne-Simpson Method

Definition 12. Predictor-Corrector (Milne-Simpson):
A multistep method using four previous points.

Milne's predictor (Adams-Bashforth type / Milne):

(Use the form prescribed in your textbook; both appear in Indian UG books.)

Simpson's corrector:

where $f_{n+1}^{(p)} = f\big(x_{n+1},\, y_{n+1}^{(p)}\big)$.
One may iterate the corrector until successive values agree.

Starting values: Milne-Simpson needs $y_0, y_1, y_2, y_3$. Compute the first few steps by a single-step method (e.g. RK4), then switch to Milne-Simpson.

Order: The corrector is related to Simpson's rule; the method is of order $4$ when used carefully.

Worked Numerical Patterns (Exam Style)

Pattern A — Lagrange (5–7 marks)

Q. Use Lagrange's formula to find $f(1.5)$ given $f(1)=1$, $f(2)=8$, $f(3)=27$ (for $f(x)=x^3$).

Solution.
Here $x_0=1$, $x_1=2$, $x_2=3$, $y_0=1$, $y_1=8$, $y_2=27$, and $x=1.5$.

(True value $1.5^3 = 3.375$; three points of a cubic need not recover the cubic exactly unless four points — here degree $\leq 2$ approx.)

Pattern B — Newton-Raphson (5–7 marks)

Q. Find a real root of $x^3 - 2x - 5 = 0$ correct to 3 decimal places by Newton-Raphson.

Solution.
Let $f(x) = x^3 - 2x - 5$, $f'(x) = 3x^2 - 2$.
$f(2) = -1 < 0$, $f(3) = 16 > 0$, so a root lies in $(2,3)$. Take $x_0 = 2$.

$x_1 = 2 - \dfrac{-1}{10} = 2.1$,
$x_2 = 2.1 - \dfrac{2.1^3 - 4.2 - 5}{3(2.1)^2-2} = 2.0945\ldots$,
Continue until $|x_{n+1}-x_n| < 0.0005$. Root $\approx 2.094$.

Pattern C — Euler (5 marks)

Q. Using Euler's method with $h = 0.2$, estimate $y(0.4)$ given $y' = y - x$, $y(0)=2$.

Solution.
$y_{n+1} = y_n + h(y_n - x_n)$.
$n=0$: $x_0=0$, $y_0=2$, $y_1 = 2 + 0.2(2-0) = 2.4$.
$n=1$: $x_1=0.2$, $y_2 = 2.4 + 0.2(2.4-0.2) = 2.84$.
Hence $y(0.4) \approx 2.84$.

Question Bank

1-Mark Questions

Definitions:

1. Define interpolation.
2. Write Lagrange's interpolating formula for $n=1$ (linear case).
3. Define first divided difference $f[x_0,x_1]$.
4. Define the forward difference operator $\Delta$.
5. Define the backward difference operator $\nabla$.
6. Write the relation $E = 1 + \Delta$.
7. State Newton's forward difference formula (up to second difference).
8. State the formula used in the bisection method.
9. Write the Regula-Falsi iteration formula.
10. Write the Newton-Raphson iteration formula.
11. What is partial pivoting in Gauss elimination?
12. When can Cholesky decomposition be applied?
13. Write one iteration formula of the Gauss-Seidel method for a $2\times 2$ system.
14. What is the relaxation factor $\omega$ in SOR?
15. Write Euler's formula for $y'=f(x,y)$.
16. Write the classical RK4 formula for $y_{n+1}$.
17. What is a predictor-corrector method?
18. Why does Milne-Simpson need starting values?

Short facts:

1. Order of convergence of Newton-Raphson (simple root).
2. Error bound after $n$ bisections on $[a,b]$.
3. If $f$ is a polynomial of degree $3$, what is $\Delta^4 y_i$?
4. Prefer Newton forward or backward when $x$ is near the end of the table?
5. Gauss-Seidel vs Jacobi — which uses updated values within a sweep?

5-Mark Questions

1. Derive Lagrange's interpolation formula for three points $(x_0,y_0)$, $(x_1,y_1)$, $(x_2,y_2)$.
2. Prove that divided differences are symmetric in their arguments.
3. Derive Newton's divided difference formula.
4. Prove $E = 1 + \Delta$ and $\Delta^n y_i = \nabla^n y_{i+n}$.
5. Derive Newton's forward difference formula from divided differences.
6. Describe the bisection method and prove $|x_n-\xi| \leq (b-a)/2^n$.
7. Derive the Regula-Falsi formula geometrically.
8. Derive the Newton-Raphson formula from the tangent line.
9. Explain Gauss elimination with a $3\times 3$ example (symbolic or numeric).
10. Write the algorithm for Cholesky decomposition.
11. Write Gauss-Seidel iteration for
    $a_{11}x + a_{12}y = b_1,\quad a_{21}x + a_{22}y = b_2.$
12. Explain the method of relaxation (SOR) and the role of $\omega$.
13. Derive Euler's method from Taylor's theorem.
14. Write RK4 formulae and explain the meaning of $k_1,k_2,k_3,k_4$.
15. Explain Milne's predictor and Simpson's corrector.

7 or 8-Mark Questions

1. State and derive Lagrange's interpolating formula. Hence find $f(2.5)$ from a given table of four values.
2. Define divided differences. Prove Newton's fundamental formula. Construct a divided difference table and interpolate.
3. Establish Newton's forward difference formula. Using a difference table, estimate $f(x)$ near the start of the table.
4. Compare bisection, Regula-Falsi, and Newton-Raphson (formula, geometry, rate of convergence, limitations).
5. Find a real root of a given transcendental equation by Newton-Raphson correct to 4 decimal places (full iteration table).
6. Solve a $3\times 3$ linear system by Gauss elimination with partial pivoting.
7. Solve $A\mathbf{x}=\mathbf{b}$ by Cholesky's method (given SPD $A$).
8. Perform 3 iterations of Gauss-Seidel for a diagonally dominant system; comment on convergence.
9. Using Euler's method, tabulate the solution of an IVP up to a given $x$; compare briefly with modified Euler.
10. Apply RK4 for one or two steps to a given IVP; show all $k_i$ calculations.
11. Describe the Milne-Simpson predictor-corrector method. Starting from RK4 values, compute one predictor-corrector step.