BE Computer Engineering (Pokhara University) Numerical Methods (PU, MTH 252) Question Paper 2078 Nepal

(a) Newton-Raphson Formula, Geometry and Failure

Derivation. Let $x_n$ be an approximation to a root of $f(x)=0$. Expand $f$ in a Taylor series about $x_n$ and retain terms up to the first order:

$$f(x_n+h)\approx f(x_n)+h\,f'(x_n)=0\quad\Rightarrow\quad h=-\frac{f(x_n)}{f'(x_n)}.$$

Thus the next approximation $x_{n+1}=x_n+h$ gives the iterative formula

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

Geometric interpretation. At the point $(x_n,f(x_n))$ draw the tangent to the curve $y=f(x)$. This tangent meets the $x$-axis at $x_{n+1}$. The method therefore successively replaces the curve by its tangent, and the point where the tangent cuts the $x$-axis is taken as the improved root.

Conditions for failure / non-convergence. The method may fail or diverge when:

1. $f'(x_n)=0$ or is very small (tangent is horizontal or nearly so), giving a huge correction.
2. The initial guess $x_0$ is far from the root or lies near a turning point/inflection point, so iterates may oscillate or diverge.
3. There is a multiple root, where convergence degrades from quadratic to linear.
4. The function is not sufficiently smooth ($f'$ not continuous) near the root.

When it does converge near a simple root the convergence is quadratic, i.e. the error satisfies $e_{n+1}\approx C\,e_n^{2}$.

(b) Root of $x^3-2x-5=0$

Here $f(x)=x^{3}-2x-5,\ f'(x)=3x^{2}-2$, with $x_0=2$.

After three iterations the root is

$$\boxed{x\approx 2.0946\ \text{(correct to 4 decimal places).}}$$

Comparison with bisection. The bisection method halves the interval each step, so its error decreases linearly: $e_{n+1}=\tfrac12 e_n$, needing about $\log_2(\text{range}/\text{tol})\approx 14$ steps for $10^{-4}$ accuracy on $[2,3]$. Newton-Raphson roughly doubles the number of correct digits each iteration (quadratic convergence) and reached 4-decimal accuracy in only 3 iterations. Hence Newton-Raphson converges much faster, though bisection is always convergent once the root is bracketed.

(a) Newton's Forward Difference Interpolation Formula

For equally spaced nodes $x_i=x_0+ih$, with forward differences $\Delta y_0,\Delta^2 y_0,\dots$, put $p=\dfrac{x-x_0}{h}$. The interpolating polynomial is

$$y(x)=y_0+p\,\Delta y_0+\frac{p(p-1)}{2!}\Delta^2 y_0+\frac{p(p-1)(p-2)}{3!}\Delta^3 y_0+\cdots$$

Outline of derivation. Using the shift operator $E$ with $E\,y_0=y_1$ and $E=1+\Delta$, we have $y_p=E^{p}y_0=(1+\Delta)^{p}y_0$. Expanding $(1+\Delta)^p$ by the binomial theorem,

$$y_p=\Big[1+p\Delta+\tfrac{p(p-1)}{2!}\Delta^2+\tfrac{p(p-1)(p-2)}{3!}\Delta^3+\cdots\Big]y_0,$$

which is the required formula.

When to prefer it. Newton's forward formula should be used to interpolate near the beginning of the table (small $p$, $x$ close to $x_0$). Newton's backward formula is preferred near the end of the table (close to $x_n$). Since 1976 is close to the first tabulated year 1971, the forward formula is appropriate here.

(b) Population in 1976

With $x_0=1971,\ h=10$, the forward difference table is:

Here $p=\dfrac{1976-1971}{10}=0.5$. Applying the formula:

$$y=46+(0.5)(20)+\frac{(0.5)(-0.5)}{2}(-5)+\frac{(0.5)(-0.5)(-1.5)}{6}(2)+\frac{(0.5)(-0.5)(-1.5)(-2.5)}{24}(-3)$$ $$y=46+10+0.625+0.125+0.1172=56.867.$$ $$\boxed{\text{Population in 1976}\approx 56.87\ \text{thousand}\ (\approx 56{,}867).}$$

(a) Fourth-Order Runge-Kutta Method

For $\dfrac{dy}{dx}=f(x,y),\ y(x_0)=y_0$, the RK4 method estimates the increment in $y$ over one step $h$ as a weighted average of four slopes:

$$k_1=h\,f(x_n,y_n),$$ $$k_2=h\,f\!\left(x_n+\tfrac{h}{2},\,y_n+\tfrac{k_1}{2}\right),$$ $$k_3=h\,f\!\left(x_n+\tfrac{h}{2},\,y_n+\tfrac{k_2}{2}\right),$$ $$k_4=h\,f\!\left(x_n+h,\,y_n+k_3\right),$$ $$\boxed{y_{n+1}=y_n+\tfrac{1}{6}\,(k_1+2k_2+2k_3+k_4).}$$

Idea of derivation. Write $y_{n+1}=y_n+\sum w_i k_i$ and expand each $k_i$ and the exact solution in Taylor series in $h$. Matching the coefficients of $h,h^2,h^3,h^4$ gives the standard weights $\tfrac16,\tfrac13,\tfrac13,\tfrac16$ and the half-step evaluation points. The local truncation error is $O(h^{5})$, so the method is fourth order accurate.

(b) Solve $\frac{dy}{dx}=x+y^{2},\ y(0)=1,\ h=0.1$

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

$$k_1=0.1\,f(0,1)=0.1(1)=0.1000$$ $$k_2=0.1\,f(0.05,1.05)=0.1(0.05+1.1025)=0.1153$$ $$k_3=0.1\,f(0.05,1.0577)=0.1(0.05+1.1187)=0.1169$$ $$k_4=0.1\,f(0.1,1.1169)=0.1(0.1+1.2475)=0.1347$$ $$y_1=1+\tfrac16(0.1+2(0.1153)+2(0.1169)+0.1347)=1.1165.$$

So $y(0.1)\approx 1.1165$.

Step 2 ($x_1=0.1,\ y_1=1.1165$):

$$k_1=0.1347,\quad k_2=0.1551,\quad k_3=0.1576,\quad k_4=0.1823$$ $$y_2=1.1165+\tfrac16(0.1347+2(0.1551)+2(0.1576)+0.1823)=1.2736.$$ $$\boxed{y(0.2)\approx 1.2736.}$$

(a) Gauss-Seidel Method and Diagonal Dominance

For a system $Ax=b$, solve the $i$-th equation for $x_i$:

$$x_i^{(k+1)}=\frac{1}{a_{ii}}\Big(b_i-\sum_{j<i}a_{ij}x_j^{(k+1)}-\sum_{j>i}a_{ij}x_j^{(k)}\Big).$$

Unlike Jacobi's method, Gauss-Seidel uses each newly computed value immediately within the same iteration, which usually speeds up convergence. Iterations continue until successive values agree to the required accuracy.

Diagonal dominance. The matrix is (strictly) diagonally dominant if for every row

$$|a_{ii}|>\sum_{j\neq i}|a_{ij}|.$$

When this holds, the magnitude of the iteration's correction terms is bounded by a factor $<1$, so the error contracts every sweep; hence the iteration converges to the unique solution regardless of the starting guess. If a system is not dominant, rearranging the equations to achieve dominance is advisable.

(b) Solve by Gauss-Seidel

The system is already diagonally dominant. Rearranged:

$$x=\frac{9-2y-z}{10},\quad y=\frac{-44-2x+2z}{20},\quad z=\frac{22+2x-3y}{10}.$$

Starting with $x=y=z=0$:

The values are stable to three decimal places:

$$\boxed{x=1.000,\quad y=-2.000,\quad z=3.000.}$$

Secant Method for $f(x)=\cos x - x e^{x}=0$

Formula: $\displaystyle x_{n+1}=x_n-f(x_n)\,\frac{x_n-x_{n-1}}{f(x_n)-f(x_{n-1})}.$

With $x_0=0,\ x_1=1$: $f(0)=1,\ f(1)=\cos 1 - e=0.5403-2.7183=-2.1780.$

Successive values agree to three decimals at $x\approx 0.5177$–$0.5178$.

$$\boxed{\text{Root }x\approx 0.518.}$$

Fitting $y=a e^{bx}$ by Least Squares

Take logarithms: $\ln y=\ln a + b x$. Let $Y=\ln y,\ A=\ln a$, so $Y=A+bx$ — a straight line fitted by least squares.

With $n=5,\ \sum x=15$:

$$b=\frac{n\sum xY-\sum x\sum Y}{n\sum x^2-(\sum x)^2}=\frac{5(27.4443)-15(7.4860)}{5(55)-225}=\frac{24.9315}{50}=0.4986,$$ $$A=\frac{\sum Y-b\sum x}{n}=\frac{7.4860-0.4986(15)}{5}=0.0013\ \Rightarrow\ a=e^{0.0013}\approx 1.0013\approx 1.00.$$

Hence the fitted curve is

$$\boxed{y\approx 1.00\,e^{0.499\,x}.}$$

Estimate at $x=6$: $y=1.0013\,e^{0.4986\times 6}=1.0013\,e^{2.9916}\approx 19.95.$

$$\boxed{y(6)\approx 19.95.}$$

Evaluating $\displaystyle\int_0^1\frac{dx}{1+x^2}$ with $n=6,\ h=\tfrac16$

$$x_i:\ 0,\ \tfrac16,\ \tfrac26,\ \tfrac36,\ \tfrac46,\ \tfrac56,\ 1.$$

(a) Trapezoidal rule.

$$I_T=\frac{h}{2}\big[(y_0+y_6)+2(y_1+y_2+y_3+y_4+y_5)\big]=\frac{1/6}{2}\big[1.5+2(3.955445)\big]=0.784241.$$

(b) Simpson's 1/3 rule.

$$I_S=\frac{h}{3}\big[(y_0+y_6)+4(y_1+y_3+y_5)+2(y_2+y_4)\big]$$ $$=\frac{1/6}{3}\big[1.5+4(2.363137)+2(1.592308)\big]=0.785398.$$

Exact value: $\displaystyle\int_0^1\frac{dx}{1+x^2}=\tan^{-1}1=\frac{\pi}{4}=0.785398.$

Comment. Simpson's rule is far more accurate: it fits parabolic arcs (error $O(h^4)$) and matches the exact value to six decimals, whereas the trapezoidal rule (error $O(h^2)$) is correct only to about two decimals.

Numerical Differentiation at $x=1.1$

The nodes are equally spaced with $h=0.2$ and $x_0=1.0$. Since $1.1$ is near the start, use Newton's forward-difference derivative formulae with $p=\dfrac{1.1-1.0}{0.2}=0.5$.

Forward difference table:

First derivative:

$$\frac{dy}{dx}=\frac{1}{h}\Big[\Delta y_0+\frac{2p-1}{2}\Delta^2 y_0+\frac{3p^2-6p+2}{6}\Delta^3 y_0+\cdots\Big].$$

With $p=0.5$ ($\Delta^4=0$):

$$\frac{dy}{dx}=\frac{1}{0.2}\Big[0.128+\frac{0}{2}(0.288)+\frac{3(0.25)-3+2}{6}(0.048)\Big]=\frac{1}{0.2}\big[0.128-0.002\big]=\frac{0.126}{0.2}=0.63.$$

Second derivative:

$$\frac{d^2y}{dx^2}=\frac{1}{h^2}\Big[\Delta^2 y_0+(p-1)\Delta^3 y_0+\frac{12p^2-36p+22}{24}\Delta^4 y_0+\cdots\Big].$$

With $p=0.5$: $\dfrac{d^2y}{dx^2}=\dfrac{1}{0.04}\big[0.288+(-0.5)(0.048)\big]=\dfrac{0.264}{0.04}=6.6.$

$$\boxed{y'(1.1)\approx 0.63,\qquad y''(1.1)\approx 6.6.}$$

(a) Truncation vs Round-off Error

Truncation error arises from approximating an infinite/exact mathematical process by a finite one — i.e. from terminating a series or replacing a derivative/integral by a finite formula. Example: using $\sin x\approx x-\dfrac{x^3}{6}$ (dropping higher terms) introduces a truncation error of order $x^5/120$.

Round-off error arises because a computer/calculator stores numbers with a finite number of digits, so each stored value or arithmetic result is rounded. Example: representing $\tfrac13=0.3333$ to 4 digits introduces a round-off error of about $0.000033$.

(b) Error in $u=\dfrac{5xy^{2}}{z^{3}}$

The maximum absolute error is

$$\delta u=\Big|\frac{\partial u}{\partial x}\Big|\delta x+\Big|\frac{\partial u}{\partial y}\Big|\delta y+\Big|\frac{\partial u}{\partial z}\Big|\delta z.$$

Partial derivatives:

$$\frac{\partial u}{\partial x}=\frac{5y^2}{z^3},\quad \frac{\partial u}{\partial y}=\frac{10xy}{z^3},\quad \frac{\partial u}{\partial z}=-\frac{15xy^2}{z^4}.$$

At $x=y=z=1$: $u=5$, and the partials are $5,\ 10,\ 15$ in magnitude. With $\delta x=\delta y=\delta z=0.001$:

$$\delta u=(5+10+15)(0.001)=30(0.001)=0.03.$$

Maximum relative error:

$$\frac{\delta u}{u}=\frac{0.03}{5}=0.006\ (=0.6\%).$$

(Equivalently, since $\dfrac{\delta u}{u}=\dfrac{\delta x}{x}+2\dfrac{\delta y}{y}+3\dfrac{\delta z}{z}=0.001+0.002+0.003=0.006.$)

$$\boxed{\delta u_{\max}=0.03,\qquad \left(\frac{\delta u}{u}\right)_{\max}=0.006.}$$

Modified Euler Method for $\frac{dy}{dx}=x+y,\ y(0)=1,\ h=0.1$

Predictor (Euler): $y_{n+1}^{(0)}=y_n+h\,f(x_n,y_n).$

Corrector: $y_{n+1}^{(k+1)}=y_n+\dfrac{h}{2}\big[f(x_n,y_n)+f(x_{n+1},y_{n+1}^{(k)})\big]$, iterated until convergence.

Step 1: find $y(0.1)$, with $x_0=0,\ y_0=1,\ f(x_0,y_0)=1$

• Predictor: $y_1^{(0)}=1+0.1(1)=1.1000.$
• Corrector 1: $y_1=1+\tfrac{0.1}{2}\big[1+(0.1+1.1)\big]=1+0.05(2.2)=1.1100.$
• Corrector 2: $y_1=1+0.05\big[1+(0.1+1.1100)\big]=1+0.05(2.2100)=1.1105.$
• Corrector 3: $y_1=1+0.05\big[1+(0.1+1.1105)\big]=1.11053.$

Converged: $\boxed{y(0.1)\approx 1.1105.}$

Step 2: find $y(0.2)$, with $x_1=0.1,\ y_1=1.1105,\ f(x_1,y_1)=1.2105$

• Predictor: $y_2^{(0)}=1.1105+0.1(1.2105)=1.2316.$
• Corrector 1: $y_2=1.1105+0.05\big[1.2105+(0.2+1.2316)\big]=1.1105+0.05(2.6421)=1.2426.$
• Corrector 2: $y_2=1.1105+0.05\big[1.2105+(0.2+1.2426)\big]=1.1105+0.05(2.6531)=1.2432.$
• Corrector 3: $y_2=1.1105+0.05\big[1.2105+(0.2+1.2432)\big]=1.2432.$

Converged: $\boxed{y(0.2)\approx 1.2432.}$

(For reference the exact solution $y=2e^{x}-x-1$ gives $y(0.1)=1.1103,\ y(0.2)=1.2428$, confirming the accuracy.)

LU (Doolittle) Decomposition

The system $x+y+z=6,\ 3x+3y+4z=20,\ 2x+y+3z=13$ has coefficient matrix

$$A=\begin{pmatrix}1&1&1\\3&3&4\\2&1&3\end{pmatrix},\qquad b=\begin{pmatrix}6\\20\\13\end{pmatrix}.$$

Note on pivoting. Direct Doolittle factorization fails at the $(2,2)$ pivot because $u_{22}=a_{22}-l_{21}u_{12}=3-3(1)=0$. We therefore interchange rows 2 and 3 (partial pivoting) and factorize the rearranged system

$$A'=\begin{pmatrix}1&1&1\\2&1&3\\3&3&4\end{pmatrix},\qquad b'=\begin{pmatrix}6\\13\\20\end{pmatrix}.$$

Doolittle factorization $A'=LU$ ($L$ unit lower-triangular):

$$L=\begin{pmatrix}1&0&0\\2&1&0\\3&0&1\end{pmatrix},\qquad U=\begin{pmatrix}1&1&1\\0&-1&1\\0&0&1\end{pmatrix}.$$

Check: $u_{11}=1,u_{12}=1,u_{13}=1;\ l_{21}=2,\ u_{22}=1-2=-1,\ u_{23}=3-2=1;\ l_{31}=3,\ l_{32}=(3-3)/(-1)=0,\ u_{33}=4-3-0=1.$

Forward substitution $Ly=b'$:

$$y_1=6,\quad y_2=13-2(6)=1,\quad y_3=20-3(6)-0=2.$$

Back substitution $Ux=y$:

$$x_3=2,\quad -x_2+x_3=1\Rightarrow x_2=1,\quad x_1+x_2+x_3=6\Rightarrow x_1=3.$$ $$\boxed{x=3,\quad y=1,\quad z=2.}$$