- Email: [email protected]
New explicit and implicit “improved euler” methods for the integration of ordinary differential equations
them. Engng, Vol. Compur. PrInted III Great Brltain. All
NEW
12. No. rights
11, pp. reserved
EXPLICIT METHODS ORDINARY
1083-1086,
1988 Copyright
(Received
53.00 + 0.00 0098-1354/88 1988 Pergamon Press plc
AND IMPLICIT “IMPROVED EULER” FOR THE INTEGRATION OF DIFFERENTIAL EQUATIONS T.
0. Department
c
HANNA
of Chemical and Nuclear Engineering, University of California, Santa Barbara, CA 93106, U.S.A.
2 October
1987; final revision
received
7 March
1988;
received for
publication
7 April
1988)
Abatrnct~The new improved Euler methods given here offer several advantages for the solution of ordinary differential equations. The explicit form is of second-order global accuracy but requires only one derivative evaluation per step. It is therefore more efficient than either simple Euler or Runge-Kutta two. This new method is conditionally stable, provides an estimate of local accuracy for error control, is well-suited to global extrapolation and is applicable to real-time problems. Potential important applications also include optimization, continuation problems and parameter estimation. The implicit form is unconditionally stable, offering second-order global accuracy with a stability which is in between backward Euler and the trapezoidal rule. This method should be valuable for stiff problems, and in particular it should serve as an improvement to the well-known Crank-Nicolson method for partial differential equations. INTRODUCTION
The purpose of this paper is to describe a new, very simple and efficient method for the integration of systems of nonstiff ordinary differential equations. The essential feature of this new single-step, explicit method is that it requires only one derivative evaluation per step but it has a second-order global accuracy. For problems where the derivative evaluation is the main computational expense, the new method requires only the cost of a first-order forward Euler calculation, but achieves qualitatively the same accuracy as a Runge-Kutta-two procedure which requires two derivative evaluations per step. Thus this new explicit “improved Euler” method can be viewed as a fundamental improvement of both the Euler and Runge-Kutta-two procedures, since it has the best attributes of each. There is one small price to be paid for this benefit. The stability of the new method will be shown to be somewhat inferior, but in nonstiff problems this is of little consequence. In addition to the explicit form referred to above, an implicit form of the new method will also be given, and applications will be discussed where each form may be useful. THE
BASIC
METHOD
The new method proposed here is developed as follows. We consider integration of a single equation Y’ =f(x, Y); extension to systems of such equations is immediate. The basic idea is to formally integrate the equation and then imagine implementing one step of a successive approximation procedure. This is illustrated in equation (1) where yO(x) is the initial approximation:
By the theory of successive approximation, if we were to evaluate (1) numerically with sufficient accuracy, then y, (x) would generally be a better approximation to the true solution than yO(x), at least for some range of x. Now suppose that we regard yE(x, h), the forward Euler approximation, to be y”(x), the first approximation to the true solution in equation (1). Since the local error of the forward Euler method is O(h*), if the integral in equation (1) is evaluated numerically be a method having local accuracy 0(h3), then the approximation y,(x) will also have local accuracy O(h3). A simple quadrature method which is most useful in this context is the trapezoidal rule. By writing equation (1) between the points x and x + h, using yE(x, h) in place of yO( x), we get:
= Y(O) +
s
‘f[t.
0
x,(t)1 dr.
(1)
(2)
In equations (1) and (2), it is supposed that the entire y,(x, h) trajectory is known wherever needed as an input to the integral. To use the above successive approximation idea in a more advantageous way, it is instead implemented by calculating the simple Euler trajectory as we go along, updating it through the quadrature improvement in equation (2) as given by the trapezoidal rule. Thus, rather than calculating a separate Euler trajectory as an input to equation (2) we implement the following recurrence relations: z(x
+ h) =_Y(x)
calculate Y(X
Initial conditions:
1083
f[x
+ hf[x,
2(X)1,
+ h, Z(X + h)],
(3a) (3b)
+h)=Y(x)+; x {_fk
Y,(X)
x+ h I-It. J-L(I, h)] dt. f V
J’,(-x +h)=y,(x)+
z(x)1 +_m z(0) = y(O) =
+ k 2(x + h)l). Y(0).
(3c)
1084
0. T. HANNA
Equation (3a) is just a simple Euler step, starting from the improved value y(x) and using a slope f[x, r(x)] evaluated using the Euler value z(x). Equation (3~) is just the trapezoidal rule applied to equation (Z), where both slopes are evaluated using the Euler values z(x), z(x + h). The key point here is that a slope at a particular x is evaluated only once, using for Y the Euler value z(x). The slopes are never subsequently corrected after the improved value y(x) has been obtained from equation (3~). By deliberately failing to compute the slopes [using y(xJJ, a very substantial “corrected” eflciency is achieved. For the “cost” of only one derivative evaluation per step, it is easily seen that the new method achieves a global accuracy of O(h2), compared to O(h) for simple Euler. Both the Euler and new improved Euler methods cost one derivative evaluation per step. It should be observed that the calculation sequence in equations (3) would be just a particular Runge-Kutta-two (R-K-2) method, if the derivatives were updated after the correction of y. However, this would then require two derivative evaluations per step. Since the R-K-2 method has a global accuracy of 0(h2), it would be expected that R-K-2 and the new method would have comparable accuracy in practice. This will be illustrated later by typical numerical comparisons. The preceding development illustrates how the new method applies to a single equation. For a system of equations of the type Y; =x(x, Y,, Y,. . Y,) equations (3) become: 2,(x + h) =Y,(x>
+ hflb-,
2, (XL
z&).
. .
,
z,(x)l.
(4a)
calculate J;[x+-h,z,(x+h),
z,(x+h)
+f;[x+h,z,(x+h)
,...,
,...
z,(x+h)],
,z,(x+h)]}.
(4b)
(4~)
Just as for a single equation, in both theory and practice it is easily shown that the new approximations y,(x) have a global error O(h’) for a “cost” of only one derivative evaluation per equation per step. To utilize the advantages of the new method proposed here, it is necessary to employ the following computation sequence. We always assume that the derivatives necessary to do the new Euler step are available from the previous step. Thus we first use these old derivatives to calculate the new Euler values zi(x + h) from equation (4a). Then we evaluate a11 derivatives once during the current step to get /;[x +&2,(x +hI, z2(x + h). . .I. Finally, equation (4~) is used to calculate the improved y, values; this completes one integration step. The basis of the efficiency of the procedure is that each derivative
evaluation (after the first step) is used twice in equations (4); once in (4~) for the previous step and once in (4a) for the current step. It can be shown that both z,(x) and vi(x) are approximations to Y,(X) having global errors of O(h2). The per-step errors of z,(x) and y,(x) are 0(h2) and 0(h3), respectively. It is therefore convenient, if desired, to estimate the local error in terms of the difference zi(x) - Y,(X). Since this difference is O(h*) for small h, one can adjust the step-size on the basis of the formula: desired error ‘J2.h (5) current error > curren“ Here the value [z&r) -y,(x)] is used for “error current”. This procedure has been investigated in the thesis of Hollan (1975), where it was shown to work quite satisfactorily. The procedure has also worked well in several proprietary large-scale simulations. The stability of the new method referred to the prototype equation y’ = -KY can be determined in the usual manner. The resulting second-order difference expression has the characteristic equation of - (Kh/2) = 0. Investigation 1’ - [l - (3/2)Kh]l this equation shows that the real stability limit for the new method is Kh < 1. This result is somewhat inferior to the requirement Kh < 2 for both the simple Euler and R-K-2 methods. However, this difference is not very significant in nonstiff problems since the accuracy in such problems is limited by truncation error, not stability. It is of interest to note that by “averaging” the improved Euler u(x) and z(x) values according to a *v(x) + (1 - a).z(x) we can produce a first-order accurate solution whose stability limit is nearly four times that of the improved Euler procedure. It should also be mentioned that the new improved Euler method can also be implemented for certain so-called “real time” problems, where the algorithm approximating Y(x + h) is not allowed to use any derivative information at (x + h). In such a situation, the new procedure can be successfully operated by using z(x + h) as the desired approximation at time (x + h) and then, after time (x + h), subsequently evaluating y(x + h) to continue the computation sequence toward z(x + 2h), etc. This is useful because, as mentioned earlier, z(x) has a global accuracy of O(h2). h new=
PERFORMANCE
AND
COMPARISON
In this section the new method will be compared to the Euler and R-K-2 procedures for a simple but typical problem. The comparison will be made for equal step-sizes as well as for the case of an equal number of derivative evaluations (equal cost). The system of two first-order equations indicated in Table 1 represents a common chemical flow reactor problem. The results which are shown are entirely typical of similar comparisons for other systems of equa-
New improved
Euler
methods for ordinary differential equations
Table 1. Error comparison of methods for the solution at x = 0.8 of the system: y; = -yy + y2, y,(O) = I; y; = -2y, + y:, y*(O) = 0. Exact value of y,(O.S) (to five figures) is 0.64623; error in ~~(0.8) for indicated method Step
Euler
size
0.4 0.2 0.1 0.05 0.025 0.0125
0.046 0.023 0.011 0.0057 0.0028 0.0014
0.00625
0.0007
New
Y = y(h)
-0.0403 -O.M)SS -0.0011 - 0.0003 -0.00006 - 0.00002 -0.ooOO1
OF
GLOBAL
EXTRAPOLATION
It is natural to inquire whether the intrinsic advantages of the new method might be extended to higher order accuracy through the use of extrapolation to the limit. It is well-known (Henrici, 1962, 1964; Hall and Watt, 1976) that in principle, extrapolation to the limit can be used in conjunction with a number of numerical methods to provide improvement of accuracy and error estimation. Here we explore the use of globat extrapolation, which involves extrapolating the results of two (or more) independent trajectories. “local” extrapolation updates the imBy contrast, proved values as the calculation proceeds (as, for midpoint rule”, Hall and example, in the “modified Watt 1976). The great advantage of global extrapolation is that it can be used to estimate global errors. Only a brief consideration of this matter will be presented here; a detailed development of this idea,
Table
3.
Comparison New
5
of new
method
(h = 0.05, 0.1, 2.0 0.301215 -0.306229 -0.004713 0.042064 -0.004320
0.2)
method
and
New
method
(h = 0.025,
+ C,hZ + C,h’ + O(V).
(6)
This formula applies under broad conditions for any dependent variable at any particular fixed value of x; the coefficients C,, C, depend on x but are independent of h. Y(x) is the true unknown solution and y(h)
tions. By studying the results in Table 1, we can see the following trends. For any given step-size, the accuracy of R-K-2 and the new method are comparable, with R-K-2 generally being a little more accurate. For the same step-size, the accuracy of the Euler method is generally much lower than that of either the R-K-2 or new methods. In terms of accuracy, for the same number of derivative evaluations, we must compare the errors of Euler and the new method for a given step-size vs those of R-K-2 at twice this value. On this fired-cosr basis (in terms of derivative evaluations), it is seen that the accuracy of the new method far exceeds that of Euler and is significantly higher than that of R-K-2. This conclusion is predicted by theory since the Euler method has first-order accuracy while the R-K-2 and new methods have secondorder accuracy. This same qualitative behavior has been observed by the author in virtually 100% of many other similar comparisons.
USE
along with practical applications, is given in Hanna (1985). The basic global extrapolation formula which applies to the new improved Euler method is:
R-K-2
method
-0.0212 -0.0101 -0.0024 -0.00@4 -0.00013 -0.00004 -0.00001
108.5
0.05, 2.0 0.301156 -0.306247 - 0.004586 0.41992 - 0.004339
is the approximate numerical solution corresponding to the fixed step size h. By writing equation (6) for two independent trajectories of step size h and 2h, we can solve for C,h2 to form the once-extrapolated approximation: Y =y(h)
+
Y(h) -Y(2h) 3
+ O(hx)
= y(h, 2h) + O(h’).
(7)
The once-extrapolated improved Euler approximation, y(h, 2h), has global accuracy 0(/z”). By repeating this process using three trajectories, we can estimate the O(h’) term in equation (7) and thereby produce the twice-extrapolated approximation: y(h, 2h, 4h) = y(h, 2h) +
y(h, 2h) - y(2&
4h)
(8)
7
Then: Y(x)=y(h,
2h,4h)-t
O(V).
(9)
The considerable advantage of using the improved Euler method as the basic building block in this process is that it only requires one derivative evaluation per step to achieve second-order accuracy. An example of how improved Euler with extrapolation (IEX) works on a typical system is shown in Table 2. Here we see that the extrapolation process does indeed produce higher accuracies at a lower cost in derivative evaluations. In Table 3 is shown a typical comparison between IEX and the famous classical R-K-4 method. The problem is from Hornbeck (I 975). The comparison shows that the IEX method Table 2. Integration of: y; = -y: +y2, y,(O) = I; y; = -2y2 +yj, y*(O) = 0 by new method with one and two global extrapolations. Table gives y,(x = 0.8) h
y(h)
0.2 0.I 0.05 0.025 0.015 0.00625
R-K-4
0.6563 0.6486 0.6468 I 0.64638 0.64627 0.64624
for y” + 2y’ + 4y R-K-4
0.1)
h =O.l 2.0 0.301137 -0.306259 -O.O@4572 0.041989 - 0.004342
= 0; y(O)
Ah,
2h)
0.6460 0.64621 0.646230 0.6462335 0.64623403
= 2, y’(O)
y(h.
2h.
0.646231 0.646233 0.64623404 0.64623409
= 0
Exact 2.0 0.301149 -0.306245 -0.004579 0.041987 -0.GO4340
4h)
0. T. HANNA
IO86
has an accuracy comparable to that of R-K-4, as theory predicts. What the above comparisons do not show (this matter is discussed in detail in Hanna, 1985), and what is more important than accuracy improvement, is the ability of the global extrapolation procedure to estimate global errors. IMPLICIT
IMPROVED
EULER
As discussed above, the explicit form of the new improved Euler method offers some very significant advantages for the solution of nonstiff differential equations. An implicit form of this new method appears to be advantageous for stiff problems. The implicit algorithm will be written (for a single equation) as: Z(X + h) = y(x) Y(X +
+ hf[x
NOMENCLATURE ET---
+ h, 2(X + h)],
+f[x + h, z(x
z(x)1 + h)]].
(10)
This extends in an obvious manner to a system of equations. Thus the implicit improved Euler is seen to be the implicit backward Euler calculation followed by trapezoidal improvement retaining the backward Euler derivatives. In this case, of course, the new method offers no advantage in terms of derivative evaluations. However, it does offer a very significant and subtle advantage in regard to a good combination of accuracy and stability. A consideration of the definition of the algorithm, equation (lo), shows that the implicit improved Euler is “in between” the classical implicit methods of backward Euler and trapezoidal rule. The backward Euler is accurate only to O(h) but is “super-stable,” while the trapezoidal rule has a better global accuracy [O(h2)], but is only unconditionally stable. As might be “just barely” hoped from its construction, the implicit improved Euler method turns out to combine the principal features of backward Euler and trapezoidal rule. it is clear that the global From its construction, accuracy of implicit improved Euler is O(h’). To we use the prototype equation analyze stability, y’ = - Ky. When this equation is used with equation (lo), a second-order difference equation is produced whose charactenstic equation is a quadratic of the form: A?-(I-B)i+B=O,
Kh 2(1 + Kh)
C, , C, = Coefficients in asymptotic expansion of approximate solution, equation (6)
h
h) = Y(X) + -{f-b, 2.
(11)
where: BZ
in terms of stability, it is much more like backward Euler than trapezoidal rule. The new implicit method therefore appears to be better suited to the solution of stiff differential equations than either backward Euler or trapezoidal rule. An application where this behavior should be of special interest is in the numerical solution of parabolic partial differential equations. From the above discussion, it would be expected that the new implicit procedure should combine the best features of the fully implicit classical and trapezoidal (Crank-Nicolson) methods. This matter will be reported in the future.
Kh 2(1 + Kh)’
The solution of this equation indicates that the magnitude of ;1 is always less than unity for any Kh > 0, so that the method is unconditionally stable. In particular, for Kh -* CT,we have B +f and the real part of L +f. Thus the implicit improved Euler has the second-order accuracy of the trapezoidal rule, but
f(x,
Y) = Derivative of Y, dTz =/(x,
Y)
h = Integration step size
K = Positive constant in prototype y’= -Ky x = Lndcpendent variable
equation,
yO(x) = Initial approximation in successive approximation scheme, equation (1) y,(x) = First iterate in successive approximation scheme, equation (1) y(x) = Final accepted value in new method y(h) = At location X, approximate numerical solution for step size h j(h, 2h) = At location X, once extrapolated approximate solution, y(h) +
y(h) - y(2h) 3
y(h, 2h, 4h) = At location X, twice extrapolated approximate solution, y(b, 2h) - y(2h, y(h, 2h) + ______ 7
4h) -
ye = Y = Z(X) = 1 =
Forward Euler approximation Unknown true solution of Y’ =f(x, Y) Temporary local Euler value in new method Characteristic root in solution of secondorder difference equation, yn = constant-l” CL= Weighting parameter used to increase stability REFERENCES
Hall G. and J. M. Watt (Ed.), Modern Numerical Merho& for Ordinary Differential Equations. Clarendon Press, Oxford (1976). Hanna 0. T., The use of integral global extrapolation in the numerical solution of ordinary differential equations, Paper 52 g, AfChE Annl Mrg, Chicago (1985). Hen&i P., Discrete Variable Melhods in Ordinary D~flerenriaf Equations. Wiley, New York (1962). Henrici P., Elements of Numerical Analysis. Wiley, New York (1964). Hollan M., Experimental analysis of low order RungeKutta numerical integration techniques with and without step-size control. MS. Thesis, University of California Santa Barbara (1975). Hombeck R., Numerical Methods. Quantum Inc., New York (1975).
NEW
12. No. rights
11, pp. reserved
EXPLICIT METHODS ORDINARY
1083-1086,
1988 Copyright
(Received
53.00 + 0.00 0098-1354/88 1988 Pergamon Press plc
AND IMPLICIT “IMPROVED EULER” FOR THE INTEGRATION OF DIFFERENTIAL EQUATIONS T.
0. Department
c
HANNA
of Chemical and Nuclear Engineering, University of California, Santa Barbara, CA 93106, U.S.A.
2 October
1987; final revision
received
7 March
1988;
received for
publication
7 April
1988)
Abatrnct~The new improved Euler methods given here offer several advantages for the solution of ordinary differential equations. The explicit form is of second-order global accuracy but requires only one derivative evaluation per step. It is therefore more efficient than either simple Euler or Runge-Kutta two. This new method is conditionally stable, provides an estimate of local accuracy for error control, is well-suited to global extrapolation and is applicable to real-time problems. Potential important applications also include optimization, continuation problems and parameter estimation. The implicit form is unconditionally stable, offering second-order global accuracy with a stability which is in between backward Euler and the trapezoidal rule. This method should be valuable for stiff problems, and in particular it should serve as an improvement to the well-known Crank-Nicolson method for partial differential equations. INTRODUCTION
The purpose of this paper is to describe a new, very simple and efficient method for the integration of systems of nonstiff ordinary differential equations. The essential feature of this new single-step, explicit method is that it requires only one derivative evaluation per step but it has a second-order global accuracy. For problems where the derivative evaluation is the main computational expense, the new method requires only the cost of a first-order forward Euler calculation, but achieves qualitatively the same accuracy as a Runge-Kutta-two procedure which requires two derivative evaluations per step. Thus this new explicit “improved Euler” method can be viewed as a fundamental improvement of both the Euler and Runge-Kutta-two procedures, since it has the best attributes of each. There is one small price to be paid for this benefit. The stability of the new method will be shown to be somewhat inferior, but in nonstiff problems this is of little consequence. In addition to the explicit form referred to above, an implicit form of the new method will also be given, and applications will be discussed where each form may be useful. THE
BASIC
METHOD
The new method proposed here is developed as follows. We consider integration of a single equation Y’ =f(x, Y); extension to systems of such equations is immediate. The basic idea is to formally integrate the equation and then imagine implementing one step of a successive approximation procedure. This is illustrated in equation (1) where yO(x) is the initial approximation:
By the theory of successive approximation, if we were to evaluate (1) numerically with sufficient accuracy, then y, (x) would generally be a better approximation to the true solution than yO(x), at least for some range of x. Now suppose that we regard yE(x, h), the forward Euler approximation, to be y”(x), the first approximation to the true solution in equation (1). Since the local error of the forward Euler method is O(h*), if the integral in equation (1) is evaluated numerically be a method having local accuracy 0(h3), then the approximation y,(x) will also have local accuracy O(h3). A simple quadrature method which is most useful in this context is the trapezoidal rule. By writing equation (1) between the points x and x + h, using yE(x, h) in place of yO( x), we get:
= Y(O) +
s
‘f[t.
0
x,(t)1 dr.
(1)
(2)
In equations (1) and (2), it is supposed that the entire y,(x, h) trajectory is known wherever needed as an input to the integral. To use the above successive approximation idea in a more advantageous way, it is instead implemented by calculating the simple Euler trajectory as we go along, updating it through the quadrature improvement in equation (2) as given by the trapezoidal rule. Thus, rather than calculating a separate Euler trajectory as an input to equation (2) we implement the following recurrence relations: z(x
+ h) =_Y(x)
calculate Y(X
Initial conditions:
1083
f[x
+ hf[x,
2(X)1,
+ h, Z(X + h)],
(3a) (3b)
+h)=Y(x)+; x {_fk
Y,(X)
x+ h I-It. J-L(I, h)] dt. f V
J’,(-x +h)=y,(x)+
z(x)1 +_m z(0) = y(O) =
+ k 2(x + h)l). Y(0).
(3c)
1084
0. T. HANNA
Equation (3a) is just a simple Euler step, starting from the improved value y(x) and using a slope f[x, r(x)] evaluated using the Euler value z(x). Equation (3~) is just the trapezoidal rule applied to equation (Z), where both slopes are evaluated using the Euler values z(x), z(x + h). The key point here is that a slope at a particular x is evaluated only once, using for Y the Euler value z(x). The slopes are never subsequently corrected after the improved value y(x) has been obtained from equation (3~). By deliberately failing to compute the slopes [using y(xJJ, a very substantial “corrected” eflciency is achieved. For the “cost” of only one derivative evaluation per step, it is easily seen that the new method achieves a global accuracy of O(h2), compared to O(h) for simple Euler. Both the Euler and new improved Euler methods cost one derivative evaluation per step. It should be observed that the calculation sequence in equations (3) would be just a particular Runge-Kutta-two (R-K-2) method, if the derivatives were updated after the correction of y. However, this would then require two derivative evaluations per step. Since the R-K-2 method has a global accuracy of 0(h2), it would be expected that R-K-2 and the new method would have comparable accuracy in practice. This will be illustrated later by typical numerical comparisons. The preceding development illustrates how the new method applies to a single equation. For a system of equations of the type Y; =x(x, Y,, Y,. . Y,) equations (3) become: 2,(x + h) =Y,(x>
+ hflb-,
2, (XL
z&).
. .
,
z,(x)l.
(4a)
calculate J;[x+-h,z,(x+h),
z,(x+h)
+f;[x+h,z,(x+h)
,...,
,...
z,(x+h)],
,z,(x+h)]}.
(4b)
(4~)
Just as for a single equation, in both theory and practice it is easily shown that the new approximations y,(x) have a global error O(h’) for a “cost” of only one derivative evaluation per equation per step. To utilize the advantages of the new method proposed here, it is necessary to employ the following computation sequence. We always assume that the derivatives necessary to do the new Euler step are available from the previous step. Thus we first use these old derivatives to calculate the new Euler values zi(x + h) from equation (4a). Then we evaluate a11 derivatives once during the current step to get /;[x +&2,(x +hI, z2(x + h). . .I. Finally, equation (4~) is used to calculate the improved y, values; this completes one integration step. The basis of the efficiency of the procedure is that each derivative
evaluation (after the first step) is used twice in equations (4); once in (4~) for the previous step and once in (4a) for the current step. It can be shown that both z,(x) and vi(x) are approximations to Y,(X) having global errors of O(h2). The per-step errors of z,(x) and y,(x) are 0(h2) and 0(h3), respectively. It is therefore convenient, if desired, to estimate the local error in terms of the difference zi(x) - Y,(X). Since this difference is O(h*) for small h, one can adjust the step-size on the basis of the formula: desired error ‘J2.h (5) current error > curren“ Here the value [z&r) -y,(x)] is used for “error current”. This procedure has been investigated in the thesis of Hollan (1975), where it was shown to work quite satisfactorily. The procedure has also worked well in several proprietary large-scale simulations. The stability of the new method referred to the prototype equation y’ = -KY can be determined in the usual manner. The resulting second-order difference expression has the characteristic equation of - (Kh/2) = 0. Investigation 1’ - [l - (3/2)Kh]l this equation shows that the real stability limit for the new method is Kh < 1. This result is somewhat inferior to the requirement Kh < 2 for both the simple Euler and R-K-2 methods. However, this difference is not very significant in nonstiff problems since the accuracy in such problems is limited by truncation error, not stability. It is of interest to note that by “averaging” the improved Euler u(x) and z(x) values according to a *v(x) + (1 - a).z(x) we can produce a first-order accurate solution whose stability limit is nearly four times that of the improved Euler procedure. It should also be mentioned that the new improved Euler method can also be implemented for certain so-called “real time” problems, where the algorithm approximating Y(x + h) is not allowed to use any derivative information at (x + h). In such a situation, the new procedure can be successfully operated by using z(x + h) as the desired approximation at time (x + h) and then, after time (x + h), subsequently evaluating y(x + h) to continue the computation sequence toward z(x + 2h), etc. This is useful because, as mentioned earlier, z(x) has a global accuracy of O(h2). h new=
PERFORMANCE
AND
COMPARISON
In this section the new method will be compared to the Euler and R-K-2 procedures for a simple but typical problem. The comparison will be made for equal step-sizes as well as for the case of an equal number of derivative evaluations (equal cost). The system of two first-order equations indicated in Table 1 represents a common chemical flow reactor problem. The results which are shown are entirely typical of similar comparisons for other systems of equa-
New improved
Euler
methods for ordinary differential equations
Table 1. Error comparison of methods for the solution at x = 0.8 of the system: y; = -yy + y2, y,(O) = I; y; = -2y, + y:, y*(O) = 0. Exact value of y,(O.S) (to five figures) is 0.64623; error in ~~(0.8) for indicated method Step
Euler
size
0.4 0.2 0.1 0.05 0.025 0.0125
0.046 0.023 0.011 0.0057 0.0028 0.0014
0.00625
0.0007
New
Y = y(h)
-0.0403 -O.M)SS -0.0011 - 0.0003 -0.00006 - 0.00002 -0.ooOO1
OF
GLOBAL
EXTRAPOLATION
It is natural to inquire whether the intrinsic advantages of the new method might be extended to higher order accuracy through the use of extrapolation to the limit. It is well-known (Henrici, 1962, 1964; Hall and Watt, 1976) that in principle, extrapolation to the limit can be used in conjunction with a number of numerical methods to provide improvement of accuracy and error estimation. Here we explore the use of globat extrapolation, which involves extrapolating the results of two (or more) independent trajectories. “local” extrapolation updates the imBy contrast, proved values as the calculation proceeds (as, for midpoint rule”, Hall and example, in the “modified Watt 1976). The great advantage of global extrapolation is that it can be used to estimate global errors. Only a brief consideration of this matter will be presented here; a detailed development of this idea,
Table
3.
Comparison New
5
of new
method
(h = 0.05, 0.1, 2.0 0.301215 -0.306229 -0.004713 0.042064 -0.004320
0.2)
method
and
New
method
(h = 0.025,
+ C,hZ + C,h’ + O(V).
(6)
This formula applies under broad conditions for any dependent variable at any particular fixed value of x; the coefficients C,, C, depend on x but are independent of h. Y(x) is the true unknown solution and y(h)
tions. By studying the results in Table 1, we can see the following trends. For any given step-size, the accuracy of R-K-2 and the new method are comparable, with R-K-2 generally being a little more accurate. For the same step-size, the accuracy of the Euler method is generally much lower than that of either the R-K-2 or new methods. In terms of accuracy, for the same number of derivative evaluations, we must compare the errors of Euler and the new method for a given step-size vs those of R-K-2 at twice this value. On this fired-cosr basis (in terms of derivative evaluations), it is seen that the accuracy of the new method far exceeds that of Euler and is significantly higher than that of R-K-2. This conclusion is predicted by theory since the Euler method has first-order accuracy while the R-K-2 and new methods have secondorder accuracy. This same qualitative behavior has been observed by the author in virtually 100% of many other similar comparisons.
USE
along with practical applications, is given in Hanna (1985). The basic global extrapolation formula which applies to the new improved Euler method is:
R-K-2
method
-0.0212 -0.0101 -0.0024 -0.00@4 -0.00013 -0.00004 -0.00001
108.5
0.05, 2.0 0.301156 -0.306247 - 0.004586 0.41992 - 0.004339
is the approximate numerical solution corresponding to the fixed step size h. By writing equation (6) for two independent trajectories of step size h and 2h, we can solve for C,h2 to form the once-extrapolated approximation: Y =y(h)
+
Y(h) -Y(2h) 3
+ O(hx)
= y(h, 2h) + O(h’).
(7)
The once-extrapolated improved Euler approximation, y(h, 2h), has global accuracy 0(/z”). By repeating this process using three trajectories, we can estimate the O(h’) term in equation (7) and thereby produce the twice-extrapolated approximation: y(h, 2h, 4h) = y(h, 2h) +
y(h, 2h) - y(2&
4h)
(8)
7
Then: Y(x)=y(h,
2h,4h)-t
O(V).
(9)
The considerable advantage of using the improved Euler method as the basic building block in this process is that it only requires one derivative evaluation per step to achieve second-order accuracy. An example of how improved Euler with extrapolation (IEX) works on a typical system is shown in Table 2. Here we see that the extrapolation process does indeed produce higher accuracies at a lower cost in derivative evaluations. In Table 3 is shown a typical comparison between IEX and the famous classical R-K-4 method. The problem is from Hornbeck (I 975). The comparison shows that the IEX method Table 2. Integration of: y; = -y: +y2, y,(O) = I; y; = -2y2 +yj, y*(O) = 0 by new method with one and two global extrapolations. Table gives y,(x = 0.8) h
y(h)
0.2 0.I 0.05 0.025 0.015 0.00625
R-K-4
0.6563 0.6486 0.6468 I 0.64638 0.64627 0.64624
for y” + 2y’ + 4y R-K-4
0.1)
h =O.l 2.0 0.301137 -0.306259 -O.O@4572 0.041989 - 0.004342
= 0; y(O)
Ah,
2h)
0.6460 0.64621 0.646230 0.6462335 0.64623403
= 2, y’(O)
y(h.
2h.
0.646231 0.646233 0.64623404 0.64623409
= 0
Exact 2.0 0.301149 -0.306245 -0.004579 0.041987 -0.GO4340
4h)
0. T. HANNA
IO86
has an accuracy comparable to that of R-K-4, as theory predicts. What the above comparisons do not show (this matter is discussed in detail in Hanna, 1985), and what is more important than accuracy improvement, is the ability of the global extrapolation procedure to estimate global errors. IMPLICIT
IMPROVED
EULER
As discussed above, the explicit form of the new improved Euler method offers some very significant advantages for the solution of nonstiff differential equations. An implicit form of this new method appears to be advantageous for stiff problems. The implicit algorithm will be written (for a single equation) as: Z(X + h) = y(x) Y(X +
+ hf[x
NOMENCLATURE ET---
+ h, 2(X + h)],
+f[x + h, z(x
z(x)1 + h)]].
(10)
This extends in an obvious manner to a system of equations. Thus the implicit improved Euler is seen to be the implicit backward Euler calculation followed by trapezoidal improvement retaining the backward Euler derivatives. In this case, of course, the new method offers no advantage in terms of derivative evaluations. However, it does offer a very significant and subtle advantage in regard to a good combination of accuracy and stability. A consideration of the definition of the algorithm, equation (lo), shows that the implicit improved Euler is “in between” the classical implicit methods of backward Euler and trapezoidal rule. The backward Euler is accurate only to O(h) but is “super-stable,” while the trapezoidal rule has a better global accuracy [O(h2)], but is only unconditionally stable. As might be “just barely” hoped from its construction, the implicit improved Euler method turns out to combine the principal features of backward Euler and trapezoidal rule. it is clear that the global From its construction, accuracy of implicit improved Euler is O(h’). To we use the prototype equation analyze stability, y’ = - Ky. When this equation is used with equation (lo), a second-order difference equation is produced whose charactenstic equation is a quadratic of the form: A?-(I-B)i+B=O,
Kh 2(1 + Kh)
C, , C, = Coefficients in asymptotic expansion of approximate solution, equation (6)
h
h) = Y(X) + -{f-b, 2.
(11)
where: BZ
in terms of stability, it is much more like backward Euler than trapezoidal rule. The new implicit method therefore appears to be better suited to the solution of stiff differential equations than either backward Euler or trapezoidal rule. An application where this behavior should be of special interest is in the numerical solution of parabolic partial differential equations. From the above discussion, it would be expected that the new implicit procedure should combine the best features of the fully implicit classical and trapezoidal (Crank-Nicolson) methods. This matter will be reported in the future.
Kh 2(1 + Kh)’
The solution of this equation indicates that the magnitude of ;1 is always less than unity for any Kh > 0, so that the method is unconditionally stable. In particular, for Kh -* CT,we have B +f and the real part of L +f. Thus the implicit improved Euler has the second-order accuracy of the trapezoidal rule, but
f(x,
Y) = Derivative of Y, dTz =/(x,
Y)
h = Integration step size
K = Positive constant in prototype y’= -Ky x = Lndcpendent variable
equation,
yO(x) = Initial approximation in successive approximation scheme, equation (1) y,(x) = First iterate in successive approximation scheme, equation (1) y(x) = Final accepted value in new method y(h) = At location X, approximate numerical solution for step size h j(h, 2h) = At location X, once extrapolated approximate solution, y(h) +
y(h) - y(2h) 3
y(h, 2h, 4h) = At location X, twice extrapolated approximate solution, y(b, 2h) - y(2h, y(h, 2h) + ______ 7
4h) -
ye = Y = Z(X) = 1 =
Forward Euler approximation Unknown true solution of Y’ =f(x, Y) Temporary local Euler value in new method Characteristic root in solution of secondorder difference equation, yn = constant-l” CL= Weighting parameter used to increase stability REFERENCES
Hall G. and J. M. Watt (Ed.), Modern Numerical Merho& for Ordinary Differential Equations. Clarendon Press, Oxford (1976). Hanna 0. T., The use of integral global extrapolation in the numerical solution of ordinary differential equations, Paper 52 g, AfChE Annl Mrg, Chicago (1985). Hen&i P., Discrete Variable Melhods in Ordinary D~flerenriaf Equations. Wiley, New York (1962). Henrici P., Elements of Numerical Analysis. Wiley, New York (1964). Hollan M., Experimental analysis of low order RungeKutta numerical integration techniques with and without step-size control. MS. Thesis, University of California Santa Barbara (1975). Hombeck R., Numerical Methods. Quantum Inc., New York (1975).