Wmtil kSm Wmi m KH&BRrSBESwi Hbb LIBRARY OF THE UNIVERSITY OF ILLINOIS AT URBANA-CHAMPAIGN 5/0.84 no. S9S-&&0 cof>2. 30 5/c *r r ./ 'UIUCDCS-R-T3-598 ytyu^l COO-2383-0001 1 ASYMPTOTIC ESTIMATION OF ERRORS AND DERIVATIVES FOR THE NUMERICAL SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS by C. W. Gear October 1973 DEPARTMENT OF COMPUTER SCIENCE UNIVERSITY OF ILLINOIS AT URBANA-CHAMPAIGN URBANA, ILLINOIS JHE LIBRARY. OF THE Mftti Digitized by the Internet Archive in 2013 http://archive.org/details/asymptoticestima598gear UIUCDCS-R-73-598 ASYMPTOTIC ESTIMATION OF ERRORS AND DERIVATIVES FOR THE NUMERICAL SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS* by C. W. Gear October 1973 DEPARTMENT OF COMPUTER SCIENCE UNIVERSITY OF ILLINOIS AT URBANA-CHAMPAIGN URBANA, ILLINOIS 6l801 * Supported in part by contract US AEC AT( 11-1)2 383 PREPRINT FOR LIMITED DISTRIBUTION. SUBMITTED TO IFIPs "jk. ABSTRACT Shampine [h] points out that the Milne device for estimating errors may give asymptotically incorrect results. This paper examines general classes of error estimates and shows how they can he used to form asymptotically correct estimates of the local error in the numerical solution of ordinary differential equations, or to estimate the derivatives of the solution of the differential equation. 1. INTRODUCTION Shampine [h] gives the following example to show that Milne's device for estimating the error is not asymptotically correct in some cases. Predictor: P n+2 = y n + 2h f(y Q+1 ) (D Corrector: Y ^ 2 = y n+1 + |(f(7 n+1 ) + f(y n+2 >] (2) (Asymptotically it does not matter whether we solve the corrector in eq. (2) exactly or iterate it one or more times.) If eqs. (l) and (2) are applied to y' = y, y(0) = 1 using the exact initial conditions for a fixed step size of y = 0, y = e , we find that p 2 -y 2 = -3jh 3 + o(h U ) (3) P m - y m = - \ h 3 + 0(h U ) m >_ 3 (U) The local truncation error in eqs. (l) and (2) is defined as -y(t n+ 2 } + yU n } + 2h y ' (t n + l } = ' 3 h V 3) (t n ) + 0(h U ) and -y(t n+2 ! + y(t n+1 ) ♦ |(y °° where the e.(t) have continuous (M+l-j)-th derivatives and satisfy J differential equations of the form el = f e. + f (10) j y j j (in eq. (10), as in all future equations, explicit dependence of functions on t is not shown unless necessary for clarity. Unless otherwise stated, f and its derivatives are evaluated on the solution y(t). Thus, f = 9f/9y| ,,\. The subscript n will y y=y ( t J refer to evaluation at the point (t ,y ) , while the superscript n n n will refer to evaluation at the point (t ,y(t )). Thus, y = y(t ).) n n n We will assume that eqs . (9) and (10) hold. cj> is called the principal error function (see Henri ci [3], particularly Theorems 1.5> h. 3, and 5. 12) since the local truncation error of a method is p+1 / P+2 N / \ h <{> + 0(h ). In general, .(t) is a function of derivatives of y and f and of e. for i < j , so that eq. (10) serves to define the e. i J recursively. In a multistep method in which sufficient corrector iterations are used or the order of the predictor is high enough, 4> = C y^ P+ . (See Henrici [3], Section 5-3-7.) The popular Milne method for estimating local errors (Henrici [3], Section 5-3-6) has been seen to be equivalent to forming a linear combination of values of y . and computed derivatives f(y .). For n+i n+i this reason we wish to examine the linear operator D, [z] = E [a. z(t+ih) + b.hz' (t+ih) ] (ll) h i=0 x (When we consider unequal steps, we must consider D, [z] = I [a. z(t ..) + b. hz'(t )] (lla) h n . „ in n+i in n+i i=0 where a. and b. may be functions of h.) in in For sufficiently differentiable z, we can write Djz] = E d,z (j) h j + OCh^ 1 ) (12) h j=q+l J (For unequal steps, we get M / » D.[z] = E d. z (J V + 0(h Mfl ) (12a) h n J-q+1 Jn where d. may depend on h, provided that a. and b. are sufficiently jn in in well behaved in a neighborhood of h = 0. This we must assume.) In fact, when we estimate the error of a numerical solution, we cannot evaluate eq. (ll) directly as we cannot assign a meaning to y'. Instead we must evaluate the associated nonlinear operator n k D [z] = Z [a. z(t+ih) + b.hf ( z+ih) ] (13) i=0 X X We are interested in the asymptotic expansion of D. [z] when z is identified with y for fixed t . We will call this D, [y ]. Note that n n h n D h [y n ] = D h [y n ] (1U) but that D, [y ] is not defined, h n 3. ASYMPTOTIC ERROR AND DERIVATIVE ESTIMATES In this section we will expand D [y ] and relate it to D [y ] and 4> . This result will show how y or 4 can be estimated. P P Substituting eq. (9) into eq. (13) we get k M M Djy]= E [a.{y n+1 + I e n+ V } + b.hf {y n+1 + E e n+ V}] h n i=o 1 j-p J 1 3-p J + 0(h Mfl ) (15) We write the Taylor's expansion of f(y ) about f(y ) as M f(y n ) = f(y n + E eV) + 0(h M+1 ) J=P = f n + Z f n eV + E rV + 0(h Mfl ) (l6) j=P y J j=P+l J where r. is that part of the coefficient of h which involves second and higher derivatives of f with respect to y. Note that r. depends also on e for k < j . Substituting eq. (l6) into eq. (15) and defining K. r = we get P W - \ [( v n+i + v^ 1 ' 1=0 M + E (a.e n+i + b.h{f n+i e n+i + r n+i })hJ] + 0(h Mfl ) (l T ) J=p 1 J x y j j We substitute f e. = e! - d> . from eq. (10) into this to get y J J J M L\[y ] = D ry n ] + Z D Je n ]h J h n h . h j J=P M-i k . - E E b.?fV +1 + 0(h Mfl ) (18) j=P i=0 X J where > . =

. to be expanded by Taylor's J series about t as n M-i-1 (t ) h J+1 ? +i = I hJ +l+1 J* JLsl ♦ oth* 1 ) (20- t=0 J I. where t . - t T = -=-= — (= i for equal intervals) i h When eq. (20) is substituted into eq. (l8) we get the desired result: M W - V^> + , E D h [e j ! h ° - V h J +1 ¥ B, ? (i J n ♦ 0(h W1 ) (21) where k b.(T.) £ \ ■ = -TT- (22) i=0 From eq. (12) and the continuity assumptions, D.[e n ]h J = Odi 1 * 11 ^ 1 ' J+q+l) ) (23) Eqs. (21) and (22) allow us to make the following observations: Observation 1 If (i) M > q < p or (ii) M > q > p and B_ = B =...= B =0 — 1 q-p then D h [y n ] = d q+1 h^ 1 y (q+l)n + 0(h <1+2 ) (2k) Observation 2 If q = p < M, then f> h [y n ] = (d y (<1+1)n - B Q fy**-* 1 + o(h* +2 ) (25) Observation 3 If q > p < M and B 4 , then W - -\ *p * p+1 + °<» p+2 > k. APPLICATIONS The result in the previous section will be applied to justify two common techniques and one new one for error estimation in this section. Note that the numerical solution is still assumed to satisfy eq. (9) Lemma 1 Differences can be applied either to the solution y or to the computed derivatives f(y ) to estimate y for q < N. Proof This follows directly from observation 1. If the differences of y are formed, then all b. and hence B n are zero. If differences n 1 £ of f(y ) are formed, it is trivial to verify that B = for n >v <_ l '-_ q-l. Since p >_ 1 for eq. (9) to be meaningful, one of the alternatives of observation 1 applies. Q.E.D. This means, for example, that the common technique of estimating the error in Adam's method is valid since backward differences of f(y ) are computed. These can be used to estimate the error in any order method, not just in the order used to compute the y . A simple extension can also be used to justify the following technique for fixed steps: Suppose we have an estimator D, [y ] for h y which h n satisfies the conditions of observation 1. Then D n [y ] - D. [y . ] h n h n-1 q+2 (q+2) is an estimator for h y ifq< M-l. This follows by considering the operator D [• ] - D n [• , ] = fi, [y ]. If the coefficients hn hn-1 hn B„ correspond to f) and B„ to D n , then B n = 0. If B . = for I h I h j <_ j £ q-p, then B. = for 1 <_ j <_ q-p+1 . This technique is J frequently used to estimate the error of a multistep method of order one larger than that currently in use. D [y ] is the predictor- corrector difference (see Lemma 2 below) and it is differenced again. 10 Lemma 2 If a P(EC) v predictor-corrector method is used and the predictor and corrector have the same order p, then the difference between the predicted and corrected value is a multiple of h ' y P + 0(h ) . Proof If the predictor is p n + k= . E <«i Vi * h »i f(y n + i" < 27) 1=0 and the corrector is U = .\ (a i Vi + hB i f{ W + hB k f( W (28 » 1=0 then D, [y ] = p - y h n n+k n+k = E a. y + hb. f(y ) . _ i n+i l n+i i=0 where a. = a. - a* and b. = 3. -3* (assuming 3 n = and a, = a * = -l) . 111111 k kk If the local errors in the predictor and corrector are defined as E [a. y(t ..) + h3. y'(t ..)] = c ._ h P+1 y ( P +l)n + (h P+2 ) (29) . ~ i n+i l n+i p+1 i=0 and k E i=0 then E [a* y(t n+ .) + h3* 7 f (t n+1 )] = c* +± h P+1 y (p+l)n + 0(h P+2 ; (30) D h ty] = (c p+1 - c* +1 ) h P+1 y (p+l) + 0(h P+2 ) (3D In this case (see Henrici [3], eq. (5-185)) JB*L * = k (32) p E 3* i=0 X so from observation 2, 11 k-1 £ 6. i=0 X p y = fi [y ] = [ C ~ C * ] h P +1 y(^)» + (h^ ) ( 3 3) ^n+k n+k h n p+1 ^ p+1 i=0 i Q.E.D. Comment This result is an extension of the Henrici result mentioned in Section 1, and also indicates the correction factor to be used when k-1 k Z B. ¥ £ 3*. i=0 x i=0 x Example For the example in c 3 - -i, Mi = 2 c* = i. ZBf = 1 The Milne estimate suggests that the predictor-corrector difference should be (c 3 - c*) h 3 y (3) + 0(h U ) = ^h 3 e t + 0(h U ) which is true for p - y only. But this is for n fixed, not for t fixed. In the latter case, eq. (33) tells us that *n + 2->V2 = -J h3 ^ (3) + 0(hH) <*> This example points out another problem. Why, one should ask, should we not set n = in eq. (3*0 to get P 2 - y 2 = - |h 3 + 0(hS which has been seen to be incorrect, t is, after all, a fixed value. The answer is that the initial hypothesis implied by eq. (9) is not valid at t =0 for any M > 1, but it is valid at t > for aiy M. n n 12 This is due to the starting "errors" which arise when additional values are used to start multistep methods. The corrector eq. (2) is a special case of a multistep method — namely a one-step method. In fact, the corrector is started at t = h. It is simple to verify that for t > 0, e = — te and e_ = - rr e (2-t) if this corrector is solved exactly. However, eq. (9) does not hold at t = with these values since e (0) f 0. Viewed as a two-step method, eq. (U) has the associated stability polynomial p(0 = -£ + £. Its nonprincipal root. is zero. In the case of a general multistep method with nonzero nonprincipal roots, starting errors may continue to have an effect. Henrici [3], Section 5.3-5 discusses this and shows that the starting errors do not decay for weakly stable methods. (See eq. (5-205), ibid.) Nonprincipal roots £. introduce terms of magnitude (£.) e into y - y where e is the magnitude of the starting error. For i i / M+-l\ I £. | < 1, these are smaller than 0(h J for any M, but remain of 0(e) for | E,. | =1. Consequently, we cannot expect weakly stable methods to satisfy eq. (9) unless we have a fortunate (and unlikely) choice of starting values. Strongly stable methods will satisfy eq. (9) for any fixed t away from the starting region. Thus, in the example above, "e (t )" is actually given by * 3 ( V " " & •*" < 2 -V + k (5 2> n . (35) for £ = 0. This fails to satisfy the differentiability properties at t = n = . This problem has serious implications since it means that the estimate of the error immediately following a change in step size or order may not be valid. 13 Our final application is to methods in which the local error d> h is a complicated, combination of derivatives which cannot P reasonably be computed. In this case, we use observation 3. If a differentiation formula EL [z] is chosen for which B ^ 0, and for which q > p, we get a direct estimate of B„ 4 . ^ r Op Example Consider the Heun method k " hf(y n» k l " hf( ^„ + V This is a two-stage Runge-Kutta method. The asymptotic local truncation error is h 3 * 2 = h 3 (y (3) /12 - f y y (2) A) (37) The numerical solution satisfies eq. (9). If we consider D [z] = -z(t+2h) - Mt+h) + 5z(t) + Uhz'(t+h) + 2hz'(t) --$hV fc) + 0(h 5 ) we expect to find that V y n ] = " B *2 h3 = " 6 ^2 h3 ' (38) Consider y' = y, y(0) = 1. It follows immediately that ,2 y; = y n = (i + h + |r) n (39) and that * 2 = -y/6 Uo) If we evaluate D. [y ] for eq. (39) , we get h n Ik D h [y n ] = h 3 y n + 0(h U ) agreeing with eqs. (38) and (Uo). Example The error in any fourth order Runge-Kutta method could be estimated using D [z] = -z(t+3h) - l8z(t+2h) + 9z(t+h) + 10z(t) + 9hz(t+2h) + l8hz(t+h) + 3hz(t) = 0(h) (Ul) It gives *U = " D h [y n ]/3 ° + ° (h6) In practice, this requires values at two preceding points to estimate the error in the step just performed. This probably limits the flexibility of eq. (Ul) as an estimator since those steps have to be of constant size. The equivalent of eq. ( Ul ) could be obtained for variable steps, but eq. (9) does not hold under most circumstances, so the estimate may not be valid. 5. CONCLUDING REMARKS Little has been said about the conditions under which eq. (9) holds, as a survey of all types of methods would occupy a considerable amount of space. We will summarize the situation briefly and indicate some open areas. In the case of constant steps, eq. (9) holds for any one step method whose local truncation error has an asymptotic expansion. (This includes all one step methods known to the writer.) Eq. (9) is also valid for strongly stable P(EC)n±l and corrector only multistep methods 15 (but note the starting problems indicated, by example in the last section). P(EC) methods present problems which need to be examined because there are two different values of y saved from step to step — the final corrected value y and the penultimate value y used in the n n last evaluation of f(y ). n Variable steps present serious problems. If there exists a sufficiently smooth function 9 (t) such that < A <_ e(t) < 1 (1+2) and h = hO(t ) n n it is possible to show that the above methods satisfy eq. (9)- For one step methods this is an extension of the approach used in Henrici [3] and Gear [l], Section U.7. For multistep methods it is necessary to verify that Stetter's [5] results can be applied. Gear and Tu [2] show that eq. (U2) is sufficient for stability for "reasonable" variable step multistep methods. These are multistep methods whose coefficients can be expressed via eq. (U2) as functions of h which are sufficiently smooth in a neighborhood of the origin. The analysis seems to have little meaning unless it is possible to express all step sizes h as smooth functions of h, and it is questionable what this means in the context of a computer program! 16 REFERENCES [l ] C. W. Gear, Numerical initial value problems in ordinary differential equations, Prentice-Hall, New Jersey, 1971. [2] C. W. Gear, K.-W. Tu, The effect of variable mesh size on the stability of multistep methods, UIUCDCS-R-73-570, University of Illinois, April 1973 (to appear SIAM Journal on Numerical Analysis). [3] P. Henri ci, Discrete variable methods in ordinary differential equations, John Wiley & Sons, New York, 1962. [h] L. F. Shampine, Local extrapolation in the solution of ordinary differential equations, Math. Comp., vol. 27, no. 121, January 1973, 91-97. [5] H. J. Stetter, Asymptotic expansions for the error of discretization algorithms for non-linear functional equations, Num. Math. , vol. 7, 1965, 18-31. Form AEC-427 (6/68) AECM 3201 U.S. ATOMIC ENERGY COMMISSION UNIVERSITY-TYPE CONTRACTOR'S RECOMMENDATION FOR DISPOSITION OF SCIENTIF'C AND TECHNICAL DOCUMENT I AEC REPORT NO. COO -2383-0001 i See Instructions on Reverse Side ) 2. TITLE ASYMPTOTIC ESTIMATION OF ERRORS AND DERIVATIVES FOR THE NUMERICAL SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS 3. TYPE OF DOCUMENT (Check one): 03 a- Scientific and technical report LJ b. Conference paper not to be published in a journal: Title of conference Date of conference Exact location of conference Sponsoring organization □ c. Other (Specify) 4. RECOMMENDED ANNOUNCEMENT AND DISTRIBUTION (Check one): [J a. AEC's normal announcement and distribution procedures may be followed. D b. Make available on.y within AEC and to AEC contractors and other U.S. Government agencies and their contractors. I I c - Make no announcement or distrubution. 5. REASON FOR RECOMMENDED RESTRICTIONS: 6. SUBMITTED BY NAME AND POSITION (Please print or type) C W. Gear Professor and Principal Investigator Organization Signature Department of Computer Science University of Illinois Urbana, Illinois 6l801 ^hM^O^fu^-. Date October 1973 FOR AEC USE ONLY AEC CONTRACT ADMINISTRATOR'S COMMENTS. IF ANY, ON ABOVE ANNOUNCEME RECOMMENDATION: NT AND DISTRIBUTION PATENT CLEARANCE: D a. AEC patent clearance has been granted by responsible AEC patent group. J b. Report has been sent to responsible AEC patent group for clearance. LJ c. Patent clearance not required. BIBLIOGRAPHIC DATA SHEET 1. Report No. UIUCDCS-R-73-598 3. Recipient's Accession No. 4. Title and Subtitle ASYMPTOTIC ESTIMATION OF ERRORS AND DERIVATIVES FOR THE NUMERICAL SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS 5. Report Date October 1973 7. Author(s) C. W. Gear 8- Performing Organization Rept. No. 9. Performing Organization Name and Address Department of Computer Science University of Illinois Urbana, Illinois 6l801 10. Project/Task/Work Unit No. 11. Contract/Grant No. US AEC AT(ll-l)2383 12. Sponsoring Organization Name and Address US AEC Chicago Operations Office 9800 South Cass Avenue Argonne, Illinois 13. Type of Report & Period Covered Preprint: technical 14. 15. Supplementary Notes 16. Abstracts Shampine points out that the Milne device for estimating errors may give asymptotically incorrect results. This paper examines general classes of error estimates and shows how they can be used to form asymptotically correct estimates of the local error in the numerical solution of ordinary differential equations, or to estimate the derivatives of the solution of the differential equation. 7. Key Words and Document Analysis. 17a. Descriptors 'b. Identifiers/Open-Ended Terms c COSATI Field/Group Availability Statement unlimited distribution 19. Security Class (This Report) UNCLASSIFIED 20. Security Class (This Page UNCLASSIFIED 21. No. of Pages 20 22. Price 'RM NTIS-35 ( 10-70) USCOMM-DC 40329-P7 1 tf* \# m m n nil! UB ■ EBHI HI HON HBHH9 r SjWKBl BBS mHBm smKasSmMWM m SWmt mhBHHHH^SHI nBaHHSH W3B3& 39 mm