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

(a) Newton-Raphson iterative formula

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

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

Since the improved estimate is $x_{n+1}=x_n+h$, the Newton-Raphson formula is:

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

Order of convergence. If $\alpha$ is a simple root, the error satisfies $e_{n+1}\approx \dfrac{f''(\alpha)}{2f'(\alpha)}e_n^{2}$, so the method is quadratically convergent (order 2).

Two situations in which it fails:

1. When $f'(x_n)=0$ (or is very small) at an iterate — the tangent is horizontal and $x_{n+1}\to\infty$.
2. When the initial guess is poor or near a turning point/inflection — iterates may oscillate or diverge (it also converges only linearly at a multiple root).

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

Here $f(x)=x^3-2x-5$, $f'(x)=3x^2-2$. With $x_0=2$:

So the root, correct to four decimals, is $\;x \approx 2.0946$, reached in 3 iterations.

Comparison with Bisection on $[2,3]$. Bisection needs $n$ such that $\dfrac{b-a}{2^{n}}\le 10^{-4}$, i.e. $\dfrac{1}{2^{n}}\le 10^{-4}\Rightarrow 2^{n}\ge 10^{4}\Rightarrow n\ge \log_2 10^{4}\approx 13.3$. Hence Bisection requires about 14 iterations.

Conclusion: Newton-Raphson reaches the same accuracy in ~3 iterations versus ~14 for Bisection, reflecting its quadratic versus linear convergence.

(a) Newton's divided difference table

Data: $x=0,1,2,4,5$, $f=2,3,12,147,326$.

Sample computations: $f[0,1]=\frac{3-2}{1-0}=1$, $f[1,2]=\frac{12-3}{2-1}=9$, $f[0,1,2]=\frac{9-1}{2-0}=4$, etc.

Newton's divided-difference polynomial:

$$f(x)=2 + 1\,(x-0) + 4\,(x-0)(x-1) + 3.875\,(x-0)(x-1)(x-2) + 0.10833\,(x-0)(x-1)(x-2)(x-4).$$

Estimate $f(3)$: with $x=3$,

$$f(3)=2 + (3) + 4(3)(2) + 3.875(3)(2)(1) + 0.10833(3)(2)(1)(-1).$$ $$f(3)=2 + 3 + 24 + 23.25 - 0.65 \approx \boxed{51.6}.$$

(b) Why divided differences are preferred when a point is added

When a new data point is appended, Newton's divided-difference formula only requires one extra term to be computed and added — all previously computed differences and coefficients remain valid:

$$P_{n+1}(x)=P_n(x)+f[x_0,\dots,x_{n+1}]\prod_{i=0}^{n}(x-x_i).$$

In contrast, Lagrange's formula must be rebuilt entirely, because every basis polynomial $L_i(x)$ contains the factor $(x-x_{new})$ and every denominator changes, forcing a full recomputation.

Illustration: Given points at $x=0,1,2$ we have a quadratic; adding $x=4$ in Newton's form just adds the term $3.875\,(x)(x-1)(x-2)$, whereas Lagrange requires recomputing all four basis polynomials from scratch. Hence Newton's form is more economical and incremental.

(a) Composite Simpson's 1/3 rule

Divide $[a,b]$ into an even number $n$ of equal sub-intervals of width $h=\dfrac{b-a}{n}$, with nodes $x_i=a+ih$. Over each pair of strips $[x_{2k},x_{2k+2}]$, approximating $f$ by a parabola gives

$$\int_{x_{2k}}^{x_{2k+2}} f\,dx \approx \frac{h}{3}\big(f_{2k}+4f_{2k+1}+f_{2k+2}\big).$$

Summing over all pairs:

$$\boxed{\int_a^b f\,dx \approx \frac{h}{3}\Big[f_0+f_n + 4\!\!\sum_{i\,\text{odd}}\!\! f_i + 2\!\!\sum_{i\,\text{even}}^{i\neq 0,n}\!\! f_i\Big].}$$

Error term. The composite error is

$$E = -\frac{(b-a)h^{4}}{180}\,f^{(4)}(\xi),\qquad a<\xi<b,$$

so the rule is $O(h^4)$ and is exact for polynomials up to degree 3.

(b) Evaluate $\displaystyle\int_0^6 \frac{dx}{1+x^2}$, $n=6$, $h=1$

Nodes and values of $f(x)=\frac{1}{1+x^2}$:

(i) Trapezoidal rule

$$I_T = h\Big[\tfrac{1}{2}(f_0+f_6) + (f_1+f_2+f_3+f_4+f_5)\Big] = 1\big[0.51351 + 0.89728\big] \approx 1.41080.$$

(ii) Simpson's 1/3 rule

$$I_S = \frac{1}{3}\big[(f_0+f_6) + 4(f_1+f_3+f_5) + 2(f_2+f_4)\big] = \frac{1}{3}\big[1.02703 + 4(0.63846) + 2(0.25882)\big] \approx 1.36617.$$

Exact value: $\int_0^6 \frac{dx}{1+x^2}=\tan^{-1}6 \approx 1.40565$.

Comment: Although Simpson's rule is normally more accurate, here the integrand $\frac{1}{1+x^2}$ varies very rapidly near $x=0$ where the coarse $h=1$ grid samples it poorly, so with only 6 strips the Trapezoidal result happens to lie closer to the exact value. Refining the mesh (larger $n$/smaller $h$) makes Simpson's $O(h^4)$ result converge far faster than the Trapezoidal $O(h^2)$ result.

(a) Fourth-order Runge-Kutta method

For $\dfrac{dy}{dx}=f(x,y)$, $y(x_0)=y_0$, with step size $h$, the RK4 method advances $y_n\to y_{n+1}$ using four slope estimates:

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

The weighted average gives the increment:

$$\boxed{\,y_{n+1} = y_n + \tfrac{1}{6}\big(k_1 + 2k_2 + 2k_3 + k_4\big)\,}$$

It is a single-step, self-starting method of local error $O(h^5)$ and global accuracy $O(h^4)$.

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

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

$$k_1=0.1(0+1)=0.10000,\quad k_2=0.1(0.05+1.05)=0.11000,$$ $$k_3=0.1(0.05+1.055)=0.11050,\quad k_4=0.1(0.1+1.1105)=0.12105.$$ $$y_1 = 1 + \tfrac{1}{6}(0.10000+0.22000+0.22100+0.12105) = 1.110342.$$

So $y(0.1)\approx 1.11034$.

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

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

(Exact: $y=2e^{x}-x-1$, $y(0.2)=1.24281$ — agreement to 4 decimals.)

Gauss elimination with partial pivoting

System (augmented matrix):

$$\left[\begin{array}{ccc|c} 2 & 1 & 1 & 10 \\ 3 & 2 & 3 & 18 \\ 1 & 4 & 9 & 16 \end{array}\right]$$

Column 1 pivot: largest $|a_{i1}|$ is 3 (row 2) → swap R1↔R2:

$$\left[\begin{array}{ccc|c} 3 & 2 & 3 & 18 \\ 2 & 1 & 1 & 10 \\ 1 & 4 & 9 & 16 \end{array}\right]$$

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

$$\left[\begin{array}{ccc|c} 3 & 2 & 3 & 18 \\ 0 & -0.3333 & -1 & -2 \\ 0 & 3.3333 & 8 & 10 \end{array}\right]$$

Column 2 pivot: $|3.3333|>|{-0.3333}|$ → swap R2↔R3:

$$\left[\begin{array}{ccc|c} 3 & 2 & 3 & 18 \\ 0 & 3.3333 & 8 & 10 \\ 0 & -0.3333 & -1 & -2 \end{array}\right]$$

Eliminate: $R_3\leftarrow R_3-\tfrac{-0.3333}{3.3333}R_2 = R_3+0.1\,R_2$:

$$\left[\begin{array}{ccc|c} 3 & 2 & 3 & 18 \\ 0 & 3.3333 & 8 & 10 \\ 0 & 0 & -0.2 & -1 \end{array}\right]$$

Back substitution:

$$z = \frac{-1}{-0.2}=5,\qquad y=\frac{10-8(5)}{3.3333}=\frac{-30}{3.3333}=-9,\qquad x=\frac{18-2(-9)-3(5)}{3}=\frac{21}{3}=7.$$ $$\boxed{x=7,\quad y=-9,\quad z=5}$$

Why partial pivoting is necessary: It places the entry of largest magnitude on the diagonal before elimination, which (i) avoids division by a zero or very small pivot, and (ii) keeps the multipliers $\le 1$ in magnitude, suppressing the growth of round-off error and improving numerical stability.

Gauss-Seidel method

Convergence condition: The iteration converges (for any starting vector) if the coefficient matrix is diagonally dominant, i.e. for every row $|a_{ii}| \ge \sum_{j\neq i}|a_{ij}|$ with strict inequality for at least one row. Here:

• Row 1: $|20| > |1|+|{-2}| = 3$ ✓
• Row 2: $|20| > |3|+|{-1}| = 4$ ✓
• Row 3: $|20| > |2|+|{-3}| = 5$ ✓

So the system is diagonally dominant and Gauss-Seidel converges.

Iteration formulas:

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

Starting from $(0,0,0)$, using the most recent values:

The values are converging to the exact solution

$$\boxed{x=1,\quad y=-1,\quad z=1.}$$

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

The normal equations are

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

With $n=5$:

Normal equations:

$$204 = 5a + 15b,\qquad 748 = 15a + 55b.$$

Solving:

$$b = \frac{n\sum xy - \sum x\sum y}{n\sum x^2-(\sum x)^2} = \frac{5(748)-15(204)}{5(55)-15^2} = \frac{3740-3060}{275-225}=\frac{680}{50}=13.6,$$ $$a = \frac{\sum y - b\sum x}{n} = \frac{204 - 13.6(15)}{5} = \frac{204-204}{5}=0.$$

Best-fit line: $\;\boxed{y = 13.6\,x}$

Estimate at $x=6$: $\;y = 13.6(6) = \boxed{81.6}.$

Derivatives by Newton's forward-difference formula

With $h=0.1$ and $x_0=1.0$, build the forward-difference table:

At the leading point ($p=0$), the formulas are

$$f'(x_0)=\frac{1}{h}\Big[\Delta f_0 - \tfrac{1}{2}\Delta^2 f_0 + \tfrac{1}{3}\Delta^3 f_0 - \tfrac{1}{4}\Delta^4 f_0\Big],$$ $$f''(x_0)=\frac{1}{h^{2}}\Big[\Delta^2 f_0 - \Delta^3 f_0 + \tfrac{11}{12}\Delta^4 f_0\Big].$$

First derivative at $x=1.0$:

$$f'(1.0)=\frac{1}{0.1}\Big[0.414-\tfrac{1}{2}(-0.036)+\tfrac{1}{3}(0.006)-\tfrac{1}{4}(-0.002)\Big]=\frac{1}{0.1}(0.4345)= \boxed{4.345}.$$

Second derivative at $x=1.0$:

$$f''(1.0)=\frac{1}{0.01}\Big[-0.036-(0.006)+\tfrac{11}{12}(-0.002)\Big]=\frac{1}{0.01}(-0.043833)= \boxed{-4.3833}.$$

(a) Round-off vs truncation error

Round-off error arises from representing numbers with a finite number of digits (limited machine/decimal precision); the exact value is replaced by its nearest representable approximation. Example: storing $1/3 = 0.333333\ldots$ as $0.3333$ introduces a round-off error of about $3.3\times10^{-5}$.

Truncation error arises from terminating an infinite mathematical process (a series, limit, or iteration) after finitely many terms — independent of machine precision. Example: approximating $\sin x \approx x - \dfrac{x^3}{6}$ drops the remaining $+\dfrac{x^5}{120}-\cdots$ terms; the omitted tail is the truncation error.

(b) Errors in rounding $\pi$

True value $\pi = 3.14159265$, approximation $\bar\pi = 3.1416$.

Absolute error:

$$E_a = |\pi-\bar\pi| = |3.14159265 - 3.1416| = 7.35\times10^{-6}\;(\approx 0.00000735).$$

Relative error:

$$E_r = \frac{E_a}{|\pi|} = \frac{7.35\times10^{-6}}{3.14159265} \approx 2.339\times10^{-6}.$$

Percentage error:

$$E_p = E_r\times 100 \approx 2.34\times10^{-4}\,\% \;(\approx 0.000234\%).$$

(a) Secant method formula and comparison

Iterative formula:

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

Difference from Newton-Raphson: Newton-Raphson uses the exact derivative $f'(x_n)$ and needs one starting point; the Secant method replaces $f'(x_n)$ by the slope of the secant through the two previous points $\dfrac{f(x_n)-f(x_{n-1})}{x_n-x_{n-1}}$, so it needs two initial approximations and does not require $f'(x)$. Its order of convergence is $\approx 1.618$ (superlinear), lower than Newton's 2, but each step is cheaper since no derivative is evaluated.

(b) Two iterations for $f(x)=\cos x - x\,e^{x}=0$

$f(x_0)=f(0)=\cos 0 - 0 = 1.000000$, $\;f(x_1)=f(1)=\cos 1 - e = 0.540302 - 2.718282 = -2.177980$.

Iteration 1:

$$x_2 = x_1 - f(x_1)\frac{x_1-x_0}{f(x_1)-f(x_0)} = 1 - (-2.177980)\frac{1-0}{-2.177980-1} = 1 - \frac{2.177980}{3.177980}\approx 0.314665.$$

$f(x_2)=\cos(0.314665) - 0.314665\,e^{0.314665} \approx 0.519871$.

Iteration 2:

$$x_3 = x_2 - f(x_2)\frac{x_2-x_1}{f(x_2)-f(x_1)} = 0.314665 - 0.519871\frac{0.314665-1}{0.519871-(-2.177980)}.$$ $$x_3 = 0.314665 - 0.519871\cdot\frac{-0.685335}{2.697851} \approx 0.314665 + 0.132063 = 0.446728.$$

$f(x_3)\approx 0.203545$.

After two iterations $\;\boxed{x \approx 0.4467}\;$ (continuing converges to the true root $\approx 0.5178$).

Population in 1925 by Newton's forward-difference interpolation

Years are equally spaced ($h=10$), so use the forward formula with $x_0=1911$. Difference table (population in thousands):

Let $p=\dfrac{x-x_0}{h}=\dfrac{1925-1911}{10}=1.4$.

Newton's forward formula:

$$f = f_0 + p\,\Delta f_0 + \frac{p(p-1)}{2!}\Delta^2 f_0 + \frac{p(p-1)(p-2)}{3!}\Delta^3 f_0 + \frac{p(p-1)(p-2)(p-3)}{4!}\Delta^4 f_0.$$

Substituting $p=1.4$, $\Delta f_0=3,\ \Delta^2 f_0=2,\ \Delta^3 f_0=0,\ \Delta^4 f_0=3$:

$$f = 12 + 1.4(3) + \frac{(1.4)(0.4)}{2}(2) + \frac{(1.4)(0.4)(-0.6)}{6}(0) + \frac{(1.4)(0.4)(-0.6)(-1.6)}{24}(3).$$ $$f = 12 + 4.2 + 0.56 + 0 + 0.0672 \approx \boxed{16.83\text{ thousand}}.$$

Assumption: The formula assumes the tabulated data can be represented by a single polynomial through all the points and that the arguments are equally spaced; it is best applied near the beginning of the table (small $p$), which is why the forward formula is appropriate for 1925.

Modified Euler (Heun's) method

For $\dfrac{dy}{dx}=f(x,y)=x^2+y$, $y(0)=1$, $h=0.1$. Predictor (Euler) then corrector iterated:

$$y^{P}_{n+1}=y_n + h f(x_n,y_n),\qquad y^{(k+1)}_{n+1}=y_n + \tfrac{h}{2}\big[f(x_n,y_n)+f(x_{n+1},y^{(k)}_{n+1})\big].$$

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

$f(0,1)=0+1=1$. Predictor: $y^P = 1 + 0.1(1) = 1.1$.

• Corrector 1: $y = 1 + 0.05[\,1 + (0.1^2+1.1)\,]=1+0.05(2.11)=1.10550$.
• Corrector 2: $y = 1 + 0.05[\,1 + (0.01+1.10550)\,]=1+0.05(2.11550)=1.10578$.
• Corrector 3: $y = 1 + 0.05[\,1 + (0.01+1.10578)\,]=1.10579$.

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

Step 2: $x_1=0.1,\,y_1=1.10579$

$f(0.1,1.10579)=0.01+1.10579=1.11579$. Predictor: $y^P = 1.10579 + 0.1(1.11579)=1.21737$.

• Corrector 1: $y = 1.10579 + 0.05[\,1.11579 + (0.2^2+1.21737)\,]=1.10579+0.05(2.37316)=1.22445$.
• Corrector 2: $y = 1.10579 + 0.05[\,1.11579 + (0.04+1.22445)\,]=1.10579+0.05(2.38024)=1.22480$.
• Corrector 3: $y = 1.10579 + 0.05[\,1.11579+(0.04+1.22480)\,]=1.22482$.

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

(Exact solution $y=3e^{x}-x^2-2x-2$ gives $y(0.1)=1.10571$, $y(0.2)=1.22461$ — close agreement.)