Applied Numerical Mathematics 31 (1999) 227238

RungeKutta methods for linear ordinary differential equations

Abstract Three new RungeKutta methods are presented for numerical integration of systems of linear inhomogeneous ordinary differential equations (ODEs) with constant coefcients. Such ODEs arise in the numerical solution of partial differential equations governing linear wave phenomena. The restriction to linear ODEs with constant coefcients reduces the number of conditions which the coefcients of the RungeKutta method must satisfy. This freedom is used to develop methods which are more efcient than conventional RungeKutta methods. A fourthorder method is presented which uses only two memory locations per dependent variable, while the classical fourth-order RungeKutta method uses three. This method is an excellent choice for simulations of linear wave phenomena if memory is a primary concern. In addition, fth- and sixth-order methods are presented which require ve and six stages, respectively, one fewer than their conventional counterparts, and are therefore more efcient. These methods are an excellent option for use with high-order spatial discretizations. 1999 Elsevier Science B.V. and IMACS. All rights reserved. Keywords: RungeKutta methods; Linear differential equations; Ordinary differential equations

1. Introduction We consider the numerical integration of large linear inhomogenous systems of ordinary differential equations in the form du = Au g(t), (1) dt where A is an M by M matrix whose elements depend on neither u nor t, and u and g(t) are vectors of length M. Such essentially autonomous systems arise in the numerical solution of partial differential equations (PDEs) governing linear wave phenomena after application of a spatial discretization such as a nite-difference, nite-volume, or nite-element method. Examples of such PDEs are the linearized Euler equations governing acoustic waves and the Maxwell equations governing electromagnetic waves. The elements of A depend on the PDE and the spatial discretization. The inhomogeneous term g(t) is
Corresponding author. E-mail: dwz@oddjob.utias.utoronto.ca.

0168-9274/99/$20.00 1999 Elsevier Science B.V. and IMACS. All rights reserved. PII: S 0 1 6 8 - 9 2 7 4 ( 9 8 ) 0 0 1 2 9 - 9

228

D.W. Zingg, T.T. Chisholm / Applied Numerical Mathematics 31 (1999) 227238

associated with either a source term or the boundary conditions. In the context of wave propagation, the system of ODEs is often mildly stiff with the eigenvalues of A typically lying near the imaginary axis. The system of ODEs arising from the application of a spatial discretization to a system of PDEs can be very large, especially in three-dimensional simulations. Consequently, the constraints on the methods used for integrating these systems are somewhat different from those which have driven much of the development of numerical methods for initial value problems. Due to their high accuracy and modest memory requirements, explicit RungeKutta methods have become popular for simulations of wave phenomena [79,18,20]. The classical fourth-order RungeKutta method requires three memory locations per dependent variable [1,6], but low-storage methods requiring only two memory locations per dependent variable can be derived [4,16,17]. This property is easily achieved by a third-order RungeKutta method [17], but an additional stage is required for a fourth-order method [4]. Since the primary cost of the integration is in the evaluation of the derivative function, and each stage requires a function evaluation, the additional stage represents a signicant increase in expense. For the same reason, error checking is generally not performed when solving very large systems of ODEs arising from the discretization of PDEs. There have been several attempts to develop efcient methods for integrating linear systems of ODEs [5,1214]. The basic premise of these methods is that the major cost in evaluating the derivative function is in forming the matrix A and the vector g(t). In the application considered here, the simulation of linear wave phenomena, the matrix A is never explicitly formed or stored. Hence the methods previously proposed for linear systems are not appropriate for this application. It is well known that a RungeKutta method with p stages has an order of accuracy not exceeding p [2, 3]. For p 4, methods of order p can be derived with p stages. However, fth-and sixth-order methods require at least six and seven stages, respectively. Nine stages are required for seventh-order accuracy and eleven for eighth-order accuracy [2]. Since the cost for our application is roughly proportional to the number of stages, this represents a signicant limitation of higher-order RungeKutta methods. Several authors have considered various approximations to reduce the number of stages and the storage requirements of high-order RungeKutta methods. Shanks [15] was able to develop schemes with a reduced number of stages by requiring only that the accuracy conditions be approximately satised. Zingg et al. [19,20] propose methods with low storage requirements which are of high order for linear homogeneous ODEs but second-order otherwise. A similar idea was proposed previously by Lorenz [11]. In this paper, we develop RungeKutta methods specically for linear ODEs with constant coefcients. By removing the constraints imposed by nonlinearity in the derivative function, high-order RungeKutta methods can be derived which are more efcient in some respect than the classical methods. In the next section, we present a fourth-order method which requires less memory than the classical fourth-order RungeKutta method. We then present fth- and sixth-order methods requiring fewer derivative function evaluations per time step than fth- and sixth-order RungeKutta methods applicable to nonlinear problems.

2. General form of an explicit RungeKutta method Without loss of generality, we consider the following scalar ODE: du = f (t, u). dt (2)

D.W. Zingg, T.T. Chisholm / Applied Numerical Mathematics 31 (1999) 227238

229

A general p-stage explicit RungeKutta method can be written as k1 = f (tn , un ),

ki = f tn + ci h, un + h
j =1 p

aij kj ,

i = 2, . . . , p,

(3)

un+1 = un + h
i=1

bi ki ,

where h =

t is the time step, tn = nh, and un is an approximation to u(tn ).

3. Low-storage fourth-order method We consider rst the case p = 4. With the constraints c2 = a21 , c3 = a31 + a32 , c4 = a41 + a42 + a43

(4)

there remain ten parameters. For fourth-order accuracy, there are eight conditions which must be satised. Four of these arise even for linear homogeneous constant-coefcient ODEs. A further three conditions must be met if the ODEs are inhomogeneous. The nal condition is associated with non-constant coefcients or nonlinearity. Therefore, fourth-order RungeKutta methods are a two-parameter family of which the classical method is a particular choice. If we restrict our attention to linear constant-coefcient ODEs, the number of conditions is reduced to seven. They are
4

bi = 1,
i=1 4

ci bi = 1 , 2
i=2

c2 a32 b3 + b4 (c2 a42 + c3 a43 ) = 1 , 6 c2 a32 a43 b4 =

(5)

bi ci2 = 1 , 3
i=2 4

bi ci3 = 1 , 4
i=2 2 2 2 b3 c2 a32 + b4 c2 a42 + c3 a43 = 1 . 12

230

D.W. Zingg, T.T. Chisholm / Applied Numerical Mathematics 31 (1999) 227238

These conditions can be derived in the same manner as the classical conditions. The eighth condition is eliminated by exploiting the fact that, for linear ODEs, 2f 2f = = 0. u2 ut The reduction in the number of conditions to be satised does not permit us to reduce the number of stages. However, we can obtain reduced storage requirements. Following the approach of Wray [17], the requirement that only two memory locations be used imposes the following three constraints: b1 = a41 = a31 , b2 = a42 . (6)

With these constraints, only two memory locations are required for both the dependent variable and the value of the time derivative. Hence the method requires minimal storage even when compact or spectral methods are used for the spatial discretization. With the memory locations denoted A and B, the method proceeds as follows: 1. Initially, un is stored in A, and B is empty. 2. The term k1 = f (tn , un ) is evaluated and stored in B. 3. The quantity un + ha 31 k1 , is calculated and stored in A. 4. The quantity un + ha 21 k1 is calculated and stored in B. 5. The term k2 = f (tn + c2 h, un + ha 21 k1 ) is evaluated and stored in B. 6. The contents of the two memory locations are linearly combined to form un + h(a31 k1 + a32 k2 ), which is stored in B. 7. With a41 = a31 , another linear combination gives un + h(a41 k1 + a42 k2 ), which is stored in A. 8. The term k3 = f [tn + c3 h, un + h(a31 k1 + a32 k2 )] is evaluated and stored in B. 9. The contents of the two memory locations are linearly combined to form un + h(a41 k1 + a42 k2 + a43 k3 ), which is stored in B. 10. With b1 = a41 and b2 = a42 , another linear combination gives un + h(b1 k1 + b2 k2 + b3 k3 ), which is stored in A. 11. The term k4 = f [tn + c4 h, un + h(a41 k1 + a42 k2 + a43 k3 )] is evaluated and stored in B. 12. The contents of the two memory locations are linearly combined to form un+1 . With the additional constraints imposed by the low-storage requirement, we are left with seven parameters to satisfy the seven conditions given in Eq. (5). Although this system may possess more than one solution, the only solution we have found is a21 = c2 = 0.69631521002413, c4 = 0.82502163765, a32 = 0.21640084013679, a43 = 0.69991725920066, b4 = 0.39507289160708. c3 = 0.29441651741, b1 = a41 = a31 = 0.07801567728325, b2 = a42 = 0.04708870117112, b3 = 0.47982272993855,

D.W. Zingg, T.T. Chisholm / Applied Numerical Mathematics 31 (1999) 227238

231

4. Five-stage fth-order method For the case p = 5, we have, in addition to the constraints given in Eq. (4), the following condition: c5 = a51 + a52 + a53 + a54 . (7) Consequently, adding the fth stage has produced ve additional parameters for a total of fteen. The coefcients must satisfy the following eleven conditions in order to produce fth-order accuracy for linear constant-coefcient ODEs:
5

bi = 1,
i=1 5

ci bi = 1 , 2
i=2

c2 a32 b3 + b4 (c2 a42 + c3 a43 ) + b5 (c2 a52 + c3 a53 + c4 a54 ) = 1 , 6 c2 a32 a43 b4 + b5 c2 (a42 a54 + a32 a53 ) + c3 a43 a54 = c2 a32 a43 a54 b5 =
5 1 , 120 1 , 24

bi ci2 = 1 , 3
i=2 5

(8)

bi ci3 = 1 , 4
i=2 2 2 2 2 2 2 b3 c2 a32 + b4 c2 a42 + c3 a43 + b5 c2 a52 + c3 a53 + c4 a54 = 5 1 , 12

bi ci4 = 1 , 5
i=2 3 3 3 (b4 a42 + b3 a32 + b5 a52 )c2 + (b5 a53 + b4 a43 )c3 + b5 a54 c4 = 2 2 2 2 b5 a54 a42 c2 + a43 c3 + a53 a32 c2 + b4 a43 a32 c2 = 1 . 60 1 , 20

Thus a four-parameter family of solutions is obtained. Several different criteria can be applied in order to choose a method from this family. The following values have been found by minimizing the L2 norm of a vector containing the coefcients of the method: a21 = c2 = 0.21, c3 = 0.43, c4 = 0.68, c5 = 0.85, a32 = 0.47418546365915, a52 = 0.26302355344001, a43 = 0.57068167533284, b3 = 0.41041645692809, b4 = 0.04092124960122, b1 = 0.09235969809721 a42 = 0.13437223603429, b2 = 0.16574368303091, a53 = 0.10434139625551, a54 = 0.39377303853165, b5 = 0.37240141154501,

232

D.W. Zingg, T.T. Chisholm / Applied Numerical Mathematics 31 (1999) 227238

with a31 = c3 a32 , a41 = c4 a42 a43 , a51 = c5 a52 a53 a54 .

5. Six-stage sixth-order method With p = 6, the following condition must be satised in addition to the constraints given in Eqs. (4) and (7): c6 = a61 + a62 + a63 + a64 + a65 . (9)

Therefore, there remain twenty-one free coefcients. The requirement of sixth-order accuracy for linear constant-coefcient ODEs produces the following sixteen conditions:
6 6

bi = 1,
i=1 i=2

ci bi = 1 , 2

c2 a32 b3 + b4 (c2 a42 + c3 a43 ) + b5 (c2 a52 + c3 a53 + c4 a54 ) + b6 (c2 a62 + c3 a63 + c4 a64 + c5 a65 ) = 1 , 6 c2 a32 a43 b4 + b5 c2 (a42 a54 + a32 a53 ) + c3 a43 a54
1 , 24 1 b6 a65 a54 (a43 c3 + a42 c2 ) + a53 a32 c2 + a64 a43 a32 c2 + b5 a54 a43 a32 c2 = 120 , 1 c2 a32 a43 a54 a65 b6 = 720 , 6 6 bi ci2 = 1 , bi ci3 = 1 , 3 4 i=2 i=2 2 2 2 2 2 2 b3 c2 a32 + b4 c2 a42 + c3 a43 + b5 c2 a52 + c3 a53 + c4 a54 2 2 2 2 1 + b6 a65 c5 + a64 c4 + a63 c3 + a62 c2 = 12 , 6

+ b6 a65 (a54 c4 + a53 c3 + a52 c2 ) + a64 (a43 c3 + a42 c2 ) + a63 a32 c2 =

(10)

bi ci4 = 1 , 5
i=2 3 3 (b4 a42 + b3 a32 + b5 a52 + b6 a62 )c2 + (b5 a53 + b4 a43 + b6 a63 )c3 3 3 + (b5 a54 + b6 a64 )c4 + b6 a65 c5 = 1 , 20

2 2 2 2 2 2 b6 a65 a52 c2 + a53 c3 + a54 c4 + a64 a42 c2 + a43 c3 + a63 a32 c2 2 2 2 2 + b5 a54 a42 c2 + a43 c3 + a53 a32 c2 + b4 a43 a32 c2 = 6 1 , 60

bi ci5 = 1 , 6
i=2 4 4 4 4 4 4 4 4 4 4 b6 (a62 c2 + a63 c3 + a64 c4 + a65 c5 ) + b5 (a52 c2 + a53 c3 + a54 c4 ) + b4 (a42 c2 + a43 c3 ) + b3 a32 c2 = 3 c2 b6 (a65 a52 + a64 a42 + a63 a32 ) + b5 (a54 a42 + a53 a32 ) + b4 a43 a32 1 , 120 2 2 2 2 1 b6 a65 a54 a43 c3 + a42 c2 + a64 a43 a32 c2 + b5 a54 a43 a32 c2 = 360 . 3 3 + c3 b6 (a65 a53 + a64 a43 ) + b5 a54 a43 + c4 b6 a65 a54 = 1 , 30

D.W. Zingg, T.T. Chisholm / Applied Numerical Mathematics 31 (1999) 227238

233

Using the same criterion as for the fth-order method, the following coefcients have been chosen from the ve-parameter family of solutions to the above conditions: a21 = c2 = 0.15, c3 = 0.36, c4 = 0.57, c5 = 0.75, a32 = 0.45818181818182, a42 = 0.09769454545455, a52 = 0.10861879806510, a62 = 0.20874226393025, b2 = 0.24971305394585, a43 = 0.48766666666667, a53 = 0.04655817933320, a63 = 0.12686271445897, b3 = 0.11278150363005, a54 = 0.44703799502007, a64 = 0.02734417934727, b4 = 0.35718962665957, a65 = 0.37591957583530, b5 = 0.00478351095633, b6 = 0.24659027402511, b1 = 0.03850905269576 with a31 = c3 a32 , a51 = c5 a52 a53 a54 , a41 = c4 a42 a43 , a61 = c6 a62 a63 a64 a65 . c6 = 0.90,

6. Stability contours The stability contours of the three new methods are shown in Fig. 1. Satisfaction of the rst four conditions in Eq. (5) ensures that the new fourth-order method has the same stability contour as the classical fourth-order RungeKutta method. Similarly, the stability contours of the ve-stage fth-order

Fig. 1. Stability contours for the fourth-order (), fth-order (- - -), and sixth-order ( ) methods.

234

D.W. Zingg, T.T. Chisholm / Applied Numerical Mathematics 31 (1999) 227238

method and the six-stage sixth-order method are uniquely dened and do not depend on which members of the respective families are selected. Although the stable regions of the fth- and sixth-order methods are somewhat larger than that of the fourth-order method, the increase is not sufcient to compensate for the cost of the additional stages. Therefore, the fourth-order method is a better choice if the time step is limited by stability considerations. It should be emphasized that the present fourth-order method has exactly the same properties with respect to stiff ODEs as the classical fourth-order RungeKutta method. The stable regions of the fth- and sixth-order methods do not include the imaginary axis. Systems with pure imaginary eigenvalues are obtained when central differencing is applied to the spatial derivatives in partial differential equations governing wave propagation phenomena with no physical dissipation, in the absence of boundary conditions. However, Zingg et al. [20] have demonstrated that by adding a small amount of numerical dissipation to the spatial discretization, stable schemes can be obtained using such methods. The amount of dissipation required is sufciently low that the overall accuracy of the scheme is not compromised. The stability contour of the method successfully used in [9] for simulations of the propagation and scattering of electromagnetic waves is identical to that of the present sixth-order method. 7. Fourier error analysis Using Fourier analysis we can determine the errors produced by an integration method when applied to a linear homogeneous ODE. Since our interest is in wave propagation, we consider a scalar ODE of the form du = iu, (11) dt where is a real constant. The RungeKutta methods developed here produce a solution in the form un = n u0 , where = 1 (i)k k! k=0
p

(12)

(13)

and p is the number of stages. The local amplitude and phase errors are determined from as follows: era = | | 1, erp =
1

(14)

tan (i /r ) + 1, (15) h where r and i denote the real and imaginary parts of . Figs. 2 and 3 show the local amplitude and phase errors produced by the three new methods. In order to account for the number of stages, the errors are plotted versus h/p. Hence the errors shown are for approximately equal computational effort. Since the time step is thus proportional to p, the amplitude error shown is | |1/p 1. The gures show that each increase in the order of the method produces an increase in accuracy even though the extra work has been accounted for. Hence the fth- and sixth-order methods can be more efcient than the fourth-order method if a sufciently accurate spatial discretization is used.

D.W. Zingg, T.T. Chisholm / Applied Numerical Mathematics 31 (1999) 227238

235

Fig. 2. Amplitude error produced by the fourth-order (), fth-order (- - -), and sixth-order ( ) methods.

Fig. 3. Phase error produced by the fourth-order (), fth-order (- - -), and sixth-order ( ) methods.

8. Application to the numerical solution of a PDE We now apply the new methods to the numerical solution of the linear convection equation given by u u + = 0, t x (16)

236

D.W. Zingg, T.T. Chisholm / Applied Numerical Mathematics 31 (1999) 227238

Fig. 4. Root-mean-square error as a function of the number of grid nodes for the numerical solution of the linear convection equation.

where u = u(x, t) and 0 x 1. This equation is commonly used as a model equation in the study of numerical algorithms for the solution of hyperbolic systems of PDEs. The boundary condition is u(0, t) = sin t with = 16 . We run until a periodic steady state is obtained. The spatial operator is seventh-order, with a nine-point stencil and suitable boundary operators. The error is displayed as a function of the number of grid nodes in Fig. 4. The ratio t/ x is unity in all cases. The three methods presented here are designated linear. In addition, the classical fourth-order RungeKutta method is shown, along with a six-stage fth-order method [3, p. 202] and an eight-stage sixth-order method [10, p. 122]. The error produced by the new low-storage fourth-order method is indistinguishable from that produced by the classical fourth-order method. The linear fth-order method produces somewhat more error than the six-stage fth-order method, but is much more accurate than the fourth-order methods. The two sixth-order methods produce comparable errors, with the linear method requiring two fewer stages. Overall, these results conrm the validity and usefulness of the new methods for use in the numerical solution of linear systems of hyperbolic PDEs, such as the time-domain Maxwell equations.

9. Conclusions Three new RungeKutta methods have been presented for the integration of linear systems of ODEs with constant coefcients. If the time step size is limited by stability, then the new fourth-order method

D.W. Zingg, T.T. Chisholm / Applied Numerical Mathematics 31 (1999) 227238

237

is the most suitable of the new methods. This method requires less memory than the classical fourthorder RungeKutta method and less computational effort than the low-storage method proposed in [4]. In any application to linear ODEs where the classical fourth-order RungeKutta method is currently the method of choice, the new fourth-order method provides the same stability characteristics and equivalent accuracy with reduced memory requirements. If the time step is limited by accuracy, and memory is a secondary concern, then the new fth- and sixth-order methods present an efcient new alternative. Since the expense of the methods is roughly proportional to the number of stages for the problems of interest here, the new fth- and sixth-order methods are signicantly more efcient than their counterparts for nonlinear ODEs. The sixth-order method is a particularly good choice for use with high-order spatial discretizations in the numerical solution of partial differential equations governing linear wave phenomena.

References
[1] E.K. Blum, A modication of the RungeKutta fourth-order method, Math. Comp. 16 (1962) 176187. [2] J.C. Butcher, The non-existence of ten stage eighth order explicit RungeKutta methods, BIT 25 (1985) 521 540. [3] J.C. Butcher, The Numerical Analysis of Ordinary Differential Equations (Wiley, New York, 1987). [4] M.H. Carpenter and C.A. Kennedy, Fourth-order 2N -storage RungeKutta schemes, NASA TM 109112 (June 1994). [5] W.H. Enright, The efcient solution of linear constant-coefcient systems of differential equations, Simulation 30 (1978) 129133. [6] S. Gill, A process for the step-by-step integration of differential equations in an automatic digital computing machine, Proc. Cambridge Philos. Soc. 47 (1951) 96108. [7] Z. Haras and S. TaAsan, Finite-difference schemes for long-time integration, J. Comput. Phys. 114 (1994) 265279. [8] F.Q. Hu, M.Y. Hussaini and J.L. Manthey, Low-dissipation and low-dispersion RungeKutta schemes for computational acoustics, J. Comput. Phys. 124 (1996) 177191. [9] H. Jurgens and D.W. Zingg, Implementation of a high-accuracy nite-difference scheme for linear wave phenomena, Houston J. Math. (July 1996). [10] J.D. Lambert, Computational Methods in Ordinary Differential Equations (Wiley, New York, 1973). [11] E.N. Lorenz, An N -cycle time-differencing scheme for stepwise numerical integration, Monthly Weather Rev. 99 (1971) 644648. [12] L. Randez, Optimizing the numerical integration of initial value problems in shooting methods for linear boundary value problems, SIAM J. Sci. Comput. 14 (1993) 860871. [13] L.F. Shampine, Cheaper integration of linear systems, Simulation 20 (1973) 17. [14] L.F. Shampine, Solution of structured non-stiff ODEs, J. Comput. Appl. Math. 15 (1986) 293300. [15] E.B. Shanks, Higher-order approximations of RungeKutta type, NASA TN D-2920 (September 1965). [16] J.H. Williamson, Low-storage RungeKutta schemes, J. Comput. Phys. 35 (1980) 4856. [17] A.A. Wray, Minimal storage time advancement schemes for spectral methods, Manuscript, NASA Ames Research Center (unpublished). [18] D.W. Zingg, A review of high-order and optimized nite-difference methods for simulating linear wave phenomena, AIAA Paper 97-2088, in: Proceedings of the 13th AIAA Computational Fluid Dynamics Conference, American Institute for Aeronautics and Astronautics (June 1997).

238

D.W. Zingg, T.T. Chisholm / Applied Numerical Mathematics 31 (1999) 227238

[19] D.W. Zingg and H. Lomax, Some aspects of high-order numerical solutions of the linear convection equation with forced boundary conditions, AIAA Paper 93-3381, in: Proceedings of the 11th AIAA Computational Fluid Dynamics Conference, American Institute for Aeronautics and Astronautics (June 1993). [20] D.W. Zingg, H. Lomax and H.M. Jurgens, High-accuracy nite-difference schemes for linear wave propagation, SIAM J. Sci. Comput. 17 (1996) 328346.