BE Computer Engineering (IOE, TU) Numerical Methods (IOE, SH 553) Question Paper 2078 Nepal

(a) Newton-Raphson method

Derivation. Let $x_n$ be an approximation to the root of $f(x)=0$. Expand $f$ in a Taylor series about $x_n$ and retain only the linear term:

$$f(x_{n+1}) \approx f(x_n) + (x_{n+1}-x_n)\,f'(x_n) = 0.$$

Solving for $x_{n+1}$ gives the Newton-Raphson iterative formula:

$$\boxed{\,x_{n+1} = x_n - \dfrac{f(x_n)}{f'(x_n)}\,}, \qquad f'(x_n)\neq 0.$$

Geometrical interpretation. At the point $(x_n, f(x_n))$ draw the tangent to the curve $y=f(x)$. The tangent has slope $f'(x_n)$; the point where this tangent crosses the $x$-axis is taken as the next approximation $x_{n+1}$. Thus each iteration replaces the curve locally by its tangent line and uses the tangent's $x$-intercept as the improved estimate. Repeating this drives the iterate towards the actual root, and the convergence is quadratic near a simple root.

Two situations of failure:

1. $f'(x_n)$ is zero or very small (a horizontal/near-horizontal tangent): the tangent is nearly parallel to the $x$-axis, so $x_{n+1}$ shoots far away and the iteration diverges or overflows.
2. Poor initial guess / unsuitable function shape (e.g. a point near a local maximum or minimum, or oscillation between two values): the successive iterates may cycle or move away from the root instead of converging.

(b) Root of $x^3-5x+1=0$ in $(0,1)$

Here $f(x)=x^3-5x+1$ and $f'(x)=3x^2-5$. Take $x_0=0.5$.

Since $x_3$ and $x_2$ agree to four decimals,

$$\boxed{x \approx 0.2016}.$$

(a) Natural cubic spline

Given $(n+1)$ points $(x_0,y_0),\dots,(x_n,y_n)$ with $x_0<x_1<\cdots<x_n$, a cubic spline $S(x)$ is a piecewise cubic: on each sub-interval $[x_i,x_{i+1}]$ it is a separate cubic $S_i(x)=a_i+b_i(x-x_i)+c_i(x-x_i)^2+d_i(x-x_i)^3$ ($i=0,\dots,n-1$), giving $4n$ unknown coefficients.

Continuity conditions at the interior knots $x_1,\dots,x_{n-1}$:

1. Interpolation: $S_i(x_i)=y_i$ and $S_i(x_{i+1})=y_{i+1}$ — the spline passes through every data point and is continuous.
2. First-derivative continuity: $S_{i-1}'(x_i)=S_i'(x_i)$ — the slopes match at each interior knot.
3. Second-derivative continuity: $S_{i-1}''(x_i)=S_i''(x_i)$ — the curvatures match, giving a smooth ($C^2$) curve.

These conditions give $4n-2$ equations. Writing $M_i=S''(x_i)$, the second-derivative continuity reduces to the tridiagonal system

$$h_{i-1}M_{i-1}+2(h_{i-1}+h_i)M_i+h_iM_{i+1}=6\left(\frac{y_{i+1}-y_i}{h_i}-\frac{y_i-y_{i-1}}{h_{i-1}}\right),\quad h_i=x_{i+1}-x_i.$$

Boundary conditions (natural spline): the curvature is set to zero at the two ends,

$$S''(x_0)=M_0=0, \qquad S''(x_n)=M_n=0.$$

These two extra conditions complete the system, which is solved (tridiagonal) for the $M_i$ and hence the cubic coefficients.

(b) Newton's divided difference, estimate $f(3)$

Data: $x=0,1,2,4,5$, $f=1,3,9,33,51$.

Divided-difference table (leading coefficients):

Leading coefficients: $a_0=1,\ a_1=2,\ a_2=2,\ a_3=0,\ a_4=0$.

$$P(x)=1+2(x-0)+2(x-0)(x-1)+0+0.$$

At $x=3$:

$$P(3)=1+2(3)+2(3)(2)=1+6+12=\boxed{19}.$$

Thus $f(3)\approx 19$. (Indeed $P(x)=2x^2+1$ fits the data exactly.)

(a) Gauss elimination with partial pivoting

For $[A]\{x\}=\{b\}$ the method reduces the augmented matrix to upper-triangular form by elementary row operations, then back-substitutes.

Procedure (for column $k=1,2,\dots,n-1$):

1. Pivot search: among the entries $a_{kk},a_{k+1,k},\dots,a_{nk}$ in the current column, find the one of largest absolute value.
2. Row interchange: swap that row with row $k$, so the largest element becomes the pivot $a_{kk}$.
3. Elimination: for each row $i>k$, compute the multiplier $m_{ik}=a_{ik}/a_{kk}$ and replace row $i$ by $\text{row}_i - m_{ik}\,\text{row}_k$, making all entries below the pivot zero.
4. After triangularisation, solve by back substitution $x_n=b_n/a_{nn}$, $x_i=\big(b_i-\sum_{j>i}a_{ij}x_j\big)/a_{ii}$.

Why partial pivoting is necessary / improves stability. Without pivoting the pivot $a_{kk}$ may be zero (division by zero) or very small. A small pivot makes the multipliers $m_{ik}$ very large, so rounding errors already present in the entries are amplified during elimination, badly corrupting the answer. By always choosing the largest-magnitude available pivot, all multipliers satisfy $|m_{ik}|\le 1$, which keeps the growth of round-off error bounded and greatly improves numerical stability.

(b) Gauss-Seidel, three iterations

The system is diagonally dominant. Rearranged iteration equations:

$$x=\frac{17-y+2z}{20},\quad y=\frac{-18-3x+z}{20},\quad z=\frac{25-2x+3y}{20}.$$

Start $(x,y,z)=(0,0,0)$, using newest values immediately:

After three iterations,

$$\boxed{x\approx 1.000,\quad y\approx -1.000,\quad z\approx 1.000}.$$

(a) Five-point formula for the Laplace equation

Discretise the region with a uniform mesh of spacing $h$ in both $x$ and $y$. Denote $u_{i,j}=u(x_i,y_j)$. Using central second differences:

$$\frac{\partial^2 u}{\partial x^2}\approx\frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{h^2},\qquad \frac{\partial^2 u}{\partial y^2}\approx\frac{u_{i,j+1}-2u_{i,j}+u_{i,j-1}}{h^2}.$$

Substituting into $u_{xx}+u_{yy}=0$ and multiplying by $h^2$:

$$u_{i+1,j}+u_{i-1,j}+u_{i,j+1}+u_{i,j-1}-4u_{i,j}=0,$$

so the standard five-point (Laplacian) formula is

$$\boxed{\,u_{i,j}=\tfrac14\big(u_{i+1,j}+u_{i-1,j}+u_{i,j+1}+u_{i,j-1}\big)\,}.$$

Each interior value equals the average of its four neighbours (left, right, top, bottom).

(b) Square plate, four interior nodes

Label the interior nodes on a $2\times2$ interior grid: top row $u_1$ (left), $u_2$ (right); bottom row $u_3$ (left), $u_4$ (right). Boundary: top edge $=100$, the other three edges $=0$.

Applying the five-point average to each node:

$$4u_1=100+u_2+u_3,\quad 4u_2=100+u_1+u_4,$$ $$4u_3=u_1+u_4+0+0,\quad 4u_4=u_2+u_3+0+0.$$

Gauss-Seidel (Liebmann) iteration, starting from $0$ and using newest values:

By symmetry $u_1=u_2$ and $u_3=u_4$. The converged interior values are

$$\boxed{u_1=u_2=37.5,\qquad u_3=u_4=12.5}.$$

Bisection method for $f(x)=xe^{x}-2=0$ on $[0,1]$

Check signs: $f(0)=-2<0$, $f(1)=e-2=0.7183>0$, so a root lies in $[0,1]$. Each step takes the midpoint $c=(a+b)/2$ and keeps the half-interval where the sign change occurs.

The iterates have settled to three decimal places:

$$\boxed{x\approx 0.853}.$$

Simpson's 1/3 rule for $\displaystyle\int_0^1\frac{dx}{1+x^2}$, $h=0.25$

Nodes ($n=4$ intervals) and ordinates $f(x)=\dfrac{1}{1+x^2}$:

Simpson's 1/3 rule:

$$I=\frac{h}{3}\Big[f_0+f_4+4(f_1+f_3)+2f_2\Big]$$ $$=\frac{0.25}{3}\big[1.0+0.5+4(0.941176+0.640000)+2(0.800000)\big]$$ $$=\frac{0.25}{3}\,[1.5+6.324706+1.6]=\frac{0.25}{3}(9.424706)=\boxed{0.785392}.$$

Comparison with exact value. The exact value is $\int_0^1\frac{dx}{1+x^2}=\tan^{-1}(1)=\dfrac{\pi}{4}=0.785398$. The error is only $\approx 6\times10^{-6}$, showing Simpson's rule is very accurate here.

Improving accuracy with Romberg integration. Romberg integration starts from the composite trapezoidal rule $R_{k,1}$ for successively halved step sizes and then applies Richardson extrapolation

$$R_{k,j}=\frac{4^{\,j-1}R_{k,j-1}-R_{k-1,j-1}}{4^{\,j-1}-1}$$

to systematically cancel the leading error terms. The first extrapolation column already reproduces Simpson's rule, and each further column removes higher-order error, so accuracy improves rapidly (effectively raising the order of the method) without computing many extra function values.

Least-squares straight line $y=a+bx$

With $n=5$ points, form the required sums:

Normal equations:

$$\sum y = na + b\sum x,\qquad \sum xy = a\sum x + b\sum x^2$$ $$30.1 = 5a+15b,\qquad 110.4 = 15a+55b.$$

Solve:

$$b=\frac{n\sum xy-\sum x\sum y}{n\sum x^2-(\sum x)^2}=\frac{5(110.4)-15(30.1)}{5(55)-15^2}=\frac{552-451.5}{275-225}=\frac{100.5}{50}=2.01.$$ $$a=\frac{\sum y-b\sum x}{n}=\frac{30.1-2.01(15)}{5}=\frac{30.1-30.15}{5}=-0.01.$$

Fitted line:

$$\boxed{y=-0.01+2.01\,x}.$$

Estimate at $x=6$:

$$y(6)=-0.01+2.01(6)=12.05.$$

Central difference formulae from Taylor's series

Expand $f$ about $x$ with step $h$:

$$f(x+h)=f(x)+hf'(x)+\tfrac{h^2}{2}f''(x)+\tfrac{h^3}{6}f'''(x)+\cdots$$ $$f(x-h)=f(x)-hf'(x)+\tfrac{h^2}{2}f''(x)-\tfrac{h^3}{6}f'''(x)+\cdots$$

First derivative: subtract the second from the first; the even terms cancel:

$$f(x+h)-f(x-h)=2hf'(x)+\tfrac{h^3}{3}f'''(x)+\cdots$$ $$\Rightarrow\ f'(x)=\frac{f(x+h)-f(x-h)}{2h}+O(h^2).$$

Second derivative: add the two expansions; the odd terms cancel:

$$f(x+h)+f(x-h)=2f(x)+h^2f''(x)+\tfrac{h^4}{12}f^{(4)}(x)+\cdots$$ $$\Rightarrow\ f''(x)=\frac{f(x+h)-2f(x)+f(x-h)}{h^2}+O(h^2).$$

Both central-difference formulae are second-order accurate, i.e. the truncation error is $O(h^2)$.

Numerical evaluation at $x=2.0$, $h=0.2$

Data: $f(1.8)=6.050,\ f(2.0)=7.389,\ f(2.2)=9.025$.

$$f'(2.0)=\frac{f(2.2)-f(1.8)}{2h}=\frac{9.025-6.050}{0.4}=\frac{2.975}{0.4}=7.4375.$$ $$f''(2.0)=\frac{f(2.2)-2f(2.0)+f(1.8)}{h^2}=\frac{9.025-2(7.389)+6.050}{0.04}=\frac{0.297}{0.04}=7.425.$$ $$\boxed{f'(2.0)\approx 7.4375,\qquad f''(2.0)\approx 7.425.}$$

(Consistent with $f=e^x$, for which $f'(2)=f''(2)=e^2\approx 7.389$.)

Fourth-order Runge-Kutta for $y'=x+y$, $y(0)=1$, $h=0.1$

The RK4 update is

$$k_1=hf(x,y),\ k_2=hf\!\big(x+\tfrac h2,y+\tfrac{k_1}2\big),\ k_3=hf\!\big(x+\tfrac h2,y+\tfrac{k_2}2\big),\ k_4=hf(x+h,y+k_3),$$ $$y_{n+1}=y_n+\tfrac16(k_1+2k_2+2k_3+k_4),\qquad f(x,y)=x+y.$$

Step 1: $x_0=0,\ y_0=1\ \to\ x=0.1$

$$k_1=0.1(0+1)=0.100000$$ $$k_2=0.1(0.05+1.05)=0.110000$$ $$k_3=0.1(0.05+1.055)=0.110500$$ $$k_4=0.1(0.1+1.1105)=0.121050$$ $$y_1=1+\tfrac16(0.1+0.22+0.221+0.12105)=1.110342.$$

Step 2: $x_1=0.1,\ y_1=1.110342\ \to\ x=0.2$

$$k_1=0.1(0.1+1.110342)=0.121034$$ $$k_2=0.1(0.15+1.170859)=0.132086$$ $$k_3=0.1(0.15+1.176385)=0.132638$$ $$k_4=0.1(0.2+1.242980)=0.144298$$ $$y_2=1.110342+\tfrac16(0.121034+0.264172+0.265277+0.144298)=1.242805.$$ $$\boxed{y(0.2)\approx 1.2428.}$$

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

Lagrange interpolation for $f(2)$

Data: $x_0=0,x_1=1,x_2=3,x_3=4$ with $f=-12,0,12,24$.

Lagrange formula:

$$f(x)=\sum_{i=0}^{3} y_i\,L_i(x),\qquad L_i(x)=\prod_{j\neq i}\frac{x-x_j}{x_i-x_j}.$$

Evaluate each basis term at $x=2$:

$$L_0(2)=\frac{(2-1)(2-3)(2-4)}{(0-1)(0-3)(0-4)}=\frac{(1)(-1)(-2)}{(-1)(-3)(-4)}=\frac{2}{-12}=-\tfrac16,$$ $$L_1(2)=\frac{(2-0)(2-3)(2-4)}{(1-0)(1-3)(1-4)}=\frac{(2)(-1)(-2)}{(1)(-2)(-3)}=\frac{4}{6}=\tfrac23,$$ $$L_2(2)=\frac{(2-0)(2-1)(2-4)}{(3-0)(3-1)(3-4)}=\frac{(2)(1)(-2)}{(3)(2)(-1)}=\frac{-4}{-6}=\tfrac23,$$ $$L_3(2)=\frac{(2-0)(2-1)(2-3)}{(4-0)(4-1)(4-3)}=\frac{(2)(1)(-1)}{(4)(3)(1)}=\frac{-2}{12}=-\tfrac16.$$

Therefore

$$f(2)=(-12)(-\tfrac16)+(0)(\tfrac23)+(12)(\tfrac23)+(24)(-\tfrac16)=2+0+8-4=\boxed{6}.$$

($f(2)\approx 6$; the data in fact lie on the line $f(x)=6x-12$.)

Shooting method for a two-point BVP

Consider a second-order BVP

$$y''=f(x,y,y'),\qquad y(a)=\alpha,\quad y(b)=\beta.$$

The difficulty is that only one condition ($y(a)=\alpha$) is known at the starting end; the slope $y'(a)$ there is unknown.

Converting BVP to IVP. Introduce the unknown initial slope as a parameter $s$, i.e. assume

$$y'(a)=s.$$

Now the problem becomes an initial value problem with both $y(a)=\alpha$ and $y'(a)=s$ known, which can be integrated forward from $a$ to $b$ using any IVP solver (RK4, Euler, etc.). The computed end value $y(b)$ depends on the guessed slope, call it $\phi(s)=y(b;s)$.

Hitting the target (the "shooting"). We need $\phi(s)=\beta$. We:

1. Choose two trial slopes $s_1,s_2$, integrate the IVP, and obtain $\phi(s_1),\phi(s_2)$.
2. Since $\phi(s)-\beta=0$ is a root-finding problem in $s$, improve the guess using linear interpolation (secant method):

$$s_3=s_2-\big(\phi(s_2)-\beta\big)\frac{s_2-s_1}{\phi(s_2)-\phi(s_1)}.$$

3. Repeat until $y(b)$ matches the prescribed $\beta$ to the desired tolerance.

Thus the method "shoots" trajectories from $x=a$ with different initial slopes and adjusts the slope (like aiming a gun) until the trajectory lands on the required boundary value at $x=b$ — converting a boundary value problem into a sequence of initial value problems combined with a root-finding step.

LU (Doolittle) decomposition

The idea is to factor the coefficient matrix as a product of a lower- and an upper-triangular matrix:

$$[A]=[L][U],$$

where in the Doolittle form $[L]$ is unit lower-triangular (1's on the diagonal) and $[U]$ is upper-triangular:

$$L=\begin{bmatrix}1&0&0\\ l_{21}&1&0\\ l_{31}&l_{32}&1\end{bmatrix},\qquad U=\begin{bmatrix}u_{11}&u_{12}&u_{13}\\ 0&u_{22}&u_{23}\\ 0&0&u_{33}\end{bmatrix}.$$

The entries are found by equating $[A]=[L][U]$ element by element (row of $U$, then column of $L$, alternately).

Solving $[A]\{x\}=\{b\}$. Write $[L][U]\{x\}=\{b\}$ and set $[U]\{x\}=\{y\}$. Then:

1. Forward substitution solve $[L]\{y\}=\{b\}$ for $\{y\}$.
2. Back substitution solve $[U]\{x\}=\{y\}$ for $\{x\}$.

Main advantage over Gauss elimination (multiple right-hand sides). The expensive part — factoring $[A]$ into $[L][U]$ — depends only on $[A]$ and is done once ($\sim n^3/3$ operations). For each new right-hand-side vector $\{b\}$, only the two cheap triangular solves (forward + back substitution, $\sim n^2$ operations each) are repeated. In contrast, plain Gauss elimination would redo the full $O(n^3)$ elimination for every new $\{b\}$. Hence LU decomposition is far more efficient when $[A]\{x\}=\{b\}$ must be solved for several different $\{b\}$ (and it also gives $\det A=\prod u_{ii}$ and $A^{-1}$ readily).