Numerical Methods for Algebraic and Transcendental Equations

Error in Numerical Methods

Definition. In numerical computation we replace the exact (true) value $x$ of a quantity by an approximate value $\bar{x}$. The difference between the two is called the error:

$E = x - \bar{x} \quad (\text{true value} - \text{approximate value})$

Kinds of error

Numerical part

Given: true value $x = 625.4835$, to be approximated to 4 significant digits.

The first four significant digits are $6,\,2,\,5,\,4$; the next digit is $8 \ge 5$, so we round up:

$\bar{x} = 625.5$

Absolute error: $E_a = \lvert 625.4835 - 625.5\rvert = 0.0165$

Relative error: $E_r = \frac{0.0165}{625.4835} = 2.638\times 10^{-5}$

Percentage error: $E_p = 2.638\times 10^{-5}\times 100\% = 2.638\times 10^{-3}\% \approx 0.00264\%$

Result: The percentage error is approximately $\boxed{0.00264\%}$.

---

OR — MATLAB program for RK-2 (Heun's) method

$y' = f(x,y),\qquad y(x_0)=y_0,\qquad \text{find } y(x_p),\ \text{step } h$

Heun's RK-2 formulas: $k_1 = h\,f(x_n,y_n),\qquad k_2 = h\,f(x_n+h,\;y_n+k_1)$ $y_{n+1} = y_n + \tfrac{1}{2}(k_1+k_2),\qquad x_{n+1}=x_n+h$

% RK-2 (Heun's) method for y' = f(x,y), y(x0) = y0, at x = xp
clc; clear;

f  = @(x,y) x + y;      % <-- enter f(x,y) here
x0 = 0;                 % initial x
y0 = 1;                 % initial y
xp = 0.2;               % point where y is required
h  = 0.05;              % step size

n = round((xp - x0)/h); % number of steps
x = x0;  y = y0;

fprintf('  x        y\n');
fprintf('%6.4f  %10.6f\n', x, y);

for i = 1:n
    k1 = h*f(x, y);
    k2 = h*f(x + h, y + k1);
    y  = y + (k1 + k2)/2;
    x  = x + h;
    fprintf('%6.4f  %10.6f\n', x, y);
end

fprintf('\ny(%.4f) = %.6f\n', xp, y);

Output (for $f=x+y,\ y(0)=1,\ h=0.05$): $y(0.2) \approx 1.242609$ (exact value $2e^{0.2}-1-0.2 = 1.242806$).

Solve by LU-decomposition (Doolittle method)

Given system:

$$\begin{aligned} x+y+z&=9\\ 2x-3y+4z&=13\\ 3x+4y+5z&=40 \end{aligned}$$

In matrix form, (Ax=b):

$$A=\begin{pmatrix} 1&1&1\\ 2&-3&4\\ 3&4&5 \end{pmatrix},\qquad b=\begin{pmatrix} 9\\13\\40 \end{pmatrix}$$

Let

$$A=LU$$

where (L) is unit lower triangular and (U) is upper triangular:

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

Step 1: Find (L) and (U)

From the first row of (U):

$$u_{11}=1,\quad u_{12}=1,\quad u_{13}=1$$

From the first column of (L):

$$l_{21}=\frac{2}{1}=2,\qquad l_{31}=\frac{3}{1}=3$$

Now compute the second row of (U):

$$u_{22}=-3-l_{21}u_{12}=-3-2(1)=-5$$ $$u_{23}=4-l_{21}u_{13}=4-2(1)=2$$

Compute (l_{32}):

$$l_{32}=\frac{4-l_{31}u_{12}}{u_{22}} =\frac{4-3(1)}{-5} =-\frac15$$

Compute (u_{33}):

$$u_{33}=5-l_{31}u_{13}-l_{32}u_{23} =5-3(1)-\left(-\frac15\right)(2) =2+\frac25=\frac{12}{5}$$

Hence,

$$L=\begin{pmatrix} 1&0&0\\ 2&1&0\\ 3&-\frac15&1 \end{pmatrix},\qquad U=\begin{pmatrix} 1&1&1\\ 0&-5&2\\ 0&0&\frac{12}{5} \end{pmatrix}$$

Step 2: Solve (Ly=b)

Let

$$y=\begin{pmatrix} y_1\\y_2\\y_3 \end{pmatrix}$$

Then

$$\begin{aligned} y_1&=9\\ 2y_1+y_2&=13\\ 3y_1-\frac15y_2+y_3&=40 \end{aligned}$$

From the first equation:

$$y_1=9$$

From the second:

$$2(9)+y_2=13 \Rightarrow y_2=-5$$

From the third:

$$3(9)-\frac15(-5)+y_3=40$$ $$27+1+y_3=40 \Rightarrow y_3=12$$

Thus,

$$y=\begin{pmatrix} 9\\-5\\12 \end{pmatrix}$$

Step 3: Solve (Ux=y)

Let

$$x=\begin{pmatrix} x\\y\\z \end{pmatrix}$$

Then

$$\begin{aligned} x+y+z&=9\\ -5y+2z&=-5\\ \frac{12}{5}z&=12 \end{aligned}$$

From the last equation:

$$z=12\cdot\frac{5}{12}=5$$

From the second equation:

$$-5y+2(5)=-5$$ $$-5y+10=-5 \Rightarrow -5y=-15 \Rightarrow y=3$$

From the first equation:

$$x+3+5=9 \Rightarrow x=1$$

Therefore,

$$\boxed{x=1,\quad y=3,\quad z=5}$$

Verification

$$2(1)-3(3)+4(5)=2-9+20=13$$ $$3(1)+4(3)+5(5)=3+12+25=40$$

Hence the solution is correct.

Using Lagrange’s interpolation formula,

[ y(x)=\sum_{i=0}^{3} y_i \prod_{\substack{j=0\j\ne i}}^{3}\frac{x-x_j}{x_i-x_j} ]

Given:

[ (x_0,y_0)=(1,10),\quad (x_1,y_1)=(4,25),\quad (x_2,y_2)=(9,40),\quad (x_3,y_3)=(11,60) ]

For (x=3),

[ L_0=\frac{(3-4)(3-9)(3-11)}{(1-4)(1-9)(1-11)} =\frac{(-1)(-6)(-8)}{(-3)(-8)(-10)} =\frac{-48}{-240}=\frac15 ]

[ L_1=\frac{(3-1)(3-9)(3-11)}{(4-1)(4-9)(4-11)} =\frac{(2)(-6)(-8)}{(3)(-5)(-7)} =\frac{96}{105}=\frac{32}{35} ]

[ L_2=\frac{(3-1)(3-4)(3-11)}{(9-1)(9-4)(9-11)} =\frac{(2)(-1)(-8)}{(8)(5)(-2)} =\frac{16}{-80}=-\frac15 ]

[ L_3=\frac{(3-1)(3-4)(3-9)}{(11-1)(11-4)(11-9)} =\frac{(2)(-1)(-6)}{(10)(7)(2)} =\frac{12}{140}=\frac{3}{35} ]

Now,

[ y(3)=10L_0+25L_1+40L_2+60L_3 ]

[ =10\left(\frac15\right)+25\left(\frac{32}{35}\right)+40\left(-\frac15\right)+60\left(\frac{3}{35}\right) ]

[ =2+\frac{160}{7}-8+\frac{36}{7} ]

[ =-6+\frac{196}{7} ]

[ =-6+28=22 ]

[ \boxed{y(3)=22} ]

Secant method for $3x=\cos x+1$

Let $f(x)=3x-\cos x-1$.

Locating the root: $f(0.6)=1.8-0.82534-1=-0.02534<0$ and $f(0.62)=1.86-0.81388-1=0.04612>0$, so a root lies in $(0.6,\,0.62)$. Take $x_0=0.6$, $x_1=0.62$ (angles in radians).

Secant formula:

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

The last two iterates agree to four decimal places ($0.6071$), so the iteration has converged.

Root: $x\approx 0.6071$

Simpson's 1/3 rule for $\displaystyle\int_0^1 e^{-x^2}dx$

Take $n=4$ sub-intervals (an even number, as the rule requires), so $h=\dfrac{1-0}{4}=0.25$. Let $y=e^{-x^2}$.

$\int_0^1 e^{-x^2}dx\approx\frac h3\left[y_0+y_4+4(y_1+y_3)+2y_2\right]$

$=\frac{0.25}{3}\left[1+0.367879+4(0.939413+0.569783)+2(0.778801)\right]$

$=\frac{0.25}{3}\left[1.367879+6.036784+1.557602\right]=\frac{0.25}{3}(8.962265)\approx0.74686$

Result: $\displaystyle\int_0^1 e^{-x^2}dx\approx0.7469$ (the exact value is $0.74682$).

Euler's method for $y'=x+\sin y,\ y(0)=1,\ h=0.1$

Formula: $y_{n+1}=y_n+h\,f(x_n,y_n)=y_n+0.1\,(x_n+\sin y_n)$, with $f(x,y)=x+\sin y$ (angles in radians).

Result: $y(1)\approx 2.375$