Numerical Method

Poisson's equation by finite differences. With a square grid of step $h$, the five-point formula for $\nabla^2f=g(x,y)$ is

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

Here $h=1$ on $0\le x,y\le3$, so the interior points are $(1,1),(2,1),(1,2),(2,2)$ and $g=3x^2y$ gives

$g(1,1)=3,\quad g(2,1)=12,\quad g(1,2)=6,\quad g(2,2)=24.$

Boundary values are all zero. Writing $f_{1,1}=a,\;f_{2,1}=b,\;f_{1,2}=c,\;f_{2,2}=d$:

$-4a+b+c=3$ $a-4b+d=12$ $a-4c+d=6$ $b+c-4d=24$

Subtracting equation 3 from equation 2: $-4b+4c=6\Rightarrow c=b+1.5$.

Let $s=b+c$. Adding equations 2 and 3 gives $2a+2d-4s=18$, so $a+d=9+2s$. Adding equations 1 and 4 gives $-4(a+d)+2s=27$. Substituting:

$-4(9+2s)+2s=27\Rightarrow-6s=63\Rightarrow s=-10.5$

Then $2b+1.5=-10.5$, so $b=-6$, $c=-4.5$, and $a+d=9-21=-12$. From equation 1: $-4a-6-4.5=3\Rightarrow a=-3.375$, hence $d=-8.625$.

$f(1,1)=-3.375,\quad f(2,1)=-6,\quad f(1,2)=-4.5,\quad f(2,2)=-8.625$

Check eq. 4: $-6-4.5+34.5=24$. Correct.

Heun's method. For $y'=f(x,y)$, $y(x_0)=y_0$, with step $h$:

$y^{p}_{i+1}=y_i+hf(x_i,y_i),\qquad y_{i+1}=y_i+\frac h2\left[f(x_i,y_i)+f(x_{i+1},y^{p}_{i+1})\right]$

Algorithm.

1. Read $x_0$, $y_0$, step size $h$ and the final value $x_n$.
2. Compute the number of steps $n=(x_n-x_0)/h$.
3. Repeat $n$ times: $k_1=f(x,y)$; $y_p=y+hk_1$; $k_2=f(x+h,y_p)$; $y=y+\frac h2(k_1+k_2)$; $x=x+h$; print $x,y$.

Program (C). The example equation is $y'=x+y$; change f for another problem.

#include <stdio.h>
#include <math.h>

double f(double x, double y) {
    return x + y;          /* dy/dx = f(x, y) */
}

int main(void) {
    double x0, y0, xn, h, x, y, k1, k2, yp;
    int i, n;

    printf("Enter x0, y0: ");
    scanf("%lf %lf", &x0, &y0);
    printf("Enter final x and step size h: ");
    scanf("%lf %lf", &xn, &h);

    n = (int)floor((xn - x0) / h + 0.5);
    x = x0;
    y = y0;
    printf("%10s %12s\n", "x", "y");
    printf("%10.4f %12.6f\n", x, y);
    for (i = 0; i < n; i++) {
        k1 = f(x, y);
        yp = y + h * k1;               /* predictor (Euler) */
        k2 = f(x + h, yp);
        y = y + h / 2.0 * (k1 + k2);   /* corrector (trapezoidal) */
        x = x + h;
        printf("%10.4f %12.6f\n", x, y);
    }
    return 0;
}

Romberg integration improves trapezoidal estimates by Richardson extrapolation: $R_{k,1}=T(h_k),\quad h_k=\frac{b-a}{2^{k-1}},\qquad R_{k,j}=R_{k,j-1}+\frac{R_{k,j-1}-R_{k-1,j-1}}{4^{\,j-1}-1}$

Here $f(x)=\dfrac{\sin x}{x}$ with $f(0)=1$ (limit), $a=0,\ b=1$.

Function values: $f(0)=1$, $f(1)=\sin 1=0.841471$, $f(0.5)=\dfrac{0.479426}{0.5}=0.958851$, $f(0.25)=\dfrac{0.247404}{0.25}=0.989616$, $f(0.75)=\dfrac{0.681639}{0.75}=0.908852$.

Trapezoidal estimates

• $h=1$: $R_{1,1}=\dfrac12(1+0.841471)=0.920735$
• $h=0.5$: $R_{2,1}=\dfrac12R_{1,1}+0.5\,f(0.5)=0.460368+0.479426=0.939793$
• $h=0.25$: $R_{3,1}=\dfrac12R_{2,1}+0.25\,[f(0.25)+f(0.75)]=0.469897+0.25(1.898468)=0.944514$
• $h=0.125$: $R_{4,1}=0.945691$ (adding the four new midpoints $f(0.125),f(0.375),f(0.625),f(0.875)$)

Romberg table

Sample extrapolations: $R_{2,2}=0.939793+\frac{0.939793-0.920735}{3}=0.946146$ $R_{3,2}=0.944514+\frac{0.944514-0.939793}{3}=0.946087$ $R_{3,3}=0.946087+\frac{0.946087-0.946146}{15}=0.946083$

Successive entries agree to six decimals, so $\int_0^1\frac{\sin x}{x}\,dx\approx 0.946083$

Condition for convergence. Gauss-Seidel converges when the coefficient matrix is diagonally dominant (sufficient condition). In the order given, the system is not diagonally dominant, and the iteration diverges. No re-ordering makes the matrix strictly diagonally dominant, but the order below makes the iteration converge (spectral radius $\approx0.85<1$): $5x+2y+2z=25,\qquad x+2y+3z=4,\qquad 6x-4y+5z=10$

Iteration formulas (latest values are used immediately): $x^{(k+1)}=\frac{25-2y^{(k)}-2z^{(k)}}{5},\qquad y^{(k+1)}=\frac{4-x^{(k+1)}-3z^{(k)}}{2},\qquad z^{(k+1)}=\frac{10-6x^{(k+1)}+4y^{(k+1)}}{5}$

Start with $x=y=z=0$.

Example, iteration 1: $x=\dfrac{25}{5}=5$; $y=\dfrac{4-5-0}{2}=-0.5$; $z=\dfrac{10-30-2}{5}=-4.4$.

The iterates oscillate and settle slowly (the convergence is slow because the system is only weakly dominant). They approach $x\approx4.79,\quad y\approx2.38,\quad z\approx-1.85$

Check by direct elimination. Eliminating gives the exact solution $x=4.7885,\ y=2.375,\ z=-1.8462$. Substitution in the original equations: $4.7885+4.75-5.5385=4.00$, $28.731-9.5-9.231=10.00$, $23.942+4.75-3.692=25.00$. All three equations are satisfied.

So $x\approx4.79,\ y\approx2.38,\ z\approx-1.85$.

Setting up. Poisson's equation is $\nabla^2u=\dfrac{\partial^2u}{\partial x^2}+\dfrac{\partial^2u}{\partial y^2}=g(x,y)$. Using the five-point central difference formula on a square grid of side $h$: $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)$

With $0\le x,y\le1.5$ and $h=0.5$, the grid points are $x,y\in\{0,\,0.5,\,1.0,\,1.5\}$, a $4\times4$ grid. The boundary nodes (12 points) all have the prescribed value $0$. The interior nodes are $(0.5,0.5),\ (1,0.5),\ (0.5,1),\ (1,1)$, which we call $u_1,u_2,u_3,u_4$:

u3  u4
u1  u2

Equations. The source term is $g=f=0$ (as given in the question), so the equation is $\nabla^2u=0$ (Laplace form). For each interior node:

• $u_1$: $0+u_2+u_3+0-4u_1=0$
• $u_2$: $u_1+0+0+u_4-4u_2=0$
• $u_3$: $0+u_4+u_1+0-4u_3=0$
• $u_4$: $u_3+0+u_2+0-4u_4=0$

By symmetry $u_1=u_2=u_3=u_4=u$, so $2u-4u=0\Rightarrow u=0$.

Check by elimination: adding all four equations gives $-2(u_1+u_2+u_3+u_4)=0$, so the sum is $0$, and then each equation $4u_i=u_j+u_k$ with the symmetric structure gives $u_i=0$.

$u_1=u_2=u_3=u_4=0$

Remark. If a non-zero source $g$ were supplied, then by symmetry each interior node would satisfy $-2u=h^2g$, i.e. $u=-h^2g/2=-0.125\,g$ for a constant $g$. With $g=0$ and zero boundary values, the solution is identically zero.

Forward difference table ($h=10$)

Newton's forward interpolation formula $f(x)=f_0+p\,\Delta f_0+\frac{p(p-1)}{2!}\Delta^2f_0+\frac{p(p-1)(p-2)}{3!}\Delta^3f_0+\frac{p(p-1)(p-2)(p-3)}{4!}\Delta^4f_0,\qquad p=\frac{x-x_0}{h}$

For $x=15$: $p=\dfrac{15-10}{10}=0.5$.

Terms:

• $f_0=0.173$
• $p\,\Delta f_0=0.5(0.169)=0.0845$
• $\dfrac{0.5(-0.5)}{2}(-0.011)=(-0.125)(-0.011)=0.001375$
• $\dfrac{0.5(-0.5)(-1.5)}{6}(-0.004)=(0.0625)(-0.004)=-0.00025$
• $\dfrac{0.5(-0.5)(-1.5)(-2.5)}{24}(-0.001)=(-0.0390625)(-0.001)=0.000039$

$f(15)=0.173+0.0845+0.001375-0.00025+0.000039=0.258664$ $\boxed{f(15)\approx0.2587}$

(The data are $\sin x^\circ$; the true value of $\sin15^\circ$ is $0.2588$.)

Romberg integration. Romberg integration improves trapezoidal estimates by Richardson extrapolation. With $T_k$ the trapezoidal value using $2^{k}$ strips:

$R_{k,2}=\frac{4T_{k+1}-T_k}{3},\quad R_{k,3}=\frac{16R_{k+1,2}-R_{k,2}}{15},\quad R_{k,4}=\frac{64R_{k+1,3}-R_{k,3}}{63}$

Problem. $f(x)=\dfrac{\sin^2x}{x}$, $f(0)=0$ (limit), $f(1)=\sin^21=0.708073$.

• $h=1$ (1 strip): $T_1=\tfrac12(0+0.708073)=0.354037$
• $h=0.5$: $f(0.5)=0.459698$, so $T_2=\tfrac12T_1+0.5(0.459698)=0.406867$
• $h=0.25$: $f(0.25)=0.244835$, $f(0.75)=0.619509$, so $T_3=\tfrac12T_2+0.25(0.864344)=0.419519$
• $h=0.125$: $f(0.125)=0.124350$, $f(0.375)=0.357748$, $f(0.625)=0.547742$, $f(0.875)=0.673283$ (sum $1.703124$), so $T_4=\tfrac12T_3+0.125(1.703124)=0.422650$

Romberg table:

For example $\dfrac{4(0.406867)-0.354037}{3}=0.424477$ and $\dfrac{16(0.423737)-0.424477}{15}=0.423688$.

$\int_0^1\frac{\sin^2x}{x}\,dx\approx0.42369$

Gauss-Jordan method. The augmented matrix is reduced to the identity (reduced row echelon form) by row operations, so the solution can be read directly from the last column.

Augmented matrix (equations written with all four unknowns):

$\left[\begin{array}{cccc|c}0&1&2&3&9\\7&6&5&4&33\\8&9&0&1&27\\2&5&4&3&23\end{array}\right]$

The first pivot is zero, so interchange $R_1\leftrightarrow R_3$ and normalise ($R_1\leftarrow R_1/8$), then eliminate column 1:

$\left[\begin{array}{cccc|c}1&\frac98&0&\frac18&\frac{27}8\\0&-\frac{15}8&5&\frac{25}8&\frac{75}8\\0&1&2&3&9\\0&\frac{11}4&4&\frac{11}4&\frac{65}4\end{array}\right]$

Normalise $R_2$ by its pivot $-\tfrac{15}{8}$ and clear column 2 in all other rows:

$\left[\begin{array}{cccc|c}1&0&3&2&9\\0&1&-\frac83&-\frac53&-5\\0&0&\frac{14}3&\frac{14}3&14\\0&0&\frac{34}3&\frac{22}3&30\end{array}\right]$

Third stage (normalise $R_3$: $R_3\leftarrow\frac3{14}R_3$ gives $(0,0,1,1,3)$, then clear column 3):

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

Fourth stage ($R_4\leftarrow R_4/(-4)$, then clear column 4):

$\left[\begin{array}{cccc|c}1&0&0&0&1\\0&1&0&0&2\\0&0&1&0&2\\0&0&0&1&1\end{array}\right]$

$x_1=1,\quad x_2=2,\quad x_3=2,\quad x_4=1$

Check in equation 2: $7+12+10+4=33$. Correct.

Gauss-Seidel method. Each unknown is solved from its own equation, always using the newest available values:

$x^{(k+1)}=\frac{\ldots}{a_{11}},\quad y^{(k+1)}\text{ uses }x^{(k+1)},\quad z^{(k+1)}\text{ uses }x^{(k+1)},y^{(k+1)}$

The method converges if the system is diagonally dominant (or if the iteration matrix has spectral radius below 1). The given equations are not diagonally dominant in any order, but taking them in the order

$5x+2y+3z=12,\qquad 6x+4y+5z=16,\qquad x+2y+3z=4$

gives a convergent iteration (spectral radius $5/6<1$):

$x=\frac{12-2y-3z}{5},\qquad y=\frac{16-6x-5z}{4},\qquad z=\frac{4-x-2y}{3}$

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

The values tend to $x=2,\;y=1,\;z=0$ (convergence in $y,z$ is slow). Verification: $2+2+0=4$, $12+4+0=16$, $10+2+0=12$.

$x=2,\quad y=1,\quad z=0$

Backward difference table ($h=10$, $\nabla f_i=f_i-f_{i-1}$)

Newton's backward formula

$f(x)=f_n+p\nabla f_n+\frac{p(p+1)}{2!}\nabla^2f_n+\frac{p(p+1)(p+2)}{3!}\nabla^3f_n+\frac{p(p+1)(p+2)(p+3)}{4!}\nabla^4f_n,\qquad p=\frac{x-x_n}{h}$

Here $x_n=50$, so $p=\dfrac{45-50}{10}=-0.5$.

• $p\nabla f_n=(-0.5)(-0.123)=0.0615$
• $\dfrac{p(p+1)}{2}\nabla^2f_n=\dfrac{(-0.5)(0.5)}{2}(-0.023)=(-0.125)(-0.023)=0.002875$
• $\dfrac{p(p+1)(p+2)}{6}\nabla^3f_n=\dfrac{(-0.5)(0.5)(1.5)}{6}(0.009)=(-0.0625)(0.009)=-0.0005625$
• $\dfrac{p(p+1)(p+2)(p+3)}{24}\nabla^4f_n=\dfrac{(-0.5)(0.5)(1.5)(2.5)}{24}(0.024)=(-0.0390625)(0.024)=-0.0009375$

$f(45)=0.643+0.0615+0.002875-0.0005625-0.0009375=0.705875\approx0.7059$

Finite-difference form. With equal step $h$ in both directions, the five-point formula for $\nabla^2f=g(x,y)$ is

$f_{i+1,j}+f_{i-1,j}+f_{i,j+1}+f_{i,j-1}-4f_{i,j}=h^2g(x_i,y_j),\qquad g=xy,\ h=1$

The grid has $x,y\in\{0,1,2,3\}$. The interior points are $(1,1),(2,1),(1,2),(2,2)$; all boundary values are $2$. Let $f_1=f(1,1)$, $f_2=f(2,1)$, $f_3=f(1,2)$, $f_4=f(2,2)$.

• At $(1,1)$, $g=1$: $f_2+f_3+2+2-4f_1=1\Rightarrow-4f_1+f_2+f_3=-3$
• At $(2,1)$, $g=2$: $f_1+f_4+2+2-4f_2=2\Rightarrow f_1-4f_2+f_4=-2$
• At $(1,2)$, $g=2$: $f_1+f_4+2+2-4f_3=2\Rightarrow f_1-4f_3+f_4=-2$
• At $(2,2)$, $g=4$: $f_2+f_3+2+2-4f_4=4\Rightarrow f_2+f_3-4f_4=0$

By symmetry $f_2=f_3$:

$-4f_1+2f_2=-3,\qquad f_1-4f_2+f_4=-2,\qquad 2f_2-4f_4=0\Rightarrow f_4=\tfrac12f_2$

Second equation: $f_1-4f_2+0.5f_2=-2\Rightarrow f_1=3.5f_2-2$. Substituting in the first: $-14f_2+8+2f_2=-3\Rightarrow f_2=\tfrac{11}{12}$. Then $f_1=\tfrac{77}{24}-2=\tfrac{29}{24}$, $f_4=\tfrac{11}{24}$.

$f(1,1)=1.2083,\quad f(2,1)=f(1,2)=0.9167,\quad f(2,2)=0.4583$

Augmented matrix:

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

$R_2\to R_2-2R_1$, $R_3\to R_3-R_1$:

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

Interchange $R_2\leftrightarrow R_3$ (to bring a non-zero pivot up), then $R_2\to R_2/(-4)$:

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

$R_1\to R_1-2R_2$:

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

The last row is $0=0$, so there are only two independent equations for three unknowns: the second equation of the system is twice the first. The system is consistent with infinitely many solutions. Taking $z=t$ as a free parameter:

$x=4-t,\qquad y=2t,\qquad z=t\qquad(t\in\mathbb R)$

Check with $t=1$: $(3,2,1)$ gives $3+4-3=4$, $6+8-6=8$, $3-4+5=4$. All three equations hold.

Non-linear equation. An equation $f(x)=0$ in which $f$ is not a linear function of the unknown (it contains terms such as $x^2$, $\sin x$, $e^{x}$, $\log x$) is called a non-linear equation. Such equations usually cannot be solved in closed form, so a root is found by an iterative numerical method.

Derivation of the secant method. Let $x_{n-1}$ and $x_n$ be two approximations of the root. Draw the secant line through $(x_{n-1}, f(x_{n-1}))$ and $(x_n, f(x_n))$ and take the point where it cuts the $x$-axis as the next approximation $x_{n+1}$.

The slope of the secant is $m=\frac{f(x_n)-f(x_{n-1})}{x_n-x_{n-1}}$

The line through $(x_n, f(x_n))$ with this slope is $y-f(x_n)=m(x-x_n)$. Putting $y=0$, $x=x_{n+1}$: $x_{n+1}=x_n-\frac{f(x_n)}{m}$ $\boxed{x_{n+1}=x_n-\frac{f(x_n)\,(x_n-x_{n-1})}{f(x_n)-f(x_{n-1})}}$

This is Newton's formula with $f'(x_n)$ replaced by the finite-difference slope, so no derivative is needed. The order of convergence is about $1.618$.

Solving $f(x)=x^2+\cos x-e^{-x}-2=0$ (angles in radians).

$f(1)=1+0.5403-0.3679-2=-0.8276$, $f(2)=4-0.4161-0.1353-2=1.4485$, so a root lies in $(1,2)$. Take $x_0=1$, $x_1=2$.

Example of one step: $x_2=2-\dfrac{1.44852\,(2-1)}{1.44852-(-0.82758)}=2-0.63640=1.36360$.

The successive values agree to five decimals, so the root is $x\approx 1.4554$

Matrix factorization. Matrix factorization (decomposition) writes a matrix $A$ as a product of simpler matrices. In LU factorization, $A=LU$ where $L$ is lower triangular and $U$ is upper triangular. Then $Ax=b$ is solved in two easy triangular steps: $Ly=b$ (forward substitution) and $Ux=y$ (back substitution).

Doolittle algorithm. $L$ has a unit diagonal ($l_{ii}=1$). For $i=1,\dots,n$: $u_{ij}=a_{ij}-\sum_{k=1}^{i-1}l_{ik}u_{kj},\quad j=i,\dots,n$ $l_{ji}=\frac{1}{u_{ii}}\left(a_{ji}-\sum_{k=1}^{i-1}l_{jk}u_{ki}\right),\quad j=i+1,\dots,n$

Factorization of A.

Row 1 of $U$: $u_{11}=2,\ u_{12}=4,\ u_{13}=-4,\ u_{14}=0$.

Column 1 of $L$: $l_{21}=1/2=0.5,\ l_{31}=2/2=1,\ l_{41}=1/2=0.5$.

Row 2 of $U$: $u_{22}=5-0.5(4)=3$, $u_{23}=-5-0.5(-4)=-3$, $u_{24}=-3-0.5(0)=-3$.

Column 2 of $L$: $l_{32}=\dfrac{3-1(4)}{3}=-\dfrac13$, $l_{42}=\dfrac{4-0.5(4)}{3}=\dfrac23$.

Row 3 of $U$: $u_{33}=1-1(-4)-(-\tfrac13)(-3)=1+4-1=4$, $u_{34}=3-1(0)-(-\tfrac13)(-3)=3-1=2$.

Column 3 of $L$: $l_{43}=\dfrac{-2-0.5(-4)-\tfrac23(-3)}{4}=\dfrac{-2+2+2}{4}=\dfrac12$.

Row 4 of $U$: $u_{44}=2-0.5(0)-\tfrac23(-3)-\tfrac12(2)=2+2-1=3$.

$L=\begin{bmatrix}1&0&0&0\\0.5&1&0&0\\1&-\frac13&1&0\\0.5&\frac23&\frac12&1\end{bmatrix},\qquad U=\begin{bmatrix}2&4&-4&0\\0&3&-3&-3\\0&0&4&2\\0&0&0&3\end{bmatrix}$

Solving Ax = b.

Forward substitution $Ly=b$:

• $y_1=12$
• $y_2=18-0.5(12)=12$
• $y_3=8-1(12)+\tfrac13(12)=8-12+4=0$
• $y_4=8-0.5(12)-\tfrac23(12)-\tfrac12(0)=8-6-8=-6$

Back substitution $Ux=y$:

• $3x_4=-6\Rightarrow x_4=-2$
• $4x_3+2x_4=0\Rightarrow x_3=1$
• $3x_2-3x_3-3x_4=12\Rightarrow 3x_2=12+3-6=9\Rightarrow x_2=3$
• $2x_1+4x_2-4x_3=12\Rightarrow 2x_1=12-12+4=4\Rightarrow x_1=2$

$x=(2,\ 3,\ 1,\ -2)^T$

Check in the first equation: $2(2)+4(3)-4(1)+0=12$. In the last: $2+12-2-4=8$. Both are correct.