De Meyer, Hans

Van Daele, Marnix

Van Hecke, Tanja

Vanden Berghe, Guido

English

Ghent University Academic Bibliography

Journal of Computational and Applied Mathematics 125 (2000) 107{115www.elsevier.nl/locate/camExponentially tted Runge{Kutta methods(G. Vanden Berghe , H. De Meyer 1, M. Van Daele, T. Van Hecke 2Department of Applied Mathematics and Computational Science, Universiteit Gent, Krijgslaan 281-S9, B-9000 Gent,BelgiumReceived 14 June 1999; received in revised form 6 January 2000AbstractExponentially tted Runge{Kutta methods with s stages are constructed, which exactly integrate dierential initial-valueproblems whose solutions are linear combinations of functions of the form fx j exp(!x); x j exp(−!x)g, (! 2 R or iR; j=0; 1; : : : ; jmax), where 06jmax6bs=2−1c, the lower bound being related to explicit methods, the upper bound applicablefor collocation methods. Explicit methods with s 2 f2; 3; 4g belonging to that class are constructed. For these methods, astudy of the local truncation error is made, out of which follows a simple heuristic to estimate the !-value. Error and steplength control is introduced based on Richardson extrapolation ideas. Some numerical experiments show the eciencyof the introduced methods. It is shown that the same techniques can be applied to construct implicit exponentially ttedRunge{Kutta methods. c© 2000 Elsevier Science B.V. All rights reserved.MSC: 65L05; 65L06; 65L20Keywords: Runge{Kutta method; Exponential tting; Ordinary dierential equations1. IntroductionIn the last decade, a lot of research has been performed in the area of the numerical solution ofinitial value problems related to systems of rst-order ordinary dierential equations, i.e.,y0 = f(x; y); y(x0) = y0: (1)Particular tuned methods have been proposed when the solution of the above problem exhibits apronounced oscillatory character. If a good estimate of the frequency is known in advance, onecan use linear multistep methods, whereby exact integration of a given set of linearly independent( Supported by the Bilateral Scientic and Technological Cooperation 1998, Flanders-Romania, Grant BIL 98=47. Corresponding author.E-mail address: guido.vandenberghe@rug.ac.be (G.V. Berghe).1 Research Director of the Fund for Scientic Research FWO { Vlaanderen.2 Postdoctoral fellow of the Fund for Scientic Research FWO { Vlaanderen.0377-0427/00/$ - see front matter c© 2000 Elsevier Science B.V. All rights reserved.PII: S 0377-0427(00)00462-3108 G. Vanden Berghe et al. / Journal of Computational and Applied Mathematics 125 (2000) 107{115functions is achieved. A rst good theoretical foundation of this technique was given by Gautschi[2] and Lyche [5]. Since then, a lot of exponentially tted linear multistep methods have beenconstructed; most of them were developed for second-order dierential equations where the rstderivative is absent, and applied to solve equations of the Schrodinger type. Also for rst-orderequations special tuned algorithms have been constructed. For an exhaustive list of references werefer to [4] and references cited therein. The study of exponentially tted Runge{Kutta (EFRK)methods is limited and of a very recent date. Paternoster [6] used the linear stage representation ofa Runge{Kutta method given in Albrecht’s approach and derived some examples of implicit Runge{Kutta{(Nystrom) methods of low algebraic order (for the denition of that property see [6]). On theother hand, Simos [7,8] constructed an explicit Runge{Kutta{(Nystrom) method of algebraic order4, which integrates certain particular rst-order initial-value problems with periodic or exponentialsolutions. In the present paper a general technique for the construction of exponentially tted Runge{Kutta methods is introduced.An s-stage Runge{Kutta method can integrate exactly at most a set of s linearly independent func-tions. This maximum is reached if the method is a collocation method. In the latter case, classicallythe set f1; x; x2; : : : ; xsg is considered; in our approach we choose the set fxm exp(!x); xm exp(−!x)gor equivalently for ! = i the set fxm sin(x); xm cos(x)g with m = 0; 1; : : : ; bs=2 − 1c. For explicitEFRK methods only two functions are present in the basis, i.e. m = 0; this is due to the fact thatclassical explicit Runge{Kutta methods have a stage order of at most 1.In order to realize this goal we dene an EFRK method asyn+1 = yn + hsXi=1bif(xn + cih; Yi) (2)withYi = iyn + hsXj=1aijf(xn + cjh; Yj) (3)or in tableau formc1 1 a11 a12 : : : a1; sc2 2 a21 a22 : : : a2; s: : :cs s as1 as2 : : : assb1 b2 : : : bs(4)This formulation is an extension of the denition of a classical Butcher tableau. A motivationfor the introduction of i-values is given in [10]. In this paper we dene that a function g(x) isintegrated exactly by a (EF)RK method if yn=g(xn) for all problems whose solution is g(x). If eachstage equation and the nal step equation exactly integrate g(x) then the (EF)RK scheme itself alsointegrates g(x) exactly. Indeed, for small h the uniqueness of the solution of system (3) is guaranteedby a proof similar to that of Theorem 341 of [1], while on the other hand this unique solution isYi = g(xn+ cih) by construction, given that each equation in (2) and (3) exactly integrates g(x) andthat the solution of the problem to be solved is g(x). The above conditions give rise to a system ofequations for the components of b and A of the following form (we present the exponential formG. Vanden Berghe et al. / Journal of Computational and Applied Mathematics 125 (2000) 107{115 109and v=!h, ! 2 R). The number of equations to retain depends on the number 2(m+1) of elementspresent in the basic set.exp(v)− 1 vsXi=1bi exp(civ) = 0;exp(v)−sXi=1bi(1 vci) exp(civ) = 0;...(5)together withexp(civ)− i  vsXj=1aij exp(cjv) = 0;ci exp(civ)−sXj=1aij(1 vcj) exp(cjv) = 0;...(6)with i = 1; : : : ; s.2. Explicit EFRK methodsIn the case of explicit EFRK methods the number of equations in (5) and (6) is restricted to two.Below we give some examples.For the two-stage method we have chosen c0 = 0 and c1 = 1. Eqs. (5) and (6) fully deneb1; b2; 1; 2 and a21 as (here we always present the trigonometric solution, i.e. != i; v= h):0 11 cos(v) sin(v)=vsin(v)v[cos(v) + 1]sin(v)v[cos(v) + 1](7)This EFRK method reduces in the limit v! 0 to the classical scheme with a21 = 1 and b1 = b2 = 12 ,which has algebraic order 2.For the three-stage explicit method we have chosen c1 = 0; c2 = 13 and c3 =23 . Eqs. (5) leaveone degree of freedom. In order to fully dene all components of b we have added the followingsupplementary equation:b1 + b2 + b3 = 1; (8)110 G. Vanden Berghe et al. / Journal of Computational and Applied Mathematics 125 (2000) 107{115which expresses that (2) reproduces exact results whenever f(x; y) = 1. Eqs. (6) again leave onedegree of freedom; we have chosen a31 = 0. The following exponentially tted method emerges:0 113cos(v=3)sin(v=3)v231 0sin(2v=3)v cos(v=3)b1 b2 b3(9)withb1 =−4 cos(v=3) sin(v=3) + v cos(v=3)− 2 sin(v=3) + v2v sin2(v=3);b2 =−2 cos(v=3) sin(v=3) + v cos(v=3)− sin(v=3)v(cos(v=3)− 1) ;b3 =−4 cos(v=3)2 sin(v=3)− 2 cos(v=3) sin(v=3) + v cos(v=3) + vv sin2(v=3):This method reduces for v! 0 to the classical scheme with a21= 13 ; a32= 23 and b1= 14 ; b2=0; b3= 34which has algebraic order 3.For the four-stage method two additional equations can be added to (5) [10]. We require that therst equation (2) reproduces exact results whenever f(x; y) is 1 or x. This delivers the followingtwo equations for the b components:4Xi=1bi = 1;4Xi=1bici = 12 :(10)Eqs. (5) and (10) fully determine the components of the b-vector. Eqs. (6) leave three degrees offreedom. Inspired by the classical fourth-order scheme we choose a31 = 0 and a42 = 0; additionally,we also set 4 = 1. The following tableau emerges:0 11=2 cos(v=2)sin(v=2)v1=21cos(v=2)0tan(v=2)v1 1 0 02 sin(v=2)vb1 b2 b3 b4(11)G. Vanden Berghe et al. / Journal of Computational and Applied Mathematics 125 (2000) 107{115 111withb1 = b4 =2 sin(v=2)− v2v(cos(v=2)− 1) ; b2 = b3 =v cos(v=2)− 2 sin(v=2)2v(cos(v=2)− 1) :For v ! 0 this method reduces to the well-known classical scheme with a21 = 12 ; a32 = 12 ; a43 = 1and b1 = b4 = 16 ; b2 = b3 =13 , which has algebraic order 4.In order to have an idea of the form of the local truncation error (LTE) for the explicit EFRKmethod, we have calculated for the three above-mentioned schemes the dierence y(xn+1) − yn+1with y(xn+1) a notation for the exact solution in the point xn+1. Note that the well-known theory ofButcher based on rooted trees cannot be applied here, since the row-sum rule is not fullled. Wepresent the LTE for a non-autonomous scalar equation. All occurring partial derivatives of f(x; y(x))with respect to x and y are collected in order to express these LTEs in terms of total derivatives ofy(x). The following results were obtained. For the two-stage method (7)LTE =− h312[y(3) − 3fyy(2) + 2(y0 − 3fyy)] + O(h4): (12) For the three-stage method (9)LTE =h4216[y(4) + 2fyy(3) + 6f2yy(2) + 2(y(2) + 2fyy0 + 6f2yy)] + O(h5): (13) For the four-stage method (11)LTE=− h52880[y(5) − 5fyy(4) − 10fxyy(3) − 10fyyfy(3) + 10f2yy(3)+15fyy(y(2))2 + 30fxyfyy(2) + 30fyyffyy(2) − 30f3yy(2)+2(y(3) − 5fyy(2) − 10fxyy0 − 10fyyfy0+10f2yy0 + 15fyyyy(2) + 30fxyfyy + 30fyyffyy − 30f3yy)+152((fy + fyyy)y(2) + 2(fy + fyyy)y)] + O(h6): (14)All functions are evaluated at x = xn and y = yn. Note that the leading order terms in the givenLTEs become zero for y = sin(x) or y = cos(x). Moreover in each of the LTE-expressions thederivative of lowest order occurring for = 0 is y(2), showing that classically only the set f1; xg isintegrated exactly.3. Local error estimation and a good choice for There exists no mathematical theory to determine  in an exact way. The only goal we can putforward is to make the LTE as small as possible by calculating  by means of a heuristically chosenalgorithm. Since for a scalar equation in each of the expressions for the LTEs (12){(14) the term112 G. Vanden Berghe et al. / Journal of Computational and Applied Mathematics 125 (2000) 107{115y(2) +2y is present we propose to calculate  in each integration interval [xn; xn+1] in the followingway:=s−y(2)(xn)y(xn); n= 0; : : : if y(xn) 6= 0 and = 0 otherwise:Note that if y(x) is a linear combination of sin(x) and cos(x), and y(xn) 6= 0, we obtain from thisrecipe  =  and y(x) will be integrated exactly (if innite precision arithmetic is assumed). For asystem of n rst-order equations, we propose to make the Euclidean norm of all (ty(2) + 2 ty; t =1; : : : n) present as small as possible. This results in=s−Pnt=1ty(xn) ty(2)(xn)Pnt=1ty(xn)2: (15)The expressions for the occurring second derivatives can be obtained analytically from the givenODEs or calculated numerically using previously derived y(xn−j) values. The -values used arethen in each integration interval taken as the positive square root of the numerically obtained 2. Ifnegative 2-values are obtained,  is replaced in the corresponding formulae by i (i2 =−1). In fact,in this case, the exponential functions instead of the trigonometric ones are integrated.Since the -values used are calculated, they are never exact due to nite precision arithmetics. Asa consequence, the leading term of the LTE does not vanish. This means that one can try to estimatenumerically, for a chosen , the LTE. A technique which can be used is Richardson extrapolation.First we x  by (15). We consider a Runge{Kutta method of order p 2 f2; 3; 4g to obtain thesolution yn+1 at xn+1. Under the usual localizing assumption that yn = y(xn) it follows from (12) to(14) that the LTE can be written in the form Tn+1 = y(xn+1)− yn+1 =C(y; f)hp+1 +O(hp+2); whereC(y; f) is some function of y, its derivatives, f(x; y) and its partial derivatives with respect to xand y, all evaluated at the point (xn; yn). Let us now compute a second numerical solution at xn+1by applying the same method twice with steplength h=2 and also starting in xn; denote the solutionso obtained by zn+1. By starting in both calculations at the same point xn one can work during thesetwo processes with the same value of . The error is nowTn+1 = y(xn+1)− zn+1 = 2C(y; f)(h=2)p+1 + O(hp+2):From these two estimates for the LTE one can derive that the error in the second calculation isgiven by: error = 2C(y; f)(h=2)p+1  (zn+1 − yn+1)=(2p − 1) If the user asks for a given tolerancetol, he can control the steplength and the error in the following way:if jerrorj<= tol accept the step and progress with the zn+1 valueif jerrorj> tol reject the stepThe step is then adapted in the following way:hnew = hold min(facmax;max(facmin; fac(tol=error)1=(1+p)));with facmax and facmin, respectively, the maximum and minimum acceptable increasing or decreas-ing factors. The symbol fac represents a safety factor in order to have an acceptable error in thefollowing step. The method for estimating the LTE and the notation used to dene hnew is given inHairer et al. [3]. In the code we have developed the following values for these factors were taken:facmax = 2, facmin = 0:5, fac = 0:9.G. Vanden Berghe et al. / Journal of Computational and Applied Mathematics 125 (2000) 107{115 1134. Numerical experimentsIn this section we solve some initial-value problems having as solution a combination of sine orcosine functions. The following cases have been considered.Example 1.y00 =−30 sin(30x); y(0) = 0; y0(0) = 1 (16)with y(0) = 0; y0(0) = 1; its exact solution is y(x) = sin(30x)=30. Eq. (16) has been solved in theinterval 06x610 with tol = 10−10.Example 2.y00 + y + y3 = B cos(x); y(0) = 0:20042672806900; y0(0) = 0 (17)with B=0:002, =1:01; accuracy is judged by comparison with a Galerkin approximation obtainedby Van Dooren [9]y(x) 0:200179477536 cos(x) + 0:246946143  10−3 cos(3x)+0:304014  10−6 cos(5x) + 0:374  10−9 cos(7x):Eq. (17) has been solved in the interval 06x6300 with tol = 10−8.Example 3.y00 + y = 0:001 cos(x); y(0) = 1; y0(0) = 0; (18)with exact solution y(x)=cos(x)+0:0005x sin(x). Eq. (18) has been solved in the interval 06x61000with tol = 10−8.Example 4.y00 + 0:2y0 + y = 0; y(0) = 1; y0(0) = 0; (19)with exact solution y(x) = exp(−0:1x)(sin(p0:99x)=p99 + cos(p0:99x)). Eq. (19) has been solvedin the interval 06x62 with tol = 10−8.After having rewritten the above equations as equivalent rst-order systems, we have solved themwith the new fourth-order EFRK method (11) and the equivalent classical fourth-order RK method;in both cases we applied the above-described Richardson extrapolation technique. We have calculatedat the endpoint of the integration interval for each example the Euclidean norm of the error vectorwith components dened as the dierence between the numerical and the exact values of the solutionand its rst derivative. These data are collected in Table 1 together with the number of acceptedand rejected steps.It is clear from Table 1 that for equations with a purely trigonometric solution (Example 1), thenew method is superior to the classical one. In this example, (15) initially gives a quite accuratevalue for , such that there is a small error and a huge increase in the stepsize in each step.114 G. Vanden Berghe et al. / Journal of Computational and Applied Mathematics 125 (2000) 107{115Table 1Comparison of the Euclidian norms of the endpoint global errors obtained by using the fourth-order EFRK and classicalRK-methods with step-length and order control based on Richardson extrapolation techniquesAccepted steps Rejected steps ErrorExample 1 Exponentially tted 33 16 4:01  10−9Classical 3495 5 1:01  10−8Example 2 Exponentially tted 1191 0 1:01  10−7Classical 2167 0 2:79  10−5Example 3 Exponentially tted 1946 1 1:14  10−5Classical 9826 0 5:75  10−5Example 4 Exponentially tted 46 3 1:96  10−7Classical 57 0 2:70  10−7However, these errors accumulate and at a certain moment such a large stepsize causes a rejectedstep. After decreasing the stepsize, this process is repeated. In cases where the solution is of a mixedtrigonometric form, although one of the terms dominates the solution (Examples 2 and 3), the newmethod is still more ecient than the classical one. In Example 4 where the solution is a mixtureof exponential and trigonometric functions both methods considered are practically equivalent.5. Implicit EFRK schemesFor a collocation RK scheme with s stages the set of functions f1; x; : : : ; xsg are exactly inte-grated. It is obvious that in the case of EFRK schemes the couples fxm exp(!x); xm exp(−!x)g;m=0; 1; : : : ; bs=2−1c have to be taken into account together. Two- and three-stage implicit methodshave been studied. For the two-stage case, where we have chosen 1 = 2 = 1 identical results havebeen obtained as the ones derived by Paternoster [6]. For general, three-stage methods the length ofthe resulting expressions is large; therefore, we present a rather simple case, the LobattoIIIA EFRKmethod; we have xed beforehand the ci-values, i.e. c1 = 0; c2 = 12 and c3 = 1. For these methodsas well the m = 0 as the m = 1 couples are taken into account in the set of basic functions. Thefollowing results emerge:a11 = 0; a12 = 0; a31 = 0; 1 = 0;a21 =−v2 − 4 + (4 + 2v2) cos2(v=2)2v2D; a23 =v2 − 4 + 4 cos2(v=2)2v2D;a22 =−(v2 − 4) cos(v=2) + 2v sin(v=2)− 2v cos2(v=2) sin(v=2) + 4 cos3(v=2)v2D;2 =−v cos(v=2) sin(v=2)− 2 + 2 cos2(v=2)D;withD =−v sin(v=2)− 2 cos(v=2) + 2 cos3(v=2);G. Vanden Berghe et al. / Journal of Computational and Applied Mathematics 125 (2000) 107{115 115b1 = a31 =2 sin(v=2)− v cos(v=2)v2 sin(v=2); b2 = a32 =−2 sin(v) + 2vv2 sin(v=2);b3 = a33 =2 sin(v=2)− v cos(v=2)v2 sin(v=2); 3 = 1:For ! 0 this method reduces to the classical LobattoIIIA method and has algebraic order 4.6. ConclusionsWe have developed a technique to construct explicit and implicit exponentially tted Runge{Kutta methods. Some examples of both types are given. A heuristic way of determining  is givenfor the case of explicit methods, which is based on the typical expression of the LTEs. In caseswhere the solutions are purely trigonometric the use of Richardson extrapolation is doubtful. Thelarge percentage of rejected steps in the rst example supports this statement. Better methods forerror estimation are welcome. It is shown that the introduced construction method can also bestraightforwardly applied to obtain exponentially tted implicit methods.AcknowledgementsThe authors are very grateful to two anonymous referees, whose comments greatly improved thispaper.References[1] J.C. Butcher, The Numerical Analysis of Ordinary Dierential Equations, Wiley, Chichester, UK, 1987.[2] W. Gautschi, Numerical integration of ordinary dierential equations based on trigonometric polynomials, Numer.Math. 3 (1961) 381{397.[3] E. Hairer, S.P. NHrsett, G. Wanner, Solving Ordinary Dierential Equations I, Nonsti Problems, 2nd Edition,Springer, Berlin, 1993.[4] L. Ixaru, Operations on oscillatory functions, Comput. Phys. Comm. 105 (1997) 1{19.[5] T. Lyche, Chebyshevian multistep methods for ordinary dierential equations, Numer. Math. 19 (1972) 65{75.[6] B. Paternoster, Runge{Kutta({Nystrom) methods for ODEs with periodic solutions based in trigonometricpolynomials, Appl. Numer. Math. 28 (1998) 401{412.[7] T.E. Simos, An exponentially tted Runge{Kutta method for the numerical integration of initial-value problems withperiodic or oscillating solutions, Comput. Phys. Comm. 115 (1998) 1{8.[8] T.E. Simos, E. Dimas, A.B. Sideridis, A Runge{Kutta{Nystrom method for the numerical integration of specialsecond-order periodic initial-value problems, J. Comput. Appl. Math. 51 (1994) 317{326.[9] R. Van Dooren, Stabilization of Cowell’s classical nite dierence methos for numerical integration, J. Comput.Phys. 16 (1974) 186{192.[10] G. Vanden Berghe, H. De Meyer, M. Van Daele, T. Van Hecke, Exponentially tted explicit Runge{Kutta methods,Comput. Phys. Comm. 123 (1999) 7{15.

Exponentially fitted Runge-Kutta methods

Abstract

Similar works

Full text

Available Versions

Ghent University Academic Bibliography

Crossref