Numerical Methods

Taking $x=1,2,3,4,5$ with $y=16,\,17.1,\,8.7,\,4.7,\,2.6$.

Linearise. From $y=ab^x$: $\log y=\log a+x\log b$, i.e. $Y=A+Bx$ with $Y=\log_{10}y$, $A=\log a$, $B=\log b$.

Normal equations: $4.4637=5A+15B$ and $11.2519=15A+55B$.

$B=\frac{5(11.2519)-15(4.4637)}{5(55)-15^2}=-0.2139,\qquad A=\frac{4.4637+15(0.2139)}{5}=1.5345$

$a=10^{1.5345}=34.24,\qquad b=10^{-0.2139}=0.611$

Fitted curve: $y=34.24\,(0.611)^x$.

At $x=7$: $y=34.24(0.611)^7=34.24\times0.0318\approx1.09$.

Taking $x=1,2,3,4$ with $y=10.5,\,2.15,\,3.85,\,8.30$.

Linearise. $y=ae^{bx}\Rightarrow\ln y=\ln a+bx$, i.e. $Y=A+bx$ with $Y=\ln y$, $A=\ln a$.

Normal equations: $6.5813=4A+10b$ and $16.3917=10A+30b$.

$b=\frac{4(16.3917)-10(6.5813)}{4(30)-10^2}=\frac{-0.2462}{20}\approx-0.0123$

$A=\frac{6.5813+10(0.0123)}{4}=1.6761,\qquad a=e^{1.6761}\approx5.34$

Fitted curve: $y\approx5.34\,e^{-0.0123x}$.

Derivation. Let $x_n$ be an approximation to a root of $f(x)=0$. Expanding $f$ about $x_n$ by Taylor's series and keeping the first two terms:

$f(x_n+h)\approx f(x_n)+h f'(x_n)$

We want $f(x_n+h)=0$, so $h=-\dfrac{f(x_n)}{f'(x_n)}$. The improved root is $x_{n+1}=x_n+h$:

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

Geometrically, $x_{n+1}$ is where the tangent to $y=f(x)$ at $(x_n,f(x_n))$ cuts the $x$-axis. The method converges quadratically near a simple root.

Numerical work. $f(x)=x^3-4x-9$, $f'(x)=3x^2-4$, $x_0=2.5$.

After four iterations the root is $x\approx 2.7065$.

Method. For $f(x)=x^3-x-1$, $f(1)=-1<0$ and $f(2)=5>0$, so a root lies in $(1,2)$. The false position (regula falsi) formula is

$c=\frac{a f(b)-b f(a)}{f(b)-f(a)}$

The interval is then replaced by $[c,b]$ if $f(c)<0$ (otherwise $[a,c]$).

After 6 iterations the root is $x\approx 1.3223$ (the true root is $1.3247$). Since $f(c)<0$ every time, $b=2$ stays fixed and the convergence is linear.

Method. Rewrite $f(x)=0$ as $x=g(x)$ and iterate $x_{n+1}=g(x_n)$. The iteration converges to a root near which $|g'(x)|<1$.

Write $x^2-5x+6=0$ as

$x=g(x)=\frac{x^2+6}{5},\qquad g'(x)=\frac{2x}{5}$

For $x<2.5$ we have $|g'(x)|<1$, so the iteration converges to the root $x=2$ (the other root $x=3$ has $g'(3)=1.2>1$ and is not reached by this form).

Take $x_0=1$:

The values approach $x=2$ steadily (linear convergence, ratio about $g'(2)=0.8$). Hence a root is $x=2$; check: $2^2-5(2)+6=0$.

Taking the tabulated data as $f(x)=\log_{10}x$ at $x=3,\,3.2,\,3.4,\,3.6$ with $f=0.4771,\,0.5051,\,0.5315,\,0.5563$.

Divided difference table

Newton's divided difference polynomial

$f(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]+\dots$

$f(x)=0.4771+0.14(x-3)-0.02(x-3)(x-3.2)$

The third difference is zero, so the polynomial is quadratic. Expanded:

$f(x)=-0.02x^2+0.264x-0.1349$

At $x=3.45$:

$f(3.45)=0.4771+0.14(0.45)-0.02(0.45)(0.25)=0.4771+0.0630-0.00225=0.53785$

So $f(3.45)\approx 0.5379$.

Taking the data as $t=1,2,3,4$ min and $T=70,83,100,124\,^\circ$C.

Linearising. Put $u=e^{t/2}$. Then $T=au+b$ is a straight line in $u$, and the normal equations are

$\sum T=a\sum u+nb,\qquad \sum uT=a\sum u^2+b\sum u$

$377=16.2377a+4b,\qquad 1705.44=84.7911a+16.2377b$

Solving: $a=9.273$, $b=56.606$, so

$T=9.273\,e^{t/2}+56.606$

Estimates

$t=1.5$: $T=9.273e^{0.75}+56.606=9.273(2.1170)+56.606\approx 76.2\,^\circ$C

$t=3.5$: $T=9.273e^{1.75}+56.606=9.273(5.7546)+56.606\approx 110.0\,^\circ$C

Taking the integral as $I=\displaystyle\int_{0.2}^{1.5}\left(\sin x-\ln x+e^x\right)dx$ with $n=6$.

$h=\frac{1.5-0.2}{6}=0.21667$

Trapezoidal rule

$I\approx\frac h2\left[y_0+y_6+2(y_1+\dots+y_5)\right]=\frac{0.21667}{2}\left[8.10323+2(17.01788)\right]\approx 4.5651$

Simpson's 1/3 rule (valid since $n=6$ is even)

$I\approx\frac h3\left[y_0+y_6+4(y_1+y_3+y_5)+2(y_2+y_4)\right]$

$=\frac{0.21667}{3}\left[8.10323+4(10.36868)+2(6.64920)\right]\approx 4.5411$

So Trapezoidal gives $4.565$ and Simpson's rule gives $4.541$ (Simpson is more accurate).

Augmented matrix:

$\left[\begin{array}{ccc|c}3&1&1&1\\2&4&1&2\\1&1&5&3\end{array}\right]$

Forward elimination. $R_2\to R_2-\tfrac23R_1$, $R_3\to R_3-\tfrac13R_1$:

$\left[\begin{array}{ccc|c}3&1&1&1\\0&\tfrac{10}{3}&\tfrac13&\tfrac43\\0&\tfrac23&\tfrac{14}{3}&\tfrac83\end{array}\right]$

$R_3\to R_3-\tfrac15R_2$:

$\left[\begin{array}{ccc|c}3&1&1&1\\0&\tfrac{10}{3}&\tfrac13&\tfrac43\\0&0&\tfrac{23}{5}&\tfrac{12}{5}\end{array}\right]$

Back substitution

$z=\frac{12/5}{23/5}=\frac{12}{23}\approx0.5217$

$y=\frac{\tfrac43-\tfrac13z}{\tfrac{10}{3}}=\frac{4-z}{10}=\frac{8}{23}\approx0.3478$

$x=\frac{1-y-z}{3}=\frac{1}{23}\approx0.0435$

Solution: $x=\tfrac1{23},\ y=\tfrac8{23},\ z=\tfrac{12}{23}$. Check in eq. 2: $\tfrac{2+32+12}{23}=2$.

The system is diagonally dominant, so Jacobi iteration converges. Solve each equation for its diagonal variable and use only the previous iterate:

$x^{(k+1)}=\frac{12-y^{(k)}-z^{(k)}}{10},\quad y^{(k+1)}=\frac{13-2x^{(k)}-z^{(k)}}{10},\quad z^{(k+1)}=\frac{14-2x^{(k)}-2y^{(k)}}{10}$

Start $(x,y,z)=(0,0,0)$.

For example iteration 2: $x=(12-1.3-1.4)/10=0.93$, $y=(13-2.4-1.4)/10=0.92$, $z=(14-2.4-2.6)/10=0.90$.

After three iterations $x\approx1.018,\ y\approx1.024,\ z\approx1.030$; the values are converging to the exact solution $(1,1,1)$.

Taking the matrix as $A=\begin{bmatrix}2&1\\3&4\end{bmatrix}$ and using the power method with $v_0=\begin{bmatrix}1\\1\end{bmatrix}$. At each step compute $Av_k$, take its largest component as the eigenvalue estimate and divide by it.

The iterates converge to

$\lambda\approx5,\qquad v\approx\begin{bmatrix}0.3333\\1\end{bmatrix}\ \text{(or }\begin{bmatrix}1\\3\end{bmatrix})$

Check: $A(1,3)^T=(5,15)^T=5(1,3)^T$. The characteristic equation $\lambda^2-6\lambda+5=0$ gives $\lambda=5,1$, so 5 is dominant.

Let $y'=z$. Then $z'=6-2z+3y$ with $y(0)=0,\ z(0)=1$, $h=0.2$.

Euler's formulas: $y_{n+1}=y_n+hz_n$, $z_{n+1}=z_n+h(6-2z_n+3y_n)$.

Step 1 ($x=0\to0.2$): $y_1=0+0.2(1)=0.2$ $z_1=1+0.2(6-2+0)=1.8$

Step 2 ($x=0.2\to0.4$): $y_2=0.2+0.2(1.8)=0.56$ $z_2=1.8+0.2(6-3.6+0.6)=2.4$

Hence $y(0.4)\approx0.56$.

$f(x,y)=\dfrac{2y}{x}$, $x_0=1,\ y_0=2$, $h=0.25$. Second-order Runge-Kutta (Heun):

$k_1=hf(x_n,y_n),\quad k_2=hf(x_n+h,y_n+k_1),\quad y_{n+1}=y_n+\tfrac12(k_1+k_2)$

Hence $y(2)\approx7.861$. (The exact solution $y=2x^2$ gives $y(2)=8$.)

Here $\nabla^2u=-10(x^2+y^2)$ on $0\le x,y\le3$ with $u=0$ on the boundary and $h=1$. The interior nodes are $(1,1),(2,1),(1,2),(2,2)$.

The five-point formula for $u_{xx}+u_{yy}=g$ gives

$u_{i+1,j}+u_{i-1,j}+u_{i,j+1}+u_{i,j-1}-4u_{i,j}=h^2g(x_i,y_j)$

so $4u_{ij}=\text{(sum of 4 neighbours)}+10(x_i^2+y_i^2)$ with right-hand sides $20,\,50,\,50,\,80$ at $(1,1),(2,1),(1,2),(2,2)$.

By symmetry $u_{21}=u_{12}$. Let $u_1=u_{11}$, $u_2=u_{21}=u_{12}$, $u_3=u_{22}$ (boundary values are 0):

$4u_1=2u_2+20,\qquad 4u_2=u_1+u_3+50,\qquad 4u_3=2u_2+80$

So $u_1=\dfrac{u_2+10}{2}$, $u_3=\dfrac{u_2+40}{2}$. Substituting:

$4u_2=\frac{u_2+10}{2}+\frac{u_2+40}{2}+50=u_2+75\ \Rightarrow\ u_2=25$

$u_1=17.5,\qquad u_3=32.5$

Hence $u(1,1)=17.5,\ u(2,1)=u(1,2)=25,\ u(2,2)=32.5$.

Cubic spline interpolation fits a separate cubic polynomial $S_i(x)$ on each subinterval $[x_i,x_{i+1}]$ of the data points $(x_0,y_0),\dots,(x_n,y_n)$, joined so the curve is smooth.

Conditions

• $S_i(x_i)=y_i$ and $S_i(x_{i+1})=y_{i+1}$ (passes through the points).
• $S_i'(x_{i+1})=S_{i+1}'(x_{i+1})$ (continuous slope).
• $S_i''(x_{i+1})=S_{i+1}''(x_{i+1})$ (continuous curvature).
• End conditions: for a natural spline $S''(x_0)=S''(x_n)=0$.

Form. With $M_i=S''(x_i)$ and $h_i=x_{i+1}-x_i$:

$S_i(x)=\frac{M_i(x_{i+1}-x)^3+M_{i+1}(x-x_i)^3}{6h_i}+\left(\frac{y_i}{h_i}-\frac{M_ih_i}{6}\right)(x_{i+1}-x)+\left(\frac{y_{i+1}}{h_i}-\frac{M_{i+1}h_i}{6}\right)(x-x_i)$

The $M_i$ come from 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]$

Advantage. It avoids the oscillation of a single high-degree polynomial and gives a smooth curve, widely used in curve design and graphics.

Newton-Cotes quadrature approximates $\int_a^b f(x)\,dx$ by integrating the Lagrange interpolating polynomial through $n+1$ equally spaced points $x_i=a+ih$, $h=\dfrac{b-a}{n}$.

$\int_a^b f(x)\,dx\approx\sum_{i=0}^{n}w_if(x_i)$

The weights $w_i$ depend only on $n$. The common cases are:

• $n=1$, Trapezoidal rule: $\displaystyle\int_a^b f\,dx\approx\frac h2\,(y_0+y_1)$
• $n=2$, Simpson's 1/3 rule: $\displaystyle\frac h3\,(y_0+4y_1+y_2)$
• $n=3$, Simpson's 3/8 rule: $\displaystyle\frac{3h}{8}\,(y_0+3y_1+3y_2+y_3)$

Composite versions apply the rule on many sub-intervals; e.g. composite Simpson: $\dfrac h3[y_0+y_n+4(\text{odd})+2(\text{even})]$ with $n$ even.

Notes. A rule with $n$ even is exact for polynomials of degree $n+1$; with $n$ odd, degree $n$. For large $n$ some weights become negative, so low-order composite rules are preferred.

PDEs of second order are classified as elliptic (Laplace, Poisson), parabolic (heat equation) and hyperbolic (wave equation). ODEs need only a step-by-step march in one variable, whereas PDEs need a grid in several variables and give large systems of equations.

Types of errors

1. Inherent (data) error: present in the input data because of measurement or earlier approximation.
2. Round-off error: caused by representing numbers with a finite number of digits (e.g. $\pi\approx3.1416$).
3. Truncation error: caused by cutting an infinite process short, e.g. keeping only a few terms of a Taylor series.
4. Absolute error: $E_a=|x_{true}-x_{approx}|$.
5. Relative error: $E_r=\dfrac{E_a}{|x_{true}|}$; percentage error $=E_r\times100$.
6. Blunders/Modelling errors: human mistakes or wrong mathematical models.

Root of $f(x)=2x^3+4x^2-4x-6$ by Newton-Raphson. $f(1)=-4<0$ and $f(2)=18>0$, so a root lies in $(1,2)$; since $f(1.5)=3.75>0$, take $x_0=1.5$. $f'(x)=6x^2+8x-4$.

The root, correct to three decimals, is $x\approx1.303$. (Indeed $f(x)=(x+1)(2x^2+2x-6)$, giving $x=\frac{-1+\sqrt{13}}{2}=1.3028$.)

Let $f(x)=\sin x+x-1$ (from $\sin x=1-x$), $f'(x)=\cos x+1$. Since $f(0)=-1<0$ and $f(1)=0.8415>0$, a root lies in $(0,1)$. Take $x_0=0.5$ (radians).

$x_{n+1}=x_n-\frac{\sin x_n+x_n-1}{\cos x_n+1}$

The last two iterates agree to three decimals, so the root is $x\approx0.511$.

Augmented matrix:

$\left[\begin{array}{ccc|c}2&2&1&6\\4&2&3&4\\1&1&1&0\end{array}\right]$

Pivoting. The largest entry in column 1 is 4, so interchange $R_1\leftrightarrow R_2$:

$\left[\begin{array}{ccc|c}4&2&3&4\\2&2&1&6\\1&1&1&0\end{array}\right]$

Eliminate $x$: $R_2\to R_2-\tfrac12R_1$, $R_3\to R_3-\tfrac14R_1$:

$\left[\begin{array}{ccc|c}4&2&3&4\\0&1&-\tfrac12&4\\0&\tfrac12&\tfrac14&-1\end{array}\right]$

In column 2 the pivot 1 is already the larger, so no swap is needed. $R_3\to R_3-\tfrac12R_2$:

$\left[\begin{array}{ccc|c}4&2&3&4\\0&1&-\tfrac12&4\\0&0&\tfrac12&-3\end{array}\right]$

Back substitution

$z=\frac{-3}{1/2}=-6,\qquad y=4+\tfrac12z=1,\qquad x=\frac{4-2y-3z}{4}=\frac{4-2+18}{4}=5$

Solution: $x=5,\ y=1,\ z=-6$. Check eq. 1: $10+2-6=6$.

Residuals: $R_1=-12-5x-2y-z$, $R_2=20+x-4y-2z$, $R_3=3-2x+3y-10z$. Start from $(0,0,0)$, so $R=(-12,\,20,\,3)$. At each step the largest residual is reduced to zero by changing the corresponding variable by $\delta=R_j/a_{jj}$; the other residuals are updated using the column of the matrix.

The residuals are now small, giving

$x\approx-4.0,\quad y\approx3.0,\quad z\approx2.0$

Check: $5(-4)+2(3)+2=-12$; $-(-4)+4(3)+2(2)=20$; $2(-4)-3(3)+10(2)=3$.

Write $A=LU$ (Doolittle, unit diagonal in $L$) for

$A=\begin{bmatrix}3&2&1\\2&3&2\\1&2&3\end{bmatrix}$

Factorisation. Row 1 of $U$: $(3,2,1)$. Multipliers: $l_{21}=\tfrac23$, $l_{31}=\tfrac13$.

$u_{22}=3-\tfrac23(2)=\tfrac53$, $u_{23}=2-\tfrac23(1)=\tfrac43$, $l_{32}=\dfrac{2-\tfrac13(2)}{5/3}=\tfrac45$, $u_{33}=3-\tfrac13(1)-\tfrac45\cdot\tfrac43=\tfrac85$.

$L=\begin{bmatrix}1&0&0\\\tfrac23&1&0\\\tfrac13&\tfrac45&1\end{bmatrix},\qquad U=\begin{bmatrix}3&2&1\\0&\tfrac53&\tfrac43\\0&0&\tfrac85\end{bmatrix}$

Forward substitution $Ly=b$, $b=(10,14,14)^T$:

$y_1=10$, $y_2=14-\tfrac23(10)=\tfrac{22}{3}$, $y_3=14-\tfrac{10}{3}-\tfrac45\cdot\tfrac{22}{3}=\tfrac{24}{5}$

Back substitution $Ux=y$:

$z=\frac{24/5}{8/5}=3,\qquad y=\frac{\tfrac{22}{3}-\tfrac43(3)}{5/3}=2,\qquad x=\frac{10-2(2)-3}{3}=1$

Solution: $x=1,\ y=2,\ z=3$.