Differential equations

NUMERICAL METHODS
derivatives beyond a particular one will vanish and there is no error in taking
the differences to obtain the derivatives.
The following is copied from the tabulation of a second-degree polynomial f(x) at values
of x from 1 to 12 inclusive:
2, 2, ?, 8, 14, 22, 32, 46, ?, 74, 92, 112.
The entries marked ? were illegible and in addition one error was made in transcription.
Complete and correct the table. Would your procedure have worked if the copying error
had been in f(6)?
Write out the entries again in row (a) below, and where possible calculate first differences
in row (b) and second differences in row (c). Denote the jth entry in row (n) by (n)j .
(a) 2
(b)
(c)
2
0
?
?
?
8
?
?
14
6
?
22
8
2
32
10
2
46
14
4
?
?
?
74
?
?
92
18
?
112
20
2
Because the polynomial is second-degree, the second differences (c)j , which are proportional
to d2 f/dx2 , should be constant, and clearly the constant should be 2. That is, (c)6 should
equal 2 and (b)7 should equal 12 (not 14). Since all the (c)j = 2, we can conclude that
(b)2 = 2, (b)3 = 4, (b)8 = 14, and (b)9 = 16. Working these changes back to row (a) shows
that (a)3 = 4, (a)8 = 44 (not 46), and (a)9 = 58.
The entries therefore should read
(a) 2, 2, 4, 8, 14, 22, 32, 44, 58, 74, 92, 112,
where the amended entries are shown in bold type.
It is easily verified that if the error were in f(6) no two computable entries in row (c)
would be equal, and it would not be clear what the correct common entry should be.
Nevertheless, trial and error might arrive at a self-consistent scheme. 27.6 Differential equations
For the remaining sections of this chapter our attention will be on the solution
of differential equations by numerical methods. Some of the general difficulties
of applying numerical methods to differential equations will be all too apparent.
Initially we consider only the simplest kind of equation – one of first order,
typically represented by
dy
= f(x, y),
dx
(27.60)
where y is taken as the dependent variable and x the independent one. If this
equation can be solved analytically then that is the best course to adopt. But
sometimes it is not possible to do so and a numerical approach becomes the
only one available. In fact, most of the examples that we will use can be solved
easily by an explicit integration, but, for the purposes of illustration, this is an
advantage rather than the reverse since useful comparisons can then be made
between the numerically derived solution and the exact one.
1020
27.6 DIFFERENTIAL EQUATIONS
x
0
0.5
1.0
1.5
2.0
2.5
3.0
0.01
(1)
0.605
0.366
0.221
0.134
0.081
0.049
0.1
(1)
0.590
0.349
0.206
0.122
0.072
0.042
0.5
(1)
0.500
0.250
0.125
0.063
0.032
0.016
h
1.0
(1)
0
0
0
0
0
0
y(exact)
1.5
(1)
−0.500
0.250
−0.125
0.063
−0.032
0.016
2
(1)
−1
1
−1
1
−1
1
3
(1)
−2
4
−8
16
−32
64
(1)
0.607
0.368
0.223
0.135
0.082
0.050
Table 27.10 The solution y of differential equation (27.61) using the Euler
forward difference method for various values of h. The exact solution is also
shown.
27.6.1 Difference equations
Consider the differential equation
dy
= −y,
dx
y(0) = 1,
(27.61)
and the possibility of solving it numerically by approximating dy/dx by a finite
difference along the lines indicated in section 27.5. We start with the forward
difference
dy
yi+1 − yi
,
(27.62)
≈
dx xi
h
where we use the notation of section 27.5 but with f replaced by y. In this
particular case, it leads to the recurrence relation
dy
= yi − hyi = (1 − h)yi .
(27.63)
yi+1 = yi + h
dx i
Thus, since y0 = y(0) = 1 is given, y1 = y(0 + h) = y(h) can be calculated, and
so on (this is the Euler method). Table 27.10 shows the values of y(x) obtained
if this is done using various values of h and for selected values of x. The exact
solution, y(x) = exp(−x), is also shown.
It is clear that to maintain anything like a reasonable accuracy only very small
steps, h, can be used. Indeed, if h is taken to be too large, not only is the accuracy
bad but, as can be seen, for h > 1 the calculated solution oscillates (when it
should be monotonic), and for h > 2 it diverges. Equation (27.63) is of the form
yi+1 = λyi , and a necessary condition for non-divergence is |λ| < 1, i.e. 0 < h < 2,
though in no way does this ensure accuracy.
Part of this difficulty arises from the poor approximation (27.62); its righthand side is a closer approximation to dy/dx evaluated at x = xi + h/2 than to
dy/dx at x = xi . This is the result of using a forward difference rather than the
1021
NUMERICAL METHODS
x
y(estim.)
−0.5 (1.648)
0
(1.000)
0.5 0.648
1.0 0.352
1.5 0.296
2.0 0.056
2.5 0.240
3.0 −0.184
y(exact)
—
(1.000)
0.607
0.368
0.223
0.135
0.082
0.050
Table 27.11 The solution of differential equation (27.61) using the Milne
central difference method with h = 0.5 and accurate starting values.
more accurate, but of course still approximate, central difference. A more accurate
method based on central differences (Milne’s method) gives the recurrence relation
dy
(27.64)
yi+1 = yi−1 + 2h
dx i
in general and, in this particular case,
yi+1 = yi−1 − 2hyi .
(27.65)
An additional difficulty now arises, since two initial values of y are needed.
The second must be estimated by other means (e.g. by using a Taylor series,
as discussed later), but for illustration purposes we will take the accurate value,
y(−h) = exp h, as the value of y−1 . If h is taken as, say, 0.5 and (27.65) is applied
repeatedly, then the results shown in table 27.11 are obtained.
Although some improvement in the early values of the calculated y(x) is
noticeable, as compared with the corresponding (h = 0.5) column of table 27.10,
this scheme soon runs into difficulties, as is obvious from the last two rows of the
table.
Some part of this poor performance is not really attributable to the approximations made in estimating dy/dx but to the form of the equation itself and
hence of its solution. Any rounding error occurring in the evaluation effectively
introduces into y some contamination by the solution of
dy
= +y.
dx
This equation has the solution y(x) = exp x and so grows without limit; ultimately
it will dominate the sought-for solution and thus render the calculations totally
inaccurate.
We have only illustrated, rather than analysed, some of the difficulties associated
with simple finite-difference iteration schemes for first-order differential equations,
1022
27.6 DIFFERENTIAL EQUATIONS
but they may be summarised as (i) insufficiently precise approximations to the
derivatives and (ii) inherent instability due to rounding errors.
27.6.2 Taylor series solutions
Since a Taylor series expansion is exact if all its terms are included, and the limits
of convergence are not exceeded, we may seek to use one to evaluate y1 , y2 , etc.
for an equation
dy
= f(x, y),
dx
(27.66)
when the initial value y(x0 ) = y0 is given.
The Taylor series is
y(x + h) = y(x) + hy (x) +
h2 h3
y (x) + y (3) (x) + · · · .
2!
3!
(27.67)
In the present notation, at the point x = xi this is written
yi+1 = yi + hyi(1) +
h2 (2) h3 (3)
y + yi + · · · .
2! i
3!
But, for the required solution y(x), we know that
dy
= f(xi , yi ),
yi(1) ≡
dx xi
(27.68)
(27.69)
and the value of the second derivative at x = xi , y = yi can be obtained from it:
yi(2) =
∂f dy
∂f
∂f
∂f
+
=
+f .
∂x ∂y dx
∂x
∂y
(27.70)
This process can be continued for the third and higher derivatives, all of which
are to be evaluated at (xi , yi ).
Having obtained expressions for the derivatives yi(n) in (27.67), two alternative
ways of proceeding are open to us:
(i) equation (27.68) is used to evaluate yi+1 , the whole process is repeated to
obtain yi+2 , and so on;
(ii) equation (27.68) is applied several times but using a different value of h
each time, and so the corresponding values of y(x + h) are obtained.
It is clear that, on the one hand, approach (i) does not require so many terms of
(27.67) to be kept, but, on the other hand, the yi (n) have to be recalculated at
each step. With approach (ii), fairly accurate results for y may be obtained for
values of x close to the given starting value, but for large values of h a large
number of terms of (27.67) must be kept. As an example of approach (ii) we
solve the following problem.
1023
NUMERICAL METHODS
x
0
0.1
0.2
0.3
0.4
0.5
y(estim.)
1.0000
1.2346
1.5619
2.0331
2.7254
3.7500
y(exact)
1.0000
1.2346
1.5625
2.0408
2.7778
4.0000
Table 27.12 The solution of differential equation (27.71) using a Taylor series.
Find the numerical solution of the equation
dy
y(0) = 1,
(27.71)
= 2y 3/2 ,
dx
for x = 0.1 to 0.5 in steps of 0.1. Compare it with the exact solution obtained analytically.
Since the right-hand side of the equation does not contain x explicitly, (27.70) is greatly
simplified and the calculation becomes a repeated application of
yi(n+1) =
∂y (n) dy
∂y (n)
=f
.
∂y dx
∂y
The necessary derivatives and their values at x = 0, where y = 1, are given below:
y(0) = 1
1
y = 2y 3/2
2
y = (3/2)(2y 1/2 )(2y 3/2 ) = 6y 2
y (3) = (12y)2y 3/2 = 24y 5/2
y (4) = (60y 3/2 )2y 3/2 = 120y 3
y
(5)
2
= (360y )2y
3/2
= 720y
6
24
7/2
120
720
Thus the Taylor expansion of the solution about the origin (in fact a Maclaurin series) is
y(x) = 1 + 2x +
6 2 24 3 120 4 720 5
x + x +
x +
x +··· .
2!
3!
4!
5!
Hence, y(estim.) = 1 + 2x + 3x2 + 4x3 + 5x4 + 6x5 . Values calculated from this are given
in table 27.12. Comparison with the exact values shows that using the first six terms gives
a value that is correct to one part in 100, up to x = 0.3. 27.6.3 Prediction and correction
An improvement in the accuracy obtainable using difference methods is possible
if steps are taken, sometimes retrospectively, to allow for inaccuracies in approximating derivatives by differences. We will describe only the simplest schemes of
this kind and begin with a prediction method, usually called the Adams method.
1024
27.6 DIFFERENTIAL EQUATIONS
The forward difference estimate of yi+1 , namely
dy
yi+1 = yi + h
= yi + hf(xi , yi ),
dx i
(27.72)
would give exact results if y were a linear function of x in the range xi ≤ x ≤ xi +h.
The idea behind the Adams method is to allow some relaxation of this and
suppose that y can be adequately approximated by a parabola over the interval
xi−1 ≤ x ≤ xi+1 . In the same interval, dy/dx can then be approximated by a linear
function:
dy
≈ a + b(x − xi ) for xi − h ≤ x ≤ xi + h.
f(x, y) =
dx
The values of a and b are fixed by the calculated values of f at xi−1 and xi , which
we may denote by fi−1 and fi :
a = fi ,
Thus
yi+1 − yi ≈
xi +h
xi
b=
fi − fi−1
.
h
(fi − fi−1 )
(x − xi ) dx,
fi +
h
which yields
yi+1 = yi + hfi + 12 h(fi − fi−1 ).
(27.73)
The last term of this expression is seen to be a correction to result (27.72). That
it is, in some sense, the second-order correction,
1 2 (2)
2 h yi−1/2 ,
to a first-order formula is apparent.
Such a procedure requires, in addition to a value for y0 , a value for either y1 or
y−1 , so that f1 or f−1 can be used to initiate the iteration. This has to be obtained
by other methods, e.g. a Taylor series expansion.
Improvements to simple difference formulae can also be obtained by using
correction methods. In these, a rough prediction of the value yi+1 is made first,
and then this is used in a better formula, not originally usable since it, in turn,
requires a value of yi+1 for its evaluation. The value of yi+1 is then recalculated,
using this better formula.
Such a scheme based on the forward difference formula might be as follows:
(i) predict yi+1 using yi+1 = yi + hfi ;
(ii) calculate fi+1 using this value;
(iii) recalculate yi+1 using yi+1 = yi + h(fi + fi+1 )/2. Here (fi + fi+1 )/2 has
replaced the fi used in (i), since it better represents the average value of
dy/dx in the interval xi ≤ x ≤ xi + h.
1025
NUMERICAL METHODS
Steps (ii) and (iii) can be iterated to improve further the approximation to the
average value of dy/dx, but this will not compensate for the omission of higherorder derivatives in the forward difference formula.
Many more complex schemes of prediction and correction, in most cases
combining the two in the same process, have been devised, but the reader is
referred to more specialist texts for discussions of them. However, because it
offers some clear advantages, one group of methods will be set out explicitly in
the next subsection. This is the general class of schemes known as Runge–Kutta
methods.
27.6.4 Runge–Kutta methods
The Runge–Kutta method of integrating
dy
= f(x, y)
dx
(27.74)
is a step-by-step process of obtaining an approximation for yi+1 by starting from
the value of yi . Among its advantages are that no functions other than f are used,
no subsidiary differentiation is needed and no additional starting values need be
calculated.
To be set against these advantages is the fact that f is evaluated using somewhat
complicated arguments and that this has to be done several times for each increase
in the value of i. However, once a procedure has been established, for example
on a computer, the method usually gives good results.
The basis of the method is to simulate the (accurate) Taylor series for y(xi + h),
not by calculating all the higher derivatives of y at the point xi but by taking
a particular combination of the values of the first derivative of y evaluated at
a number of carefully chosen points. Equation (27.74) is used to evaluate these
derivatives. The accuracy can be made to be up to whatever power of h is desired,
but, naturally, the greater the accuracy, the more complex the calculation, and,
in any case, rounding errors cannot ultimately be avoided.
The setting up of the calculational scheme may be illustrated by considering
the particular case in which second-order accuracy in h is required. To second
order, the Taylor expansion is
h2 df
,
(27.75)
yi+1 = yi + hfi +
2 dx xi
where
df
dx
=
xi
∂f
∂f
+f
∂x
∂y
≡
xi
∂fi
∂fi
+ fi
,
∂x
∂y
the last step being merely the definition of an abbreviated notation.
1026
27.6 DIFFERENTIAL EQUATIONS
We assume that this can be simulated by a form
yi+1 = yi + α1 hfi + α2 hf(xi + β1 h, yi + β2 hfi ),
(27.76)
which in effect uses a weighted mean of the value of dy/dx at xi and its value at
some point yet to be determined. The object is to choose values of α1 , α2 , β1 and
β2 such that (27.76) coincides with (27.75) up to the coefficient of h2 .
Expanding the function f in the last term of (27.76) in a Taylor series of its
own, we obtain
f(xi + β1 h, yi + β2 hfi ) = f(xi , yi ) + β1 h
∂fi
∂fi
+ β2 hfi
+ O(h2 ).
∂x
∂y
Putting this result into (27.76) and rearranging in powers of h, we obtain
∂fi
∂fi
+ β2 fi
yi+1 = yi + (α1 + α2 )hfi + α2 h2 β1
.
∂x
∂y
(27.77)
Comparing this with (27.75) shows that there is, in fact, some freedom remaining
in the choice of the α’s and β’s. In terms of an arbitrary α1 (= 1),
α2 = 1 − α1 ,
β1 = β2 =
1
.
2(1 − α1 )
One possible choice is α1 = 0.5, giving α2 = 0.5, β1 = β2 = 1. In this case the
procedure (equation (27.76)) can be summarised by
yi+1 = yi + 12 (a1 + a2 ),
(27.78)
where
a1 = hf(xi , yi ),
a2 = hf(xi + h, yi + a1 ).
Similar schemes giving higher-order accuracy in h can be devised. Two such
schemes, given without derivation, are as follows.
(i) To order h3 ,
yi+1 = yi + 16 (b1 + 4b2 + b3 ),
where
b1 = hf(xi , yi ),
b2 = hf(xi + 12 h, yi + 12 b1 ),
b3 = hf(xi + h, yi + 2b2 − b1 ).
1027
(27.79)