Quadrature and Runge-Kutta formulas

Error estimation and step size adjustment procedures for marry @$ order Runge-Kutta methods fail dramaticaily when appiied to an equation of he form y’(x) = f(r), A practical remedy is proposed anditsUSA exemplified with he F&berg [7,8] ff~~las. The management of storage for this important example is &O discussed.

INTRODUCTION Many Runge-Kutta formulas for the solution of ahe ordinary differential equation !j’ = f(~, y) are more accurate than usual when applied to equations of the form y’ =f(~), that is, to quadrature. This is veq common MGth respect to high order formulas-one might almost sav typical----for to reduce the number of parameters defining such formulas, some of them are cammonly chosen so that the formulas reduce to standard quadrature form&s when f does not depend on y. Examples of methods of orders five and six which reduce to quadrature formulas of ordersG.Xand seven respectively are found in Luther [ 1, 2J and Luther and Konen [3]. Methods of orders five through eight (along with error estimators) which have this property are given by Fehlberg in [4]. Generally this property has been regarded as advantageous (cf. Lambert [5, p. 14% and King [Gj has even constructed formulas with the constraint deliberately imposed that they be more accurate for qua&ature. The reasons for requiring this property in generating Hunge-Kutta formulas of high order are discussed in [y, 81. @nis work wu @ally hcsewh

supported by tht Ikision of Physicd Research of the U.S. Ewg

and Development .ridministrutio~~.

APPJJELI MATKEMATZCSAND CoMPUTATlOh’ 2%161 171 (1976) @ American Elsevier Publishing Company, Inc., 19%

161

L. F. qHAMPINE

162

The difficulty with such formulas dots not appear until one considers the estimation of their error and a ‘justmcrrt elf the step size king used to attain a specified accuracy. We have constructed a simple example to show clearly the difficulty and how it originstes. If we let

we have a formula of order t\vo,

and a formula of order three.

The latter is for the purpose of estimating second order formula 1)~ est =

fj” +

1

-’

y,,

+

, =

the local truncation

error of the

~[.r*-.f,,]~

W’e set that if .f does nr7t depend on y. t 1wn thr estnnated error is Lero. According to this estimate, any step size could he successhdly used, and a code based on this formula ancl error estimator will get into serious trouble when presented a quadrature problem. Despite appearances, the estimator is doing precisely what it night, and the formula is quite a reasonable one: we shall now look at them more closrly. When dealing with a formula of order I> the local truncation error has the form

The principal error function Q;(x, y) is eatinnrtcd hcbrt>l)~~lsing the of a formula of order p -I- 1. Becawc “a

rtwlt

tj,,+ ,

Qua&&m we

and Runge-Kutta For&m

163

have

~~t=!h=-K+1”I Yk+,,-!/r&+,1-[ Yb”+,)--ljn+ll =q?(X,,,y”)AP+‘+ Cl(ltP’“,L 731~s we sktain a computable estimate of the principal part of the local truncation error (1). This is used to decide whether to accept the step to x,,, 1 or to try again with a smaller step, and also to adjust the step size after both successful and unsuccessful steps. In the case of the family of simplified Runge-Kutta

formulas of the form

which are two stage, second order formuias, it is easy to show, as in Henrici 19, p, 67 and p. 781, that

It is naturai to choose cy by

SO 3~ to make thy error coefficients small, and this is what we have done in our example. However. we see from the remaining term that ‘pjx,y)~0 if fY ~0, that is, if we he,ve a quadrature problem. In making the formula as accurate as we could, we actually made it of higher order for quadrature. Our error estimator is correctly returning the value of zero. The difficulty lies in the fundamental assumpt’on of (1) that the local truncation error is of order p + 1, when in fact for thi< special kind of problem it is of order p + 2. Two kinds of remedy for this difficulty suggest themselves. One is to derive formulas that do not have it. This is not the easier choice, because high order formulas ‘ire extremely tedious to derive and the simphfying assumptions which lead to the difficulty are very important. Furthermore, one mi@t fairly argue that it hardlv stems reasonable to expend considerable effort deriving lt?ss accurate ‘formulas. Another remedy is to spot quadrature pro&ms and to IISP a special estimator for them. Of course, for implementation in convenient software it will be necessary to devise a way to spot quadratures automatically. /\ccording to the experience of IlO, IL

164

L, F. SHAMPlNE

121 the (i’,8) order pair of formulas of Fchlbcrg [4] iire the best high order Runge-Kutta forr-mlas known, Hull and his associates [la, 131 have emphasized the danger of quadrature and describe it as a major defect in the general use of these formulas. \Ve shalt treat this pair of formulas USthe most important example of a gencrat tcchniqutb for o\,tninitl~ a practical remedy to the quadfaturc difficulty. It is generally held that Rung+Kutta methods USCmuch less storage than competing schemes, but this is not strictly true. The determining factor es mainly the order of the scheme, and high order Runge-Kutta metho& du require a lot of storage. Because of this we shall consider in some detail ho\1 to economize on storage Lvith the Fehlherg (7,s) pair white at the same time dealing with quadrature; \ve arc able to significantly reduce the storage used by the codes of [ 10, 1 I]. There are two general questions NY shnlt i\ddr
some equation has the form

at least locally. \Vhile ive shall nlPMi[jn some cxamplh3 below, they seem to be quite !inusual except among contrived test prob1erl.i. When integrating a syste,m of equations, a code uses a step size adequate for all components. Because the method is presumed more accurate for quadrature than not, the step sin,e is very likely to be determined by a non-quadrature equation and so yield more accuracy than desired for the quadrature problem. The only possiblfa source of trouble is when the quadrature equation ought to bc integrated with a smaller step size than the non-quadrature equations. All robust codes limit their step size increases to cope with principal error functior s ~(x, y) which have an isolated zero. One simply must be prepared for a peculiar correlation causing the error to be of higher order or cvcn to vanish altogether. This kind of correlation is easily exhibited with our earlitbr example. From est = 0.2% [ .& - fo] it is obvious how to construct no11quadrature problems with est = 0 at SOIW argument x,~.The point we wish to make is that any robust code has a degree of protection in the event WC make an occasional mistake in diagnosing (ltladr;\t\lre.

Qwdtnture

and Runge-Kuttu Fonnuku

165

DETECTING QUADRATCJRE For our purposes it will not be ~reccssary to give the Fehlberg (7,8j formulas in detail. a.. it is only certain qualitative properties that arc employed. The form is

where

and the error is estimated by (j, + i - y,,, 1, Fr;, & detectio:, of quadrature WC use Only the property that the formula does more thm ace evaluation at the same due of the independent variable. In po,nt of fact this is done three times, since crt = 0 = alp, a4 = i = a,, a11 = 1 = a,3. In our approach to the automatic detection of quadrature we are limited to such formulas. We note that it is extremely common to have a, = 1 for some r. Then because CC,=O for all explicit Runge-Kutta formulas, we have a doubie evaluation at .\; for all steps except the veF first. Notwithstanding this observation, we do not recommend our scheme unless there are several such points within each step. The idea is simply that if cr,= ni, then w+e must have ki= kr for any component fr of the differential equation w,hich is a function of x only. With proper programming we must actually have the same machine numbers for ki and k,, POwe have a very delicate test for any dependence on y,, . . . , y,. This test is a necessary condition for quadrature. With the Fehlberg formulas we suggest proceeding as follows. After forming k,, compare it with k,. If no paar of components are exact]\!. equal, no equation is a quadrature, and we need corisider the matter no further. Because quadrature is very unusual, proceeding in t&s way ad&, a negligible testing overhead to the code when applied to more typical prolAms. \$‘hcn ftjrming the error estimate one should carefully group the tcnns I’\’

166

L. F. SHAMPINE

By ~0 doing he will get an exactly zero component in est for each c~uildraturC’ equation. If quadrature was not precluded by the first test, WCtest for xero components in est. if one is found, we further test if, say, k,, = k,,. If all these tests we passed for some equ&ion, we have found that at three argr:n(!nts, AT,,, A-”+ h/6, and .I-,,+ k, the equation is to machine accurw,v a function only of the independent \MiM~. In SO& a CUSSWC shA1 diagnose this equation as a quadrature for this step. Clearly if some equation is a quadmture \W shidl detect it by the ~&r>nm described. What are the possibilities of mist&* in the other direction?

Problems ::crrr near quadrature would appear to 1)~’the only fertile sourc’c of mistakes via underflows in the routine for eviiluating f(x, !I). Suppose f(x, 11) = F(x)+ cg(x,lj) for “small” 6. For a formula of order p we have cst = 0 (c/r *‘+‘), If a mistake is Itrade so that the problem is believed to lie quadrature, the quadrature err‘or is always taken to be 0 (h r’+‘). Without peculiar correlation it hardly seems credible that \ve could gcht the siuric machine numbers for fix, y) and F(x) unless c = o(h). In such a case the error estimate we ought to use, est, is small compared to the quadrature estimate we do use, with the consequence that we choose an unnecessarily small, but safe, step size. Other mi:stakes might arise from cstraordinary correlation, but as we have already noted, any code is vulnerable to such accidents, and diagnosing quadrature is just another source of possibilities. Considering all factors, we believe our scheme for diagnosing quadrature for the Fehlberg (7.8) formulas to be entirely satisfactory in practice, though in principle we would prefer formulas for which this uwbvardncss simpl>, did not arise.

Having diagnosed quadrature, there arc several wavs we might estimirte the truncation error properly. The essence of the matter is that we must produce a new value ri,,+ 1 which is of higher order than IJ,,+, for a quadrature problem. There are general principles such as halving and reversal of direction which can be lrsed, but they are quite expensive. For the Fehlberg (7,8) pair !I,,+ l is normally of order 7, but is of order 8 for quadrature. In the course of a step the function f,(x) is evaluated at a total of 10 distinct points. .4ccording to the general theory of interpolatory quadrature formulas 114, p. 8.31 jve can construct a suitable formula using irriy 0 distinct points of evaluation of h(x). We arc fortunate here that no csstra function evaluations are required to form a quadrature formula of srlfficxnt accuracy, and we even have some choice:. With other formulas it might be necessary to do extra evahlations, and one might even be willing to make

Qwdmturc

and Rttnge-Kutta I;on~A.~

it37

some extra evaluations in order to u\e a fall&r cjrla&atllre form&. We have preferred to keep the cost PSlow as possible. The formulas for y,+ , and &+ I involve k, and k, through k,, which correspond to th(: arguments Q, 4 I i, J, d.i,l. %ncc WCmust retain these values anyway, we ce,rtdniy w mt to usB &em in our qrlihuture form& so as to minimize storage. We must choose at least two Of the three remaining arguments, $, i, 9, to get a sufficiently accurate formula. Ali possihilitic~ k’crc investigated and the choice & aud -,t was found to ~kici the integration coefficients of smallest magnitude; it ic also convenient for onr schemr of managing storage. the coefficients are rational numhcrs when the argurner~ts arc:, as in our situation, and we have oh!ained them in this form after. tediorls c-omptltations which we shall not &tail. The rcranltint; (‘rror estimltw is

STORAGE

FOR THE FEIILBERC

(7,8) FORMULAS

If On;: were to program the Fehlberg (7,s) form&; in the most straightf0nvard way, he lvould use 1Sn works of storage to sohe a system of n equations. They are used to hold 9, and k, through ki3 ,Lnd to provide a vector in which to form the dependent variable y,, + hC @,,ikIat each stage just prior to calling a subroutine for the ev.Juation of k, =f(x,.-t a,h, yn + 112 /3,,,k,). The structure of the formulas which permits 1he storage to he reduced arises in the zero coefficients, which we no\v list: /3, 2 for Y2 4; /jp,3 for r 3 6; &,. /3,, ,.,, /319,5, /3,: ,,,, /I,n,,i, and also u’~ for 2’G j <4: 6~~for 1< j<4; C+,. We aisrt use PI1 j=,C21:I,4,/3,,,5= 1;,o,5. Noiin~ thd after the formation of k5, the values k, and k,3are never used again, pfuii and Enright [IO] overwrite k,, and k,, on these \.ectnrs to reduce the storage to 13~. t’nfortunateiy WC have t0 retain an additiollai vector k, for the estimation of err0r if one of the components is a quadrature. Sonetheicss, hy ;unning, we of other stnlcture we fiery accompiidl this estimation and still reduce the total storage to 12n. Our scheme uws a vector for y,, and 11 vectors 9f 1l) fnrrn tiw dependent variable in working storage I-,, . . . , r1 ,. For I’= 1,. .\t this point k, \vill 11ot 1)~ Tad further, COthe vet’tor t-11 and store k, in ‘;. ~qdrature error eslirnaie. r3 is free. The qnantitv k, will l)C Xit~d I-+ for ii (j..s =(J, .i~+ = I$= (3, thr quantity k, \vill ix* Because /j12,4= j& = 0. I(‘,~= t(‘s needed oni\i in t>e fornlatio1~ of thti depc‘ndtxnt !%ri;ihie for the e~‘~~uation of k,, and k,,;, al,tl A, wiii i>t* nc&d for the
L. F. SHAMPINE

t68

the problem has no quadrature components. form and store in c2 the quarrtity

which is part of the quadrature the quantity PI1 . $4

If qurdrature

is possible, we

error estimate. We then form and stow in t:s

+ &,5h

= P13.4~4

+ P1).&

which is part of the dependent variable needed for the w~Juation of both kIl and k,,. This overwrites k, and A;,, but they are no longer needed. Since k3 and k, will not be needed further, we have two free vectors v3 and c,. Complete the formation of the dtrpendent variable for k,, in t‘4 and store the result k,, in tirr. Form the dependent variable for k,, in c\~ and store the result k,, in 03. Complete the formation of the dependent variable for k,, in c’s and store the result kr3 in Q. All relevant quantities are now computed and stored. Then, for each component, one can compute the usual error compute the quadrature error estimate if estimate, test for quadrature, required, and accumulate the scalar norm of the estimated error. EXA,LIPLES After an extensive search of the literature we have found few esamples of differential equations involving quadrature. Some authors prefer to treat their equations in autonomous form. If they do not arise in this form. one adjoins an equation for the independent variable, which is a quadrature, to achieve the form. In some cases it is convenient to adjoin a quadrature for the evaluation of terms in the differential equations. For example if one is solving a system of equations which involve the error function, one may well adjoin the equation 2 Vi+r = -exl3( \‘7,

- x2),

Yr&+ r(O) =(I

to compute the values y” + r(x) = erf(x) which are needed. A related reyuircment is to provide term s which are integrals of experimental data. Kalaba and Scott [15, 161 provide examples in which quadratures are naturally added to the system. Shanrpine and Gordon’s text [17, pp. 24-251 discusses why and how one might modify his equations to prevent numeric,al difficulties such as singularities or overflows. The technique suggested may very

Qux!mkre mad Runge-Kuttu Formulas

169

well lead to quadrature during a portion of the integration. WC shall treat such an exmnple because it show the need to detect quadrature automatically. As a first numerical example we integrated the family of problems y’= (k+

1)X$

y(O.1)

= Jo-‘k”l!

over one slcp of length jr = 0. I for k = 0.1 , . . . ,9. In every case quadrature was diagncxed. For k = 0,l , ,. . ,7 the formula is eiact and the answers are correct to machine accuracy; the estimated errors are of the same size as the true errors. The case k = 8 is interesting, for the formula is not exact, the true relative error being ahout 4.52 X IO- 7, hit the quadrattire formula is, so that the estimated error is equal to the true error to within machine accuracy. For k=Q, neither is exact, and the true relative error of about 3.39~ 10e6 is approximated by 3.26~ 10 .‘. If’e did the same computations -with y’=

jk + 1) y/x,

y(o.1) = 10-ik+”

for k=O, 1, . . . ,Q. The solu;iqns of both families are exactly the same, namely xL+l, but this second family cf equations is not quadrature. Our scheme correctly says they are not quadrature in every case. IVith methods like the Adams family nf formula and the backurard differentiation formulas, the error depends on the solution of the prob& but not the form of the equation. This is not true of Hunge-Kutta furmulas, ;ince heir errors depend OE partial derivatives of the equation, and this example il!rrstrates the point. Only for k= 0 is the Kunge-Kutta procedure even close t;, exact for this second family; e.g., with k= 7 the relative error is dljout 4.30X X0-” and is estimated hy the standard formula to he about 1.34 x lo-“. For some kinds of two point boundary value problems

L. F. SHAMPINE

170 where

can be shown to hi\ve a trniquc sohltion [18, p. 116 ff.] Since thy prOblems share non-nq+tive solutions, it may be conwnient to study the modified problem. Solving this problem by a shooting method regards the slop y’(a)=8 as an unknown, integrates the initial value problem with Y(O) and ~‘(a\ given to compute a value y (b; s), and then varies s to find a value for which yj$;s) has the specified value at x = h. We note that if a slope s is tried which causes the computed solution to go negative, the modified problem becomes a quadrature. We have solved Example 7.1 of [ 181, which

has f(x,y)=

- y- y’,if, u= - l.y(n)= I, b= + 1. y(f>)= 1, bv. a shouting

method. With y’( - 1) = -20 the numerical solution go:cs negative, so that the problem becomes a quadrature. Using the Fchlberg (7,8) formulas \vc computed ~(1; -20) with an initial step size of 10e4 and a pure absolute error tolerance of lo-’ This initial step size is filr too small, but even so, the integration only takes 10 steps. The local truncation error WPSalso computed at each step Because the initial step size was so small and because the formulas are so accurate, the true local tnm&ion error was no greater than 6 x lo-” before the problem changed to quadrature. After the problem changed, the formulas were exact and the local truncation error was of the order of a unit of roundoff. The taquation changes its character in the course of a step, aild tbc behavior of t . 9 Ir.tegrator is illmninating. .-4 1ilrgClocal !run~: ;:-Iw :*;rc’ is incurred .t .,,I~ q Y 10m4 in the first derivative, hut this has n..,thing :o 4~ with ~~:~&.i;~tre ?! - prob?i-m 411 ~rld be diagnosed qW!rdiUre 0111~if this is tr-rlr for the whoi<: step. l’he c(jrle does not diagnose thi.i ;tep as quadrature. so the tarror arises from the proper use of the usual error estimator. The reason for the large error is that the solution is not smoclth across the change (its tll;,.r 1 &. 3 derivative has a jump) and a high order formula is inappropriate. This is &,nportant in context; the important thing 1s that our quadrature detection and error estimation works perfcctlv to allow the integration to proceed normally. iIS

RU?ERE!‘X::ES

4

5 6

7 8 9 10

11

12

13

14 i5

16 1; 1x