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

Gauss Elimination with Partial Pivoting

Gauss elimination converts the system $Ax=b$ into an upper-triangular form by row operations, then solves by back substitution.

Partial pivoting: before eliminating column $k$, search column $k$ (rows $k$ to $n$) for the entry of largest absolute value and swap that row into the pivot position. This avoids division by small (or zero) pivots, reducing round-off error and improving numerical stability.

Algorithm

1. Form the augmented matrix $[A\,|\,b]$.
2. For each column $k$: pick the row with max $|a_{ik}|$ ($i\ge k$), swap to pivot row.
3. Eliminate entries below the pivot: $R_i \leftarrow R_i - \frac{a_{ik}}{a_{kk}}R_k$.
4. Back-substitute from the last equation upward.

Solving the System

Augmented matrix:

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

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

Pivot column 2: $|\tfrac{10}{3}|>|{-\tfrac13}|$, swap $R_2\leftrightarrow R_3$:

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

Back substitution:

(Check: $2(7)-9+5=10$ ✓, $3(7)+2(-9)+3(5)=18$ ✓, $7+4(-9)+9(5)=16$ ✓.)

Derivation of Simpson's 1/3 Rule

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

Using the Newton forward formula with $x=x_0+ph$ and integrating from $p=0$ to $p=2$:

Extending over $n$ (even) strips gives the composite rule:

Evaluation

Apply the rule:

Exact value $=\ln 7=1.9459$; the small error arises from curvature in $f$.

Fourth-Order Runge-Kutta (RK4) Method

For $\dfrac{dy}{dx}=f(x,y)$, $y(x_0)=y_0$, RK4 advances one step of size $h$ using a weighted average of four slopes:

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

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

Step 1: $x_0=0,\;y_0=1$

So $y(0.1)\approx 1.11034$.

Step 2: $x_1=0.1,\;y_1=1.110342$

(Exact solution $y=2e^x-x-1$ gives $y(0.2)=1.24281$, confirming the result.)

Lagrange's Interpolation Formula

For $(n+1)$ data points $(x_0,y_0),\dots,(x_n,y_n)$ with unequally spaced $x$-values, the interpolating polynomial is

Each basis polynomial $L_i$ equals $1$ at $x_i$ and $0$ at all other nodes, so the formula passes exactly through every point. Its main advantage is it needs no difference table; its drawback is that adding a new point requires full recomputation.

Example

Given $(1,1),(2,8),(3,27)$, estimate $y$ at $x=2.5$:

So $y(2.5)\approx 16$ (true $2.5^3=15.625$).

Derivation of the Trapezoidal Rule

Approximate $\int_{x_0}^{x_1} f(x)\,dx$ over one strip of width $h=x_1-x_0$ by a straight line through $(x_0,y_0)$ and $(x_1,y_1)$. Using the Newton forward formula $f=y_0+p\,\Delta y_0$ with $x=x_0+ph$:

Composite rule over $n$ strips ($a=x_0,\;b=x_n,\;h=\tfrac{b-a}{n}$):

Error Term

For a single strip the truncation error is

For the composite rule over $[a,b]$:

The error is $O(h^2)$ and vanishes when $f$ is linear; it is exact for polynomials of degree $\le 1$.

Gauss-Seidel Iterative Method

An iterative technique to solve a diagonally dominant system $Ax=b$. Each equation is rearranged to express one unknown in terms of the others:

It differs from Jacobi by immediately using the newly computed values within the same iteration, so it generally converges faster.

Steps

1. Rearrange so the largest coefficients lie on the diagonal (ensures diagonal dominance $|a_{ii}|\ge\sum_{j\ne i}|a_{ij}|$, the sufficient condition for convergence).
2. Assume initial guesses (often all zero).
3. Compute $x_1,x_2,\dots,x_n$ in turn, using the latest available values.
4. Repeat until $|x_i^{(k+1)}-x_i^{(k)}|<\varepsilon$ for all $i$.

Example form for $10x+y+z=12,\;2x+10y+z=13,\;2x+2y+10z=14$:

which converges quickly to $x=y=z=1$. Gauss-Seidel saves memory and is well suited to large sparse systems.

Numerical Differentiation

When $f$ is known only at equally spaced points $x_i=x_0+ih$ (spacing $h$), derivatives are estimated from finite differences of the tabulated values $y_i$.

Forward Difference Formula

Used near the beginning of a table. With Newton's forward polynomial and $p=\tfrac{x-x_0}{h}$:

First-order approximation at $x_0$:

Backward Difference Formula

Used near the end of a table:

First-order approximation at $x_n$:

Second derivative

Forward differences use values ahead of the point; backward differences use values behind it. Both are $O(h)$ accurate; a smaller $h$ improves accuracy but excessive reduction amplifies round-off error.

Euler's Method

The simplest one-step method for $\dfrac{dy}{dx}=f(x,y),\;y(x_0)=y_0$. It uses the tangent line at the current point to estimate the next value:

Geometrically it follows the slope $f(x_n,y_n)$ over a step $h$. It is first-order accurate (global error $O(h)$), simple but requires small $h$ for acceptable accuracy.

Example

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

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

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

The exact value is $1.2428$; the gap shows Euler's first-order error.

Gauss Elimination vs Gauss-Jordan

Both are direct methods using elementary row operations on $[A\,|\,b]$, but they differ in how far the reduction is carried.

Summary: Gauss elimination is computationally cheaper but needs back substitution; Gauss-Jordan does more work to yield the solution (or inverse) directly without back substitution.

Newton's Divided Difference Interpolation

For $(n+1)$ points with unequally spaced $x$-values, the interpolating polynomial is

where the divided differences are defined recursively:

The leading (top) divided differences become the polynomial coefficients.

Advantages: works for unequal spacing, and a new data point can be added by computing just one more difference column without redoing the whole table (unlike Lagrange).

Example

For $(1,1),(2,8),(4,64)$:

At $x=3$: $f(3)=1+7(2)+7(2)(1)=1+14+14=29$.

Power Method (Largest Eigenvalue)

An iterative method to find the dominant (largest-magnitude) eigenvalue $\lambda_1$ and its eigenvector of a matrix $A$, provided $|\lambda_1|>|\lambda_2|\ge\cdots$.

Idea: repeatedly multiplying any starting vector by $A$ amplifies the component along the dominant eigenvector fastest.

Algorithm

1. Choose an initial nonzero vector $x^{(0)}$ (e.g. $[1,0,0]^T$).
2. Compute $y^{(k+1)}=A\,x^{(k)}$.
3. The largest-magnitude component of $y^{(k+1)}$ is the current eigenvalue estimate $\lambda^{(k+1)}$.
4. Normalize: $x^{(k+1)}=\dfrac{y^{(k+1)}}{\lambda^{(k+1)}}$.
5. Repeat until $\lambda$ and $x$ stabilize within tolerance.

Then $\lambda^{(k)}\to\lambda_1$ and $x^{(k)}\to$ the corresponding eigenvector.

Example (one iteration): $A=\begin{bmatrix}2&1\\1&2\end{bmatrix}$, $x^{(0)}=\begin{bmatrix}1\\0\end{bmatrix}$.

Continuing converges to $\lambda_1=3$, eigenvector $[1,1]^T$.

Note: convergence rate depends on the ratio $|\lambda_2/\lambda_1|$; the smallest eigenvalue is found by applying the method to $A^{-1}$ (inverse power method).

Least-Squares Fit of a Parabola $y=a+bx+cx^2$

Given $n$ data points $(x_i,y_i)$, choose $a,b,c$ to minimize the sum of squared residuals

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

Procedure

1. Tabulate $x,\;y,\;x^2,\;x^3,\;x^4,\;xy,\;x^2y$ for all points and form the required sums.
2. Substitute the sums into the three normal equations.
3. Solve the resulting $3\times3$ linear system (e.g. by Gauss elimination) for $a,b,c$.
4. Write the fitted parabola $y=a+bx+cx^2$.

Example: for $(0,1),(1,6),(2,17),(3,34)$ ($n=4$), $\sum x=6,\sum x^2=14,\sum x^3=36,\sum x^4=98,\sum y=58,\sum xy=128,\sum x^2 y=346$. Solving the normal equations yields $a=1,\,b=2,\,c=3$, i.e. $y=1+2x+3x^2$, which fits the data exactly.