4 An implicit method

Another approach that can be used to address the initial value problem

$\frac{dy}{dt}=f\left(t,y\right),y\left(0\right)=y_0$

is to consider integrating the differential equation

$\frac{dy}{dt}=f\left(t,y\right)$

from $t=nh$ to $t=nh+h$ . This leads to

${\left[y\left(t\right)\right]}_{t=nh}^{t=nh+h}={\int \; <!--nolimits\; -->}_{nh}^{\left(n+1\right)h}f\left(t,y\right)dt$

that is,

$y\left(nh+h\right)-y\left(nh\right)={\int \; <!--nolimits\; -->}_{nh}^{\left(n+1\right)h}f\left(t,y\right)dt$

and the problem now becomes one of approximating the integral on the right-hand side.

If we approximate the integral using the simple trapezium rule and replace the terms by their approximations we obtain the numerical method

$y_{n+1}-y_n=\frac{1}{2}h\left(f_n+f_{n+1}\right)$

The procedure for time stepping with this method is much the same as that used for Euler’s method, but with one difference. Let us imagine applying the method, we are given $y_0$ as the initial condition and now aim to find $y_1$ from

And here is the problem: the unknown $y_1$ appears on both sides of the equation. We cannot, in general, find an explicit expression for $y_1$ and for this reason the numerical method is called an implicit method.

In practice the particular form of $f$ may allow us to find $y_1$ fairly simply, but in general we have to approximate $y_1$ for example by using the bisection method, or Newton-Raphson. (Another approach that can be used involves what is called a predictor-corrector method, in other words, a “guess and improve" method, and we will discuss this again later in this Workbook.)

And then, of course, we encounter the problem again in the second time step, when calculating $y_2$ . And again for $y_3$ and so on. There is, in general, a genuine cost in implementing implicit methods, but they are popular because they have desirable properties, as we will see later in this Workbook.

In Example 2 the implicit nature of the method is not a problem because $y$ does not appear on the right-hand side of the differential equation. In other words, $f=f\left(t\right)$ .

Example 2

Suppose that $y=y\left(t\right)$ is the solution to the initial value problem

$\frac{dy}{dt}=1∕\left(t+1\right),y\left(0\right)=1$

Carry out two time steps of the trapezium method with a step size of $h=0.2$ so as to obtain approximations to $y\left(0.2\right)$ and $y\left(0.4\right)$ .

Solution

For the first time step we require $f_0=f\left(0\right)=1$ and $f_1=f\left(0.2\right)=0.833333$ and therefore

$y_1=y_0+\frac{1}{2}h\left(f_0+f_1\right)=1+0.1\times 1.833333=1.183333$

For the second time step we also require $f_2=f\left(2h\right)=f\left(0.4\right)=0.714286$ and therefore

$y_2=y_1+\frac{1}{2}h\left(f_1+f_2\right)=1.183333+0.1\times 1.547619=1.338095$

We conclude that

$y\left(0.1\right)\approx 1.183333y\left(0.2\right)\approx 1.338095$

Example 3 has $f$ dependent on $y$ , so the implicit nature of the trapezium method could be a problem. However in this case the way in which $f$ depends on $y$ is simple enough for us to be able to rearrange for an explicit expression for $y_{n+1}$ .

Example 3

Suppose that $y=y\left(t\right)$ is the solution to the initial value problem

$\frac{dy}{dt}=1∕\left(t^2+1\right)-2y,y\left(0\right)=2$

Carry out two time steps of the trapezium method with a step size of $h=0.1$ so as to obtain approximations to $y\left(0.1\right)$ and $y\left(0.2\right)$ .

Solution

The trapezium method is $y_{n+1}=y_n+\frac{h}{2}\left(f_n+f_{n+1}\right)$ and in this case $y_{n+1}$ will appear on both sides because $f$ depends on $y$ . We have

where $g\left(t\right)\equiv \frac{1}{\left(t^2+1\right)}$ which is the part of $f$ that depends on $t$ . On rearranging to get all $y_{n+1}$ terms on the left, we get

$\left(1+h\right)y_{n+1}=y_n+\frac{1}{2}h\left\{g\left(t_n\right)-2y_n+g\left(t_{n+1}\right)\right\}$

In this case $h=0.1$ .

For the first time step we require $g\left(0\right)=1$ and $g\left(0.1\right)=0.990099$ and therefore

$1.1y_1=2+0.05\left(1-2\times 2+0.990099\right)$

Hence $y_1=1.726823$ , to six decimal places.

For the second time step we also require $g\left(2h\right)=g\left(0.2\right)=0.961538$ and therefore

$1.1y_2=1.726823+0.05\left(0.990099-2\times 1.726823+0.961538\right)$

Hence $y_2=1.501566$ . We conclude that $y\left(0.1\right)\approx 1.726823$ and $y\left(0.2\right)\approx 1.501566$ to 6 d.p.

Task!

Suppose that $y=y\left(t\right)$ is the solution to the initial value problem

$\frac{dy}{dt}=t-y,y\left(0\right)=2$

Carry out two time steps of the trapezium method with a step size of $h=0.125$ so as to obtain approximations to $y\left(0.125\right)$ and $y\left(0.25\right)$ .

Answer

Example 4

The current $i$ in a simple circuit involving a resistor of resistance $R$ and an inductance loop of inductance $L$ with applied voltage $E$ satisfies the differential equation

$L\frac{di}{dt}+Ri=E$

Consider the case where $L=1$ , $R=100$ and $E=1000$ . Given that $i\left(0\right)=0$ use a value of $h=0.001$ in implementation of the trapezium method to approximate the current $i$ at times $t=0.001$ and $t=0.002$ .

Solution

The current $i$ satisfies

$\frac{di}{dt}=1000-100i$

and the trapezium approximation to this is

$i_{n+1}-i_n=\frac{h}{2}\left(2000-100i_{n+1}-100i_n\right)$

Rearranging this for $i_{n+1}$ gives

$i_{n+1}=0.904762i_n+0.952381$

It follows that

where these approximations are given to 6 decimal places.

4.1 Accuracy of the trapezium method

Let us now consider an example with a known solution and consider just how accurate the trapezium method is. Suppose that we look at the same test problem we considered when looking at Euler’s method

$\frac{dy}{dt}=y,y\left(0\right)=1.$

We know that the solution to this problem is $y\left(t\right)=e^t$ , and we now compare exact values with the values given by the trapezium method. For the sake of argument, let us consider approximations to $y\left(t\right)$ at $t=1$ . The exact value is $y\left(1\right)=2.718282$ to 6 decimal places. The following table shows results to 6 decimal places obtained on a spreadsheet program for a selection of choices of $h$ .

$\begin{array}{ccc}\hfill h\hfill & \hfill \text{Trapezium approximation}\hfill & \hfill \text{Difference between exact}\hfill \\ \hfill \hfill & \hfill \text{to }y\left(1\right)=2.718282\hfill & \hfill \text{value and trapezium approximation}\hfill \\ ̲& ̲& ̲\\ \hfill 0.2\hfill & \hfill y_5=2.727413\hfill & \hfill 0.009131\hfill \\ \hfill 0.1\hfill & \hfill y_{10}=2.720551\hfill & \hfill 0.002270\hfill \\ \hfill 0.05\hfill & \hfill y_{20}=2.718848\hfill & \hfill 0.000567\hfill \\ \hfill 0.025\hfill & \hfill y_{40}=2.718423\hfill & \hfill 0.000142\hfill \\ \hfill 0.0125\hfill & \hfill y_{80}=2.718317\hfill & \hfill 0.000035\hfill \\ \hfill \hfill \end{array}$

Notice that each time $h$ is reduced by a factor of $\frac{1}{2}$ , the error reduces by a factor of (approximately) $\frac{1}{4}$ . This observation verifies something we will see in Section 32.2, that is that the error in the trapezium approximation is (approximately) proportional to $h^2$ . This sort of behaviour is called second-order .

Exercises

1. Suppose that $y=y\left(t\right)$ is the solution to the initial value problem

   $\frac{dy}{dt}=t+yy\left(0\right)=3$

   Carry out two time steps of Euler’s method with a step size of $h=0.05$ so as to obtain approximations to $y\left(0.05\right)$ and $y\left(0.1\right)$ .

2. Suppose that $y=y\left(t\right)$ is the solution to the initial value problem

   $\frac{dy}{dt}=1∕\left(t^2+1\right)y\left(0\right)=2$

   Carry out two time steps of the trapezium method with a step size of $h=0.1$ so as to obtain approximations to $y\left(0.1\right)$ and $y\left(0.2\right)$ .

3. Suppose that $y=y\left(t\right)$ is the solution to the initial value problem

   $\frac{dy}{dt}=t^2-yy\left(0\right)=1.5$

   Carry out two time steps of the trapezium method with a step size of $h=0.125$ so as to obtain approximations to $y\left(0.125\right)$ and $y\left(0.25\right)$ .

4. The current $i$ in a simple circuit involving a resistor of resistance $R$ , an inductance loop of inductance $L$ with applied voltage $E$ satisfies the differential equation

   $L\frac{di}{dt}+Ri=E$

   Consider the case where $L=1.5$ , $R=120$ and $E=600$ . Given that $i\left(0\right)=0$ use a value of $h=0.0025$ in implementation of the trapezium method to approximate the current $i$ at times $t=0.0025$ and $t=0.005$ .

Answer