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

(a) Newton-Raphson Method

Derivation. Let $x_n$ be an approximation to a root of $f(x)=0$ and let $h$ be the correction so that $f(x_n+h)=0$. Expanding by Taylor's series and retaining only first-order terms:

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

Hence the iterative formula is

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

Order of convergence. Writing the error $\varepsilon_{n+1}=x_{n+1}-\alpha$ about the simple root $\alpha$ gives

$$\varepsilon_{n+1}\approx \frac{f''(\alpha)}{2f'(\alpha)}\,\varepsilon_n^{2}.$$

Thus the method converges quadratically (order $= 2$) for a simple root.

Two situations in which the method fails:

1. $f'(x_n)\to 0$ (near a turning point / stationary point): the correction $f/f'$ becomes very large and the iterate is thrown far away, so the method diverges.
2. Poor initial guess / oscillation: if $x_0$ is not close enough to the root, or where the curve has a near-horizontal tangent or an inflection, the iterates may oscillate between two values or diverge. (It also fails / loses its quadratic rate at a multiple root, where $f'(\alpha)=0$.)

(b) Bisection Method for $x^3-4x-9=0$

Let $f(x)=x^3-4x-9$. Since $f(2)=-9<0$ and $f(3)=6>0$, a root lies in $(2,3)$.

Continuing until the interval $<10^{-3}$, the iterates settle at

$$\boxed{x \approx 2.706}$$

correct to three decimal places.

(a) Newton's Divided Difference Formula

For data $(x_0,f_0),(x_1,f_1),\dots,(x_n,f_n)$ the divided differences are defined recursively by

$$f[x_i]=f_i,\qquad f[x_i,x_{i+1}]=\frac{f[x_{i+1}]-f[x_i]}{x_{i+1}-x_i},$$ $$f[x_i,\dots,x_{i+k}]=\frac{f[x_{i+1},\dots,x_{i+k}]-f[x_i,\dots,x_{i+k-1}]}{x_{i+k}-x_i}.$$

Interpolating polynomial. Assume

$$P(x)=a_0+a_1(x-x_0)+a_2(x-x_0)(x-x_1)+\cdots+a_n(x-x_0)\cdots(x-x_{n-1}).$$

Putting $x=x_0$ gives $a_0=f[x_0]$; putting $x=x_1$ gives $a_1=f[x_0,x_1]$; in general $a_k=f[x_0,x_1,\dots,x_k]$. Hence Newton's divided difference formula:

$$\boxed{P(x)=f[x_0]+(x-x_0)f[x_0,x_1]+(x-x_0)(x-x_1)f[x_0,x_1,x_2]+\cdots}$$

Symmetry. By induction it can be shown that

$$f[x_0,x_1,\dots,x_n]=\sum_{i=0}^{n}\frac{f_i}{\prod_{j\neq i}(x_i-x_j)}.$$

This expression is a symmetric function of $x_0,x_1,\dots,x_n$ — interchanging any two arguments leaves it unchanged. Therefore the divided differences are symmetric with respect to their arguments, i.e.

$$f[x_0,x_1]=f[x_1,x_0],\quad f[x_0,x_1,x_2]=f[x_2,x_0,x_1]=\dots$$

(b) Lagrange Interpolation at $x=10$

Data: $(5,12),(6,13),(9,14),(11,16)$. Lagrange's formula:

$$f(x)=\sum_{i}f_i\prod_{j\neq i}\frac{x-x_j}{x_i-x_j}.$$

At $x=10$:

$$L_0=\frac{(10-6)(10-9)(10-11)}{(5-6)(5-9)(5-11)}=\frac{(4)(1)(-1)}{(-1)(-4)(-6)}=\frac{-4}{-24}=\tfrac{1}{6}$$ $$L_1=\frac{(10-5)(10-9)(10-11)}{(6-5)(6-9)(6-11)}=\frac{(5)(1)(-1)}{(1)(-3)(-5)}=\frac{-5}{15}=-\tfrac{1}{3}$$ $$L_2=\frac{(10-5)(10-6)(10-11)}{(9-5)(9-6)(9-11)}=\frac{(5)(4)(-1)}{(4)(3)(-2)}=\frac{-20}{-24}=\tfrac{5}{6}$$ $$L_3=\frac{(10-5)(10-6)(10-9)}{(11-5)(11-6)(11-9)}=\frac{(5)(4)(1)}{(6)(5)(2)}=\frac{20}{60}=\tfrac{1}{3}$$ $$f(10)=12\cdot\tfrac16+13\cdot(-\tfrac13)+14\cdot\tfrac56+16\cdot\tfrac13 = 2-\tfrac{13}{3}+\tfrac{35}{3}+\tfrac{16}{3}=2+\tfrac{38}{3}.$$ $$\boxed{f(10)=\tfrac{44}{3}\approx 14.667}$$

(a) Gauss-Seidel Iterative Method

For a system $a_{11}x_1+\dots=b_1,\dots$, each equation is solved for its diagonal unknown:

$$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).$$

The distinguishing feature is that each newly computed value is used immediately in the same sweep, so it generally converges faster than Jacobi iteration. Iteration continues until successive iterates agree to the required tolerance.

Condition of convergence (diagonal dominance). A sufficient condition is that the coefficient matrix be diagonally dominant:

$$|a_{ii}|\ge\sum_{j\neq i}|a_{ij}|\quad\text{for all }i,\text{ with strict inequality for at least one }i.$$

Under this condition the iteration converges for any starting vector.

(b) Three Iterations, start $(0,0,0)$

The system is already diagonally dominant. Rearranging:

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

Iteration 1 (with $x=y=z=0$ initially):

$$x=\tfrac{17-0+0}{20}=0.8500,\;\; y=\tfrac{-18-3(0.85)+0}{20}=-1.0275,\;\; z=\tfrac{25-2(0.85)+3(-1.0275)}{20}=1.0109.$$

Iteration 2:

$$x=\tfrac{17-(-1.0275)+2(1.0109)}{20}=1.0025,\;\; y=\tfrac{-18-3(1.0025)+1.0109}{20}=-0.9998,\;\; z=\tfrac{25-2(1.0025)+3(-0.9998)}{20}=0.9998.$$

Iteration 3:

$$x=\tfrac{17-(-0.9998)+2(0.9998)}{20}=1.0000,\;\; y=\tfrac{-18-3(1.0000)+0.9998}{20}=-1.0000,\;\; z=\tfrac{25-2(1.0000)+3(-1.0000)}{20}=1.0000.$$ $$\boxed{x=1,\; y=-1,\; z=1}$$

(a) Five-Point Formula for the Laplace Equation

Use a uniform square mesh of spacing $h$ ($\Delta x=\Delta y=h$); let $u_{i,j}=u(x_i,y_j)$. Central second differences give

$$\frac{\partial^2 u}{\partial x^2}\Big|_{i,j}\approx\frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{h^2},\qquad \frac{\partial^2 u}{\partial y^2}\Big|_{i,j}\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,$$ $$\boxed{\,u_{i,j}=\tfrac14\big(u_{i+1,j}+u_{i-1,j}+u_{i,j+1}+u_{i,j-1}\big)\,}.$$

This standard five-point (cross) formula states that the value at any interior node equals the average of its four nearest neighbours (Liebmann's method solves the resulting set iteratively).

(b) Classification of $A\,u_{xx}+B\,u_{xy}+C\,u_{yy}+f=0$

The type is fixed by the discriminant $\Delta=B^2-4AC$:

Worked check: For Laplace, $A=1,B=0,C=1\Rightarrow B^2-4AC=-4<0$ (elliptic). For the heat equation $u_{xx}-u_t=0$ written in two variables, the second-order part has $A=1,B=0,C=0\Rightarrow B^2-4AC=0$ (parabolic). For the wave equation $u_{tt}-c^2u_{xx}=0$, $B^2-4AC=4c^2>0$ (hyperbolic).

Secant Method for $\cos x = 3x-1$

Write $f(x)=\cos x-3x+1=0$. The secant iteration is

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

With $x_0=0.5,\;x_1=0.6$:

$$f(0.5)=\cos 0.5-1.5+1=0.87758-0.5=0.37758,$$ $$f(0.6)=\cos 0.6-1.8+1=0.82534-0.8=0.02534.$$

Sample step (n=1):

$$x_2=0.6-\frac{0.025336\,(0.6-0.5)}{0.025336-0.37758}=0.6-\frac{0.0025336}{-0.35224}=0.60719.$$

The successive iterates agree to four significant figures, so

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

Numerical Integration of $\displaystyle\int_0^6\frac{dx}{1+x^2}$

With 6 subintervals, $h=1$. Tabulate $y=\dfrac{1}{1+x^2}$:

(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}{2}\big[(1.0+0.027027)+2(0.5+0.2+0.1+0.058824+0.038462)\big]=\frac{1}{2}[1.027027+1.794571]=\boxed{1.4108}.$$

(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}{3}\big[1.027027+4(0.5+0.1+0.038462)+2(0.2+0.058824)\big]=\frac{1}{3}[1.027027+2.553848+0.517648]=\boxed{1.3662}.$$

Comparison with exact value

$$\int_0^6\frac{dx}{1+x^2}=[\tan^{-1}x]_0^6=\tan^{-1}6=1.40565.$$

The trapezoidal value is closer here because the integrand changes rapidly near $x=0$; Simpson's rule, though normally more accurate, is more affected by the coarse mesh over this steep region.

Numerical Differentiation by Newton's Forward Formula

Here $x_0=1.0$, $h=0.2$. The forward difference table:

Third differences are constant ($0.048$) and higher differences vanish.

For $x=1.1$: $p=\dfrac{x-x_0}{h}=\dfrac{1.1-1.0}{0.2}=0.5$.

First derivative

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

With $p=0.5$: $\dfrac{2p-1}{2}=0$, $\dfrac{3p^2-6p+2}{6}=\dfrac{0.75-3+2}{6}=-0.041667$.

$$\frac{dy}{dx}=\frac{1}{0.2}\big[0.128+0+(-0.041667)(0.048)\big]=\frac{1}{0.2}(0.126)=\boxed{0.63}.$$

Second derivative

$$\frac{d^2y}{dx^2}=\frac{1}{h^2}\Big[\Delta^2y_0+(p-1)\Delta^3y_0\Big]=\frac{1}{0.04}\big[0.288+(0.5-1)(0.048)\big]=\frac{1}{0.04}(0.264)=\boxed{6.6}.$$

Thus at $x=1.1$, $f'(x)\approx 0.63$ and $f''(x)\approx 6.6$.

Least-Squares Straight Line $y=a+bx$

Normal equations:

$$\sum y=na+b\sum x,\qquad \sum xy=a\sum x+b\sum x^2.$$

with $\sum x=15,\;n=5$. The normal equations:

$$17.1=5a+15b,\qquad 64.2=15a+55b.$$

Solving (e.g. $b=\dfrac{n\sum xy-\sum x\sum y}{n\sum x^2-(\sum x)^2}=\dfrac{5(64.2)-15(17.1)}{5(55)-225}=\dfrac{321-256.5}{50}=\dfrac{64.5}{50}=1.29$):

$$b=1.29,\qquad a=\frac{\sum y-b\sum x}{n}=\frac{17.1-1.29(15)}{5}=\frac{17.1-19.35}{5}=-0.45.$$

Fitted line:

$$\boxed{y=-0.45+1.29x}.$$

Estimate at $x=2.5$:

$$y=-0.45+1.29(2.5)=-0.45+3.225=\boxed{2.775}.$$

Fourth-Order Runge-Kutta Method

Given $\dfrac{dy}{dx}=f(x,y)=x+y^2$, $y(0)=1$, $h=0.1$. RK4:

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

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

$$k_1=0.1(0+1^2)=0.100000$$ $$k_2=0.1\big(0.05+(1.05)^2\big)=0.1(0.05+1.1025)=0.115250$$ $$k_3=0.1\big(0.05+(1.057625)^2\big)=0.1(0.05+1.118570)=0.116857$$ $$k_4=0.1\big(0.1+(1.116857)^2\big)=0.1(0.1+1.247370)=0.134737$$ $$y_1=1+\tfrac16(0.1+2(0.11525)+2(0.116857)+0.134737)=1+0.116492=1.116492.$$

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

$$k_1=0.1\big(0.1+(1.116492)^2\big)=0.134655$$ $$k_2=0.1\big(0.15+(1.183819)^2\big)=0.155143$$ $$k_3=0.1\big(0.15+(1.194063)^2\big)=0.157579$$ $$k_4=0.1\big(0.2+(1.274071)^2\big)=0.182326$$ $$y_2=1.116492+\tfrac16(0.134655+2(0.155143)+2(0.157579)+0.182326)=1.116492+0.157071=1.273563.$$ $$\boxed{y(0.2)\approx 1.2736}.$$

LU (Doolittle) Decomposition

$$A=\begin{pmatrix}2&3&1\\1&2&3\\3&1&2\end{pmatrix},\quad b=\begin{pmatrix}9\\6\\8\end{pmatrix}.$$

In Doolittle's method $A=LU$ with unit diagonal in $L$:

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

Row 1 of U: $u_{11}=2,\;u_{12}=3,\;u_{13}=1.$ Column 1 of L: $l_{21}=1/2=0.5,\;l_{31}=3/2=1.5.$ Row 2 of U: $u_{22}=2-0.5(3)=0.5,\;u_{23}=3-0.5(1)=2.5.$ $l_{32}$: $l_{32}=\dfrac{1-1.5(3)}{0.5}=\dfrac{1-4.5}{0.5}=-7.$ $u_{33}$: $u_{33}=2-1.5(1)-(-7)(2.5)=2-1.5+17.5=18.$

$$L=\begin{pmatrix}1&0&0\\0.5&1&0\\1.5&-7&1\end{pmatrix},\qquad U=\begin{pmatrix}2&3&1\\0&0.5&2.5\\0&0&18\end{pmatrix}.$$

Forward substitution $Lz=b$:

$$z_1=9,\quad z_2=6-0.5(9)=1.5,\quad z_3=8-1.5(9)-(-7)(1.5)=8-13.5+10.5=5.$$

Back substitution $Ux=z$:

$$x_3=\frac{5}{18}=0.2778,$$ $$x_2=\frac{1.5-2.5(0.2778)}{0.5}=\frac{1.5-0.6944}{0.5}=1.6111,$$ $$x_1=\frac{9-3(1.6111)-1(0.2778)}{2}=\frac{9-4.8333-0.2778}{2}=1.9444.$$ $$\boxed{x=1.9444,\;y=1.6111,\;z=0.2778}$$

(equivalently $x=\tfrac{35}{18},\,y=\tfrac{29}{18},\,z=\tfrac{5}{18}$).

Bender-Schmidt Explicit Scheme for the Heat Equation

The explicit finite-difference form of $u_t=u_{xx}$ is

$$u_{i}^{j+1}=u_i^{j}+r\big(u_{i-1}^{j}-2u_i^{j}+u_{i+1}^{j}\big),\qquad r=\frac{k}{h^2},$$

where $h=\Delta x$, $k=\Delta t$. Stability requires $r\le\tfrac12$. Choosing the Bender-Schmidt value $r=\tfrac12$ gives the simplified average formula

$$\boxed{u_i^{j+1}=\tfrac12\big(u_{i-1}^{j}+u_{i+1}^{j}\big)}.$$

With $h=0.25$, $r=\tfrac12\Rightarrow k=r\,h^2=0.5(0.0625)=0.03125.$ Mesh nodes $x=0,0.25,0.5,0.75,1$.

Initial values $u(x,0)=\sin\pi x$, with boundaries $u_0=u_4=0$ for all $t$:

$$u_0=0,\;u_1=\sin\tfrac\pi4=0.7071,\;u_2=\sin\tfrac\pi2=1.0000,\;u_3=\sin\tfrac{3\pi}4=0.7071,\;u_4=0.$$

Level $j=1$ ($t=0.03125$):

$$u_1=\tfrac12(u_0+u_2)=\tfrac12(0+1.0)=0.5000$$ $$u_2=\tfrac12(u_1+u_3)=\tfrac12(0.7071+0.7071)=0.7071$$ $$u_3=\tfrac12(u_2+u_4)=\tfrac12(1.0+0)=0.5000$$

Level $j=2$ ($t=0.0625$):

$$u_1=\tfrac12(0+0.7071)=0.3536,\quad u_2=\tfrac12(0.5+0.5)=0.5000,\quad u_3=\tfrac12(0.7071+0)=0.3536.$$

The solution decays symmetrically, consistent with the exact solution $u=e^{-\pi^2 t}\sin\pi x$.