1 Three point stencil

Let us consider the boundary value problem defined by Equations (1):

$\left.\begin{array}{ccc}\hfill p\left(x\right)y^{″}\left(x\right)+q\left(x\right)y^{\prime }\left(x\right)+r\left(x\right)y\left(x\right)& \hfill =\hfill & s\left(x\right)\left(0<x<\mathrm{\&ell;}\right)\hfill \\ \hfill a_0{y}^{\prime }\left(0\right)+b_{0}y\left(0\right)& \hfill =\hfill & c_0\hfill \\ \hfill a_1{y}^{\prime }\left(\mathrm{\&ell;}\right)+b_{1}y\left(\mathrm{\&ell;}\right)& \hfill =\hfill & c_1\hfill \end{array}\right\}\left(1\right)$

The first line is the differential equation , and the second and third lines are the boundary conditions which can involve derivatives.

It is our aim to approximate the solution of this problem numerically, and we adopt an approach similar to that seen in HELM booklet  32.

We divide the interval $0<x<\mathrm{\&ell;}$ into a number, $J$ say, of subintervals each of equal width $h=\mathrm{\&ell;}∕J$ . Our numerical solution will provide an approximation to $y=y\left(x\right)$ at each value of $x$ where two subintervals meet (see Figure 1).

Figure 1

A numerical approximation to Equations (1) can be found by approximating the derivatives by finite differences. Here we will approximate $y^{″}\left(x\right)$ and $y^{\prime }\left(x\right)$ by central differences to obtain

$p\left(x_j\right)\frac{y\left(x_j+h\right)-2y\left(x_j\right)+y\left(x_j-h\right)}{h^2}+q\left(x_j\right)\frac{y\left(x_j+h\right)-y\left(x_j-h\right)}{2h}+r\left(x_j\right)y\left(x_j\right)\approx s\left(x_j\right)$

that is, on using the abbreviations established in Key Point 1,

$p_j\frac{y\left(x_j+h\right)-2y\left(x_j\right)+y\left(x_j-h\right)}{h^2}+q_j\frac{y\left(x_j+h\right)-y\left(x_j-h\right)}{2h}+r_{j}y\left(x_j\right)\approx s_j$

which we use as the motivation for the numerical method

$p_j\frac{y_{j+1}-2y_j+y_{j-1}}{h^2}+q_j\frac{y_{j+1}-y_{j-1}}{2h}+r_j{y}_j=s_j.$

This last equation can be rearranged, gathering together all the like $y$ -terms. It neatens things further to multiply through by $h^2$ as well, and the result of these manipulations appears in the following Key Point.

This approximation to the differential equation can be thought of as a three-point stencil linking three of the approximate $y$ -values. The expression

$y_{j-1}\left(p_j-\frac{h}{2}{q}_j\right)+y_j\left(h^2{r}_j-2p_j\right)+y_{j+1}\left(p_j+\frac{h}{2}{q}_j\right)=h^2{s}_j.$

is centred around $x=x_j$ and involves $y_{j-1}$ , $y_j$ and $y_{j+1}$ . The general rule when dealing with a numerical stencil like this is to centre the stencil at every point where $y$ is unknown . (This general rule will appear again, for example on page 13.)

In the following Example, matters are simplified because the functions $p$ , $q$ , $r$ and $s$ are all constant.

Example 1

Let $y=y\left(x\right)$ be a solution to the boundary value problem

$3y^{″}\left(x\right)+4y^{\prime }\left(x\right)+5y\left(x\right)=70<x<1.5$

$y\left(0\right)=2$ , $y\left(1.5\right)=2$ .

Using a mesh width of $h=0.5,$ obtain a central difference approximation to the differential equation and hence find $y_j\approx y\left(jh\right)$ , $j=1,2$ .

Solution

In general, the central difference approximation to

$p\left(x\right)y^{″}\left(x\right)+q\left(x\right)y^{\prime }\left(x\right)+r\left(x\right)y\left(x\right)=s\left(x\right)$

is

$y_{j-1}\left(p_j-\frac{h}{2}{q}_j\right)+y_j\left(h^2{r}_j-2p_j\right)+y_{j+1}\left(p_j+\frac{h}{2}{q}_j\right)=h^2{s}_j.$

In this case the coefficients are

$p_j-\frac{h}{2}{q}_j=2,h^2{r}_j-2p_j=-4.75,p_j+\frac{h}{2}{q}_j=4,h^2{s}_j=1.75.$

These values will be the same for all $x$ because $p$ , $q$ , $r$ and $s$ are constants in this Example. Hence the general stencil is

$2y_{j-1}-4.75y_j+4y_{j+1}=1.75$

In this case $h=0.5,$ and our numerical solution consists of the values

So there are two unknowns, $y_1$ and $y_2$ . We centre the stencil at each of the corresponding $x$ values. Putting $j=1$ in the numerical stencil gives

$2y_0-4.75y_1+4y_2=1.75$

Moving the stencil one place to the right, we put $j=2$ so that

$2y_1-4.75y_2+4y_3=1.75$

In these two equations $y_0$ and $y_3$ are known from the boundary conditions and we move terms involving them to the right-hand side. This leads to the system of equations

$\left[\begin{array}{cc}\hfill -4.75\hfill & \hfill 4\hfill \\ \hfill 2\hfill & \hfill -4.75\hfill \end{array}\right]\left[\begin{array}{c}\hfill y_1\hfill \\ \hfill y_2\hfill \end{array}\right]=\left[\begin{array}{c}\hfill -2.25\hfill \\ \hfill -6.25\hfill \end{array}\right]$

Solving this pair of simultaneous equations we find that

$y_1=2.45,y_2=2.35$

to 2 decimal places.

This approximation is shown in Figure 2 in which the numerical approximations to point values of $y$ are shown as circles.

The question remains how close to the exact solution these approximations are. (Of course for Example 1 above it is possible to find the analytic solution fairly easily, but this will not usually be the case.)

A pragmatic way to deal with this question is to recompute the results with a smaller value of $h$ . We know from HELM booklet  31 that the central difference approximations get closer and closer to the derivatives which they approximate as $h$ decreases. In Figure 3 the results for $h=0.5$ are given again as circles, and a computer has been used to find more accurate approximations to $y$ using $h=\frac{1.5}{7}$ (shown as squares) and yet more accurate results (shown as dots) from using $h=\frac{1.5}{10}$ . (This involves solving larger systems of equations than are manageable by hand. The methods seen in HELM booklet  30 can be used to deal with these larger systems.)

We can now have some confidence that the results we calculated using $h=0.5$ tended to overestimate the true values of $y$ .

Example 2

Let $y=y\left(x\right)$ be a solution to the boundary value problem

$y^{″}\left(x\right)+2y^{\prime }\left(x\right)-2y\left(x\right)=-30<x<2$

$y\left(0\right)=1$ , $y\left(2\right)=-2$ .

Using a mesh width of $h=0.5$ obtain a central difference approximation to the differential equation and hence find simultaneous equations satisfied by the unknowns $y_1$ , $y_2$ and $y_3$ .

Solution

In this case the coefficients are

$p_j-\frac{h}{2}{q}_j=0.5,h^2{r}_j-2p_j=-2.5,p_j+\frac{h}{2}{q}_j=1.5,h^2{s}_j=-0.75.$

These values will be the same for all $x$ because $p$ , $q$ , $r$ and $s$ are constants in this Example.

Hence the general stencil is

$0.5y_{j-1}-2.5y_j+1.5y_{j+1}=-3$

Here we have $h=0.5$ and our numerical solution consists of the values

So there are three unknowns, $y_1$ , $y_2$ and $y_3$ . We centre the stencil at each of the corresponding $x$ values. Putting $j=1$ in the numerical stencil gives

$0.5y_0-2.5y_1+1.5y_2=-0.75.$

Moving the stencil one place to the right, we put $j=2$ so that

$0.5y_1-2.5y_2+1.5y_3=-0.75$

and finally we let $j=3$ so that

$0.5y_2-2.5y_3+1.5y_4=-0.75$

In these three equations $y_0$ and $y_4$ are known from the boundary conditions and we move terms involving them to the right-hand side. This leads to the system of equations

$\left[\begin{array}{ccc}\hfill -2.5\hfill & \hfill 1.5\hfill & \hfill 0\hfill \\ \hfill 0.5\hfill & \hfill -2.5\hfill & \hfill 1.5\hfill \\ \hfill 0\hfill & \hfill 0.5\hfill & \hfill -2.5\hfill \end{array}\right]\left[\begin{array}{c}\hfill y_1\hfill \\ \hfill y_2\hfill \\ \hfill y_3\hfill \end{array}\right]=\left[\begin{array}{c}\hfill -1.25\hfill \\ \hfill -0.75\hfill \\ \hfill 2.25\hfill \end{array}\right]$

We can find (using methods from HELM booklet  30, for example) that the solution to the system of equations in Example 2 is

$y_1=0.39y_2=-0.18y_3=-0.94$ to 2 decimal places.

Task!

Let $y=y\left(x\right)$ be a solution to the boundary value problem

$y^{″}\left(x\right)+4y^{\prime }\left(x\right)=40<x<1$

$y\left(0\right)=-2$ , $y\left(1\right)=3$ .

Using a mesh width of $h=0.25$ obtain a central difference approximation to the differential equation and hence find a system of equations satisfied by $y_j\approx y\left(jh\right)$ , $j=1,2,3$ .

Answer

The temperature $y$ of an electrically heated wire of length $\mathrm{\&ell;}$ is affected by local air currents. This situation may be modelled by

$\frac{d^{2}y}{d{x}^2}=-a-b\left(Y-y\right),\left(0<x<\mathrm{\&ell;}\right).$

Consider the case where $\mathrm{\&ell;}=3$ , $a=50$ , $b=0.1$ and $Y=20^{\&compfn;}$ C and suppose that the ends of the wire are known (to 1 decimal place) to be at temperatures $y\left(0\right)=15.0^{\&compfn;}$ C and $y\left(3\right)=25.0^{\&compfn;}$ C.

Using a central difference to approximate the derivative and using 3 subintervals obtain approximations to the temperature $\frac{1}{3}$ and $\frac{2}{3}$ of the length along the wire.

Solution

This Example falls into the general case given at the beginning of this Section if we choose $p=1$ , $q=0$ , $r=-0.1$ and $s=-52$ . In this case $h=1$ and our numerical solution consists of the values

So there are two unknowns, $y_1$ and $y_2$ . We centre the stencil at each of the corresponding $x$ values. Putting $j=1$ in the numerical stencil gives

$y_0-2.1y_1+y_2=-52$

Moving the stencil one place to the right, we put $j=2$ so that

$y_1-2.1y_2+y_3=-52$

In these two equations $y_0$ and $y_3$ are known from the boundary conditions and we move terms involving them to the right-hand side. This leads to the system of equations

$\left[\begin{array}{cc}\hfill -2.1\hfill & \hfill 1\hfill \\ \hfill 1\hfill & \hfill -2.1\hfill \end{array}\right]\left[\begin{array}{c}\hfill y_1\hfill \\ \hfill y_2\hfill \end{array}\right]=\left[\begin{array}{c}\hfill -67\hfill \\ \hfill -77\hfill \end{array}\right]$

Solving this pair of simultaneous equations we find that

$y_1=63.8,y_2=67.1$

to 1 decimal place.

We conclude that the temperature $\frac{1}{3}$ of the wire’s length from the cooler end is approximately $63.8^{\&compfn;}$ C and the temperature the same distance from the hotter end is approximately $67.1^{\&compfn;}$ C, where we have rounded these numbers to the same number of places as the given boundary conditions.

The Examples and Task above were such that $p$ , $q$ , $r$ and $s$ were each equal to a constant for all values of $x$ . More realistic engineering applications may involve coefficients that vary, and the next Example is of this type.

Example 4

Let $y=y\left(x\right)$ be a solution to the boundary value problem

$ln\left(2+x\right)y^{″}\left(x\right)+x{y}^{\prime }\left(x\right)+{\left(x+1\right)}^{2}y\left(x\right)=cos\left(x\right)0<x<1.2$

$y\left(0\right)=0$ , $y\left(1.2\right)=2$ .

Using a mesh width of $h=0.4$ obtain a central difference approximation to the differential equation and hence find $y_j\approx y\left(jh\right)$ , $j=1,2$ .

Solution

In general, the central difference approximation to

$p\left(x\right)y^{″}\left(x\right)+q\left(x\right)y^{\prime }\left(x\right)+r\left(x\right)y\left(x\right)=s\left(x\right)$

is

$y_{j-1}\left(p_j-\frac{h}{2}{q}_j\right)+y_j\left(h^2{r}_j-2p_j\right)+y_{j+1}\left(p_j+\frac{h}{2}{q}_j\right)=h^2{s}_j.$

The coefficients will vary with $j$ in this Example because the functions $p$ , $q$ , $r$ and $s$ are not all constants. In this case $h=0.4$ and our numerical solution consists of the values

So there are two unknowns, $y_1$ and $y_2$ . We centre the stencil at each of the corresponding $x$ values. Putting $j=1$ in the numerical stencil gives

$0.795469y_0-1.437337y_1+0.955469y_2=0.147370$

Moving the stencil one place to the right, we put $j=2$ so that

$0.869619y_1-1.540839y_2+1.189619y_3=0.111473$

In these two equations $y_0$ and $y_3$ are known from the boundary conditions and we move terms involving them to the right-hand side. This gives the pair of equations

$\left[\begin{array}{cc}\hfill -1.437337\hfill & \hfill 0.955469\hfill \\ \hfill 0.869619\hfill & \hfill -1.540839\hfill \end{array}\right]\left[\begin{array}{c}\hfill y_1\hfill \\ \hfill y_2\hfill \end{array}\right]=\left[\begin{array}{c}\hfill 0.147370\hfill \\ \hfill -2.267766\hfill \end{array}\right]$

Solving this pair of simultaneous equations we find that

$y_1=1.40,y_2=2.26$

to 2 decimal places.

Once again we can monitor the accuracy of the results obtained in the Example above by recomputing for a smaller value of $h$ . In Figure 4 the values calculated are shown as circles and a computer has been used to obtain the more accurate results (shown as dots) obtained from choosing $h=\frac{1.2}{20}$ .

In the next Example we see that the derivative $y^{\prime }$ appears in the boundary condition at $x=0$ . This means that $y$ is not given at $x=0$ and we use the general rule given earlier in Key Point 3:

So this implies that we must centre the stencil at $x=0$ and this will cause the value $y_{-1}$ to appear. This is a fictitious value that plays no part in the solution we seek and we use the derivative boundary condition to get $y_{-1}$ in terms of $y_1$ . This is done with the central difference

$y^{\prime }\left(0\right)\approx \frac{y_1-y_{-1}}{2h}.$

The following Example implements this idea.

Example 5

Let $y=y\left(x\right)$ be a solution to the boundary value problem

$ln\left(2+x\right)y^{″}\left(x\right)+x{y}^{\prime }\left(x\right)+2y\left(x\right)=cos\left(x\right)0<x<1.2$

$y^{\prime }\left(0\right)=-1$ , $y\left(1.2\right)=2$

(Note the derivative boundary condition at $x=0$ .)

Using a mesh width of $h=0.4$ obtain a central difference approximation to the differential equation and hence find the system of equations satisfied by $y_j\approx y\left(jh\right)$ , $j=0,1,2$ .

Solution

In general, the central difference approximation to

$p\left(x\right)y^{″}\left(x\right)+q\left(x\right)y^{\prime }\left(x\right)+r\left(x\right)y\left(x\right)=s\left(x\right)$

is

$y_{j-1}\left(p_j-\frac{h}{2}{q}_j\right)+y_j\left(h^2{r}_j-2p_j\right)+y_{j+1}\left(p_j+\frac{h}{2}{q}_j\right)=h^2{s}_j.$

The coefficients vary with $j$ in this Example because the functions $p$ , $q$ , $r$ and $s$ are not all constants. In this case $h=0.4$ and our numerical solution consists of the values

So there are three unknowns, $y_0$ , $y_1$ and $y_2$ . We centre the stencil at each of the corresponding $x$ values. Putting $j=0$ in the numerical stencil gives

$0.693147y_{-1}-1.066294y_0+0.6931472y_1=0.16$

which introduces the fictitious quantity $y_{-1}$ . We attach a meaning to $y_{-1}$ on using the boundary condition at $x=0$ . Approximating the derivative in the boundary condition by a central difference gives

$\frac{y_1-y_{-1}}{2h}=-1\Rightarrow y_{-1}=y_1+0.8$

and we use this to remove $y_{-1}$ from the equation where it first appeared. Hence

$-1.066294y_0+1.386294y_1=0.7145177$

The remaining steps are similar to previous Examples. Putting $j=1$ in the stencil gives

$0.795469y_0-1.430937y_1+0.955469y_2=0.147370.$

Moving the stencil one place to the right, we put $j=2$ so that

$0.869619y_1-1.739239y_2+1.189619y_3=0.111473$

In the last equation $y_3$ is known from the boundary conditions and we move the term involving it to the right-hand side. This leaves us with the system of equations

$\left[\begin{array}{ccc}\hfill -1.066294\hfill & \hfill 1.386294\hfill & \hfill 0\hfill \\ \hfill 0.795469\hfill & \hfill -1.430937\hfill & \hfill 0.955469\hfill \\ \hfill 0\hfill & \hfill 0.869619\hfill & \hfill -1.739239\hfill \end{array}\right]\left[\begin{array}{c}\hfill y_0\hfill \\ \hfill y_1\hfill \\ \hfill y_2\hfill \end{array}\right]=\left[\begin{array}{c}\hfill 0.714518\hfill \\ \hfill 0.147370\hfill \\ \hfill -2.267766\hfill \end{array}\right]$

where the components are given to 6 decimal places.

Exercises

1. Let $y=y\left(x\right)$ be a solution to the boundary value problem

   $y^{″}\left(x\right)-2y^{\prime }\left(x\right)+3y\left(x\right)=60<x<0.75$

   $y\left(0\right)=2$ , $y\left(0.75\right)=1$

   Using a mesh width of $h=0.25$ obtain a central difference approximation to the differential equation and hence find $y_1\approx y\left(0.25\right)$ and $y_2\approx y\left(0.5\right)$ .

2. Let $y=y\left(x\right)$ be a solution to the boundary value problem

   $2y^{″}\left(x\right)+3y^{\prime }\left(x\right)=50<x<1.2$

   $y\left(0\right)=-2$ , $y\left(1.2\right)=3$

   Using a mesh width of $h=0.3$ obtain a central difference approximation to the differential equation and hence find a system of equations satisfied by $y_1\approx y\left(0.3\right)$ , $y_2\approx y\left(0.6\right)$ and $y_3\approx y\left(0.9\right)$ .

3. Let $y=y\left(x\right)$ be a solution to the boundary value problem

   $y^{″}\left(x\right)+x^2{y}^{\prime }\left(x\right)+3y\left(x\right)=x0<x<1.5$

   $y^{\prime }\left(0\right)=2$ , $y\left(1.5\right)=1$ . (Note the derivative boundary condition at $x=0$ .)

   Using a mesh width of $h=0.5$ obtain a central difference approximation to the differential equation and hence find the system of equations satisfied by $y_j\approx y\left(jh\right)$ , $j=0,1,2$ .

Answer