BSc CSIT (TU) Science Numerical Method (BSc CSIT, CSC207) Question Paper 2075 Nepal

Gauss Elimination with Partial Pivoting

Idea. Convert the augmented matrix $[A\,|\,b]$ into an upper-triangular form by elementary row operations, then solve by back substitution.

Partial pivoting. Before eliminating column $k$, search rows $k, k+1, \dots, n$ for the entry with the largest absolute value in that column and swap it into the pivot (diagonal) position. This avoids division by a small/zero pivot and controls round-off error, improving numerical stability.

Algorithm.

1. Form augmented matrix $[A\,|\,b]$.
2. For each column $k$: pick the row with $\max|a_{ik}|$ and swap to row $k$ (partial pivot).
3. Eliminate entries below the pivot using $R_i \to R_i - \frac{a_{ik}}{a_{kk}}R_k$.
4. After triangularization, back-substitute: $x_n = b_n/a_{nn}$, then $x_i = \frac{1}{a_{ii}}\left(b_i - \sum_{j>i} a_{ij}x_j\right)$.

Solving the System

Augmented matrix:

Step 1 — Pivot column 1. Largest $|a_{i1}|$ is $3$ (row 2). Swap $R_1 \leftrightarrow R_2$:

Eliminate: $R_2 \to R_2 - \tfrac{2}{3}R_1$, $R_3 \to R_3 - \tfrac{1}{3}R_1$:

Step 2 — Pivot column 2. $|\tfrac{10}{3}| > |\tfrac{1}{3}|$, so swap $R_2 \leftrightarrow R_3$:

Eliminate: $R_3 \to R_3 - \frac{-1/3}{10/3}R_2 = R_3 + \tfrac{1}{10}R_2$:

Back substitution.

Solution: $\boxed{x = 7,\; y = -9,\; z = 5}$ (verified in all three original equations).

Derivation of Simpson's 1/3 Rule

Approximate $\int_{x_0}^{x_2} f(x)\,dx$ over two equal subintervals of width $h$ (so three points $x_0, x_1, x_2$) by fitting a second-degree (parabolic) polynomial through $(x_0,y_0),(x_1,y_1),(x_2,y_2)$.

Using Newton's forward interpolation with $x = x_0 + sh$:

Integrating from $x_0$ to $x_2$ (i.e. $s = 0$ to $2$), with $dx = h\,ds$:

Substituting $\Delta y_0 = y_1 - y_0$ and $\Delta^2 y_0 = y_2 - 2y_1 + y_0$:

Composite rule (for $n$ even subintervals): the pattern is $1,4,2,4,\dots,4,1$:

Evaluation of $\int_0^6 \frac{dx}{1+x}$, $n = 6$

Here $a=0,\ b=6,\ n=6 \Rightarrow h = \frac{6-0}{6} = 1$, with $f(x)=\frac{1}{1+x}$.

Apply the rule:

Check. Exact value $= \ln 7 = 1.94591$, so the approximation is accurate to about $0.013$.

Fourth-Order Runge-Kutta (RK4) Method

For the IVP $\dfrac{dy}{dx} = f(x,y),\ y(x_0)=y_0$, the RK4 method advances the solution by one step $h$ using a weighted average of four slope estimates:

The local truncation error is $O(h^5)$ and the global error $O(h^4)$, giving high accuracy without computing derivatives of $f$.

Solve $\dfrac{dy}{dx} = x + y,\ y(0)=1$, find $y(0.2)$, $h=0.1$

Here $f(x,y) = x + y$.

Step 1: $x_0 = 0,\ y_0 = 1 \to y(0.1)$

Step 2: $x_1 = 0.1,\ y_1 = 1.110342 \to y(0.2)$

Result: $y(0.2) \approx 1.2428$. (Exact $y = 2e^x - x - 1$ gives $y(0.2)=1.242806$ — excellent agreement.)

Lagrange's Interpolation Formula

Given $n+1$ data points $(x_0,y_0),(x_1,y_1),\dots,(x_n,y_n)$ with distinct, not necessarily equally spaced $x_i$, the unique interpolating polynomial of degree $\le n$ is:

Each Lagrange basis polynomial $L_i(x)$ has the property $L_i(x_k) = 1$ if $k=i$ and $0$ otherwise, so $P(x_k)=y_k$ at every node. Its main advantage is that it needs no divided-difference table and works for unequal spacing; its drawback is that adding a new point requires recomputing all basis polynomials.

Example

Given $(1,1),(2,4),(3,9)$, estimate $f(2.5)$:

At $x=2.5$: $L_0 = \frac{(0.5)(-0.5)}{2} = -0.125$, $L_1 = \frac{(1.5)(-0.5)}{-1} = 0.75$, $L_2 = \frac{(1.5)(0.5)}{2} = 0.375$.

This matches $f(x)=x^2 \Rightarrow (2.5)^2 = 6.25$.

Derivation of the Trapezoidal Rule

To approximate $\int_{x_0}^{x_1} f(x)\,dx$ over one subinterval of width $h = x_1 - x_0$, fit a straight line (first-degree polynomial) through $(x_0,y_0)$ and $(x_1,y_1)$. Using Newton's forward formula with $x = x_0 + sh$:

Integrating over $s\in[0,1]$, $dx = h\,ds$:

Composite trapezoidal rule over $n$ subintervals ($h=(b-a)/n$):

Error Term

For a single strip, the error is

For the composite rule the total error is

Thus the trapezoidal rule is exact for linear functions and the error is proportional to the second derivative, so it is zero when $f''=0$.

Gauss-Seidel Iterative Method

The Gauss-Seidel method solves $Ax = b$ iteratively. For a system

each equation is solved for its diagonal variable:

Key feature: unlike Jacobi's method, it uses the most recently updated values immediately within the same iteration, so it generally converges faster.

Convergence condition: the method is guaranteed to converge if the coefficient matrix is diagonally dominant, i.e. $|a_{ii}| \ge \sum_{j\ne i}|a_{ij}|$ for each row (strict for at least one). Rearrange equations to make $A$ diagonally dominant before iterating.

Procedure.

1. Rearrange for diagonal dominance; assume initial guess (e.g. all zeros).
2. Compute $x_1, x_2, \dots, x_n$ in order, immediately reusing new values.
3. Repeat until $\max_i |x_i^{(k+1)} - x_i^{(k)}| < \varepsilon$ (tolerance).

Example. For $10x + y + z = 12,\ 2x + 10y + z = 13,\ 2x + 2y + 10z = 14$ (diagonally dominant), iterating from $(0,0,0)$ converges quickly to $x=y=z=1$.

Numerical Differentiation by Finite Differences

Given tabulated values $y_i = f(x_i)$ at equally spaced points with step $h$, the derivative of $f$ is approximated using a difference table built from an interpolating polynomial.

Forward Difference Formula (near the start of the table)

Using Newton's forward interpolation $f(x) = y_0 + s\Delta y_0 + \frac{s(s-1)}{2!}\Delta^2 y_0 + \cdots$ with $s=\frac{x-x_0}{h}$ and differentiating w.r.t. $x$ ($\frac{ds}{dx}=\frac1h$):

Simplest two-point form: $f'(x_0) \approx \dfrac{y_1 - y_0}{h}$, error $O(h)$.

Backward Difference Formula (near the end of the table)

Using Newton's backward interpolation and differentiating at $x_n$:

Simplest form: $f'(x_n) \approx \dfrac{y_n - y_{n-1}}{h}$, error $O(h)$.

Usage. Use the forward formula to estimate derivatives near the beginning of the table and the backward formula near the end; both are first-order accurate (error $\propto h$), while a central difference $\frac{y_{i+1}-y_{i-1}}{2h}$ is second-order accurate. Second derivatives follow by differentiating the series twice.

Euler's Method

Euler's method is the simplest single-step method for the IVP $\dfrac{dy}{dx}=f(x,y),\ y(x_0)=y_0$. It approximates the curve by its tangent line over each step $h$:

Geometrically, starting at $(x_n,y_n)$ we move along the slope $f(x_n,y_n)$ for a horizontal distance $h$. It is a first-order method with local truncation error $O(h^2)$ and global error $O(h)$, so small $h$ is needed for accuracy.

Example

Solve $\dfrac{dy}{dx} = x + y,\ y(0)=1$ for $y(0.2)$ with $h=0.1$ (here $f(x,y)=x+y$):

Step 1 ($x_0=0, y_0=1$):

Step 2 ($x_1=0.1, y_1=1.1$):

So $y(0.2) \approx 1.22$. (Exact value $1.2428$; the gap reflects Euler's first-order error, reducible with smaller $h$.)

Gauss Elimination vs. Gauss-Jordan

Both are direct methods for solving $Ax=b$ via elementary row operations on the augmented matrix, but they differ in how far the reduction is carried.

Summary: Gauss elimination triangularizes then back-substitutes; Gauss-Jordan fully diagonalizes so the solution appears directly. Gauss-Jordan does more work but is the standard route for computing the inverse of a matrix.

Newton's Divided Difference Interpolation

For $n+1$ points $(x_0,y_0),\dots,(x_n,y_n)$ with distinct, possibly unequally spaced $x_i$, the interpolating polynomial is

Divided differences are defined recursively:

Properties / advantages.

• Works for unequally spaced data (unlike Newton's forward/backward formulae).
• Symmetric in arguments — the result is independent of the order of points.
• A new data point can be added by appending one extra term without recomputing previous ones.

Example. Given $(1,1),(2,8),(4,64)$: $f[x_0,x_1]=\frac{8-1}{2-1}=7$, $f[x_1,x_2]=\frac{64-8}{4-2}=28$, $f[x_0,x_1,x_2]=\frac{28-7}{4-1}=7$. Hence $P(x)=1 + 7(x-1) + 7(x-1)(x-2)$, which evaluates the function at any required $x$.

Power Method (Largest Eigenvalue)

The power method is an iterative scheme to find the dominant eigenvalue (largest in magnitude) and its eigenvector of a square matrix $A$, provided one eigenvalue strictly dominates: $|\lambda_1| > |\lambda_2| \ge \cdots$.

Algorithm.

1. Choose a non-zero initial vector $x^{(0)}$ (e.g. $[1,0,0]^T$ or all ones).
2. Iterate: $y^{(k+1)} = A\,x^{(k)}$.
3. The dominant eigenvalue estimate is the largest-magnitude component of $y^{(k+1)}$:

4. Normalize to control growth: $x^{(k+1)} = \dfrac{y^{(k+1)}}{\lambda^{(k+1)}}$.
5. Repeat until $\lambda^{(k+1)}$ and $x^{(k+1)}$ stabilize within tolerance.

At convergence, $\lambda^{(k)} \to \lambda_1$ (largest eigenvalue) and $x^{(k)} \to$ corresponding eigenvector.

Why it works. Expressing $x^{(0)} = \sum c_i v_i$ in the eigenbasis, $A^k x^{(0)} = \sum c_i \lambda_i^k v_i = \lambda_1^k\big(c_1 v_1 + \sum_{i\ge2} c_i (\lambda_i/\lambda_1)^k v_i\big)$. Since $|\lambda_i/\lambda_1|<1$, all terms except the dominant one vanish, leaving the direction of $v_1$.

Note. Convergence speed depends on the ratio $|\lambda_2/\lambda_1|$; the inverse power method ($A^{-1}$) gives the smallest eigenvalue.

Fitting a Second-Degree Parabola by Least Squares

We fit $y = a + bx + cx^2$ to $n$ data points $(x_i,y_i)$ by minimizing the sum of squared residuals

Setting $\dfrac{\partial S}{\partial a} = \dfrac{\partial S}{\partial b} = \dfrac{\partial S}{\partial c} = 0$ gives the three normal equations:

Procedure.

1. Build a table computing $x, x^2, x^3, x^4, y, xy, x^2y$ for each point.
2. Form the column sums $\sum x,\ \sum x^2,\ \dots,\ \sum x^2 y$.
3. Substitute into the three normal equations.
4. Solve the resulting $3\times3$ linear system (e.g. by Gauss elimination) for $a, b, c$.

Example. For $(1,1),(2,5),(3,10),(4,22)$, $n=4$: $\sum x=10,\ \sum x^2=30,\ \sum x^3=100,\ \sum x^4=354,\ \sum y=38,\ \sum xy=125,\ \sum x^2y=437$. Substituting and solving the three normal equations yields the best-fit coefficients $a, b, c$, giving the required parabola $y = a + bx + cx^2$.