a 
 
CENTRAL CIRCULATION BOOKSTACKS 
 
 The person charging this material is re- 
 sponsible for its renewal or its return to 
 the library from which it was borrowed 
 on or before the Latest Date stamped 
 below. You may be charged a minimum 
 fee of $75.00 for each lost book. 
 
 Theft, mutilation, end underlining of book* are reasons 
 for disciplinary action and may result In dismissal from 
 the University. 
 
 TO RENEW CALL TELEPHONE CENTER, 333-8400 
 UNIVERSITY OF ILLINOIS LIBRARY AT URBANA-CHAMPAIGN 
 
 MA V { - 
 
 DEC 14 1998 
 JAN 1 a 2Q( 
 
 fti 
 
 ktv 
 
 When renewing by phone, write new due date below 
 previous due date, L162 
 

 YK*J$ 
 
 UIUCDCS-R-78-931 
 
 UILU-ENG 78 1723 
 C00-2383-0048 
 
 STABILITY OF VARIABLE- STEP METHODS 
 FOR ORDINARY DIFFERENTIAL EQUATIONS 
 
 by 
 
 C. W. Gear 
 
 July 1978 
 
Digitized by the Internet Archive 
 in 2013 
 
 http://archive.org/details/stabilityofvaria931gear 
 
UIUCDCS-R-78-931 
 
 STABILITY OF VARIABLE-STEP METHODS 
 FOR ORDINARY DIFFERENTIAL EQUATIONS 
 
 by 
 
 C. W. Gear 
 
 July 1978 
 
 DEPARTMENT OF COMPUTER SCIENCE 
 UNIVERSITY OF ILLINOIS AT URB ANA- CHAMPAIGN 
 URBANA, ILLINOIS 61801 
 
 Supported in part by the Department of Energy under grant ENERGY/EY-76-S-02-2383, 
 
1. Introduction 
 
 In a recent report, Zlatev [4] extends the known class of 
 
 unconditionally stable variable-step multistep methods (that is, methods 
 
 which are stable for any stepsize sequence such that the ratio of the 
 
 largest to smallest steps is bounded) to include methods of the form 
 
 (1.1) y = 0L [{h}]y . 4- a.[{h}]y „ + derivative terms, 
 n i n-1 2. n-Z 
 
 By the notation a.[{h}] we simply indicate that the a. may depend on the 
 
 ratio of one or more past steps. Zlatev shows (1.1) to be stable as long as 
 
 | Oi,, | < 1, a property that will hold for any reasonable strategy for 
 
 computing a, and a_ from {h}. (Note that |a„| < 1 for a strongly stable 
 
 constant-step method.) 
 
 The proof of convergence presented in Zlatev' s report is 
 particularly interesting for two reasons: (a) it yields a proof of 
 fixed-step convergence that is intuitively more pleasing than the 
 "standard" proof, as found in Henrici [3], being more direct; and 
 (b) it permits an extension to all strongly stable methods, which, while 
 not yielding unconditional stability (shown to be untrue in general by a 
 counter example in [1]), leads to a bound on stepsize ratios which can 
 be computed for any particular formula, and hence simply implemented in 
 practical programs. 
 
 This short paper discusses an extension of Zlatev' s proof, 
 first at an intuitive level which this writer finds particularly 
 appealing, then at a technical level. Finally, the constructive 
 application of the result to the calculation of stable-step ratios 
 (the bounds b and B) is discussed. 
 

 
 2. Zlatev's Proof Technique 
 
 Consider the strongly stable, constant-step method 
 
 (2.1) y = Z a. y . + Z h3. f(y .) 
 
 n . .. 1 n-i . _ i n-i 
 
 i=l i=0 
 
 for the initial value problem 
 
 (2.2) y' = f(y, t), y(0) = y Q . 
 
 The equation for the global error, e = y - y(t ), takes the form 
 
 n n n 
 
 (2.3) e = Z a.e .+ I h3. g .e . + d , 
 
 n . , l n-i . _ i n-i n-i n 
 i=l i=0 
 
 where g. = f evaluated at a suitable point, and d is the local truncation 
 : y n 
 
 error. The "standard" proof of convergence involves bounding the rate of 
 
 growth of solutions of this difference equation which is an order h 
 
 perturbation of a difference equation with a root at one. Zlatev considers, 
 
 instead, a difference equation for c = e - e „ . This satisfies 
 
 n n n-1 
 
 the equation 
 
 k-1 k 
 
 (2.4) c =e -e ,= Z y. c .+ Z hS . g .e , + d 
 
 n n n-1 . ' i n-i . „ l n-i n-i n 
 i=l i=0 
 
 If we write p(£) = -£ + £ & . K , P(£) = ~K + Z Y.C we find 
 
 i=l 1 i=l 1 
 
 that p(£) = K P(?)/(?-D, that is, p(?) is a "deflated" form of p(0 
 
 with the unit root replaced by zero. If the method is strongly stable, 
 
 p(£) has all of its zeros strictly inside the unit circle, hence the 
 
 growth of solutions of eq. (2.4) can be bounded by a constant times the 
 
 __ i -i 
 
 largest of the added terms ( Z h$ . g . e . + d ). If d is 0(h ) 
 
 . _ i n-i n-i n n 
 i=0 
 
 and e is 0(h P ), this term is 0(h P+1 ), hence c is 0(h P ). The global 
 
n 
 
 error e is simply e n + Z (c.), and, since there are order of 1/h terms 
 n U . t 1 
 
 i=l 
 
 in the sum, e is 0(h /h) = 0(h). Of course, this argument is 
 
 circular, so there are a few technical details to work out. This is 
 
 the subject of the next section. 
 
 3. Variable-step Methods 
 
 Consider a method given by eq. (2.1) where h depends on n, and 
 the a. and $. may vary with the ratios of stepsizes. Following the 
 notation in [1], the error equation (2.3) may be written in the form 
 (3.1) 
 
 e =(S +hS)e n +d 
 — n n n n — n-1 — n 
 
 where e is the error in the collection of saved information 
 
 *n (e ' g ' 
 
 ^n = [ V y n- 
 
 1' V-2 
 
 , . . . , h y ' , h 
 
 n y n' n-1 'n-l 
 
 , . . . ] ) , S contains the 
 
 dependence of (2.3) on g = 3f/8y, and S depends only on the coefficients 
 
 of the method. In [1] it is shown that if S (where 
 
 m 
 
 S , n ) and I S can be bounded independently 
 
 m+1 ' ' n' ' 
 
 n 
 
 s n = s =s s ,. 
 
 m . , _ 1 n n-1 
 
 j=m+l 
 
 of n > m, then convergence can be proved for consistent methods, so that 
 bounded | S and | S is a reasonable stability condition . For the 
 
 m ' ' n 
 multistep method above, 
 

 
 s = 
 n 
 
 \ a 2 ... a k j3 1 6 2 ... 
 
 K 
 
 10... 
 
 1... ( ) 
 
 
 10 
 
 _ _ 1 
 
 
 r 
 
 
 
 
 
 1 ... 
 
 O |° — 
 
 
 
 
 • 
 
 
 
 1 0_ 
 
 We now consider the similarity transformation of S to Q S Q where 
 
 n n 
 
 Q - 
 
 l o ... 
 -l l ... 
 
 O 
 
 -l i 
 
 o 
 
 o 
 
 This corresponds to a change of variables, from e,e ,,e O ,...to 
 
 n n-1 n-z 
 
 e , -c , -c ,,... and is in the spirit of Zlatev's technique. We find that 
 n n n-1 
 
QS n Q 
 
 -1 
 
 1 
 
 " Y l " Y 2 ••' " Y k-1 
 
 3, 3« ... 
 
 \ 
 
 
 
 Y l Y 2 •'• Y k-1 
 
 -3 , -3 _ ... 
 
 A 
 
 
 
 10 
 1 
 
 L 
 
 O 
 
 
 1 
 
 ... 
 
 
 
 o 
 
 1 
 1 
 
 
 
 
 
 
 
 1 0_ 
 
 n 
 
n 
 
 s n q 1 = n (Q s. q" 1 ) 
 
 111 Li J 
 
 j=m+l 
 
 1 ! y ! 
 
 n,m 
 — t , 
 
 ^n 1 
 I S B 
 
 m n,m 
 
 C^) \ cn ~ m 
 
 where 
 (3.2) 
 
 y = y . + x S with y 
 n , m n— 1 , m n m mm 
 
 = . 
 
 Suppose, for a moment, that S < r where r < 1. 
 
 n — 
 
 Then eq. (3.2) implies 
 
 (3.3) 
 
 n,m' ' — 1-r 
 
 max 
 m<i<n 
 
 Thus, if the y . are bounded, the matrix 
 
 i y 
 o 
 
 n,m 
 
 n 
 
 m 
 
 ,n-m 
 
 is bounded. Because C =0 for n > m + k and is trivially bounded, it 
 
 then follows that B is bounded if the 3. are bounded. Thus, it remains 
 
 n,m x 
 
 to consider conditions under which Is < r < 1. Note that S is the 
 
 ii n i i _ n 
 
 companion matrix of the polynomial p(£)/(£-l), all °f whose zeros are 
 inside the unit circle. Consider the coefficients for eq. (2.1) with 
 constant steps. If the largest zero of p (£)/(?-!) is r, there exists a 
 
norm | | • | | such that 
 
 S| | < r + E = r 
 
 for any £ > where S is S for constant stepsizes. (Let H be a similarity 
 
 transformation S to Jordan form with its nonzero off-diagonal elements 
 
 equal to e. Let |x| I be the max norm of Hx and si I the matrix norm 
 
 n H 
 
 I. /\ . . 
 S TT . If the 
 1 n ' ' H 
 
 coefficients a. of a variable-step method are continuous functions of the 
 1 
 
 |S is a continuous 
 n H 
 
 function of the step ratios. When all the ratios are unity, 
 
 |§ | | < r < 1. Hence, for an r such that r < r < 1, there exists a pair 
 n H 
 
 of numbers b < 1 < B such that if the step ratios lie between b and B 
 
 |S < r < 1. 
 ii u i i _ 
 
 4. Computing b and B 
 
 This can be done once and for all for any given method, it 
 does not have to be computed during an integration. First, the similarity 
 
 transform H of S to Jordan form must be computed. (If the zeros of p(£) 
 
 i i ^ ii 
 are distinct, the Jordan form is diagonal.) Next, the function |S | |„, 
 
 which depends on the step ratios, must be formed. For any pair of values 
 
 b and B. the maximum of S I L T in the hypercube of step ratios between b 
 / ' ' n' 'H 
 
 and B can be computed (approximately) . The value of b and B can 
 
 be decreased and increased respectively until this maximum 
 
 is any selected r such that 1 > r > f. (If r is chosen closer to 1, the 
 
 range between b and B will be larger, but the error constants will be 
 
 greater because of the 1/(1 - r) factor in eq. (3.3).) For many methods 
 
 one will find that b can be made very close to zero without problem as 
 
 rapid step reductions in variable step methods yield methods that are 
 
essentially one-step methods of order one or two for which the y. are 
 zero. This is fortunate because it is desirable to be able to reduce 
 the step suddenly when the solution exhibits a sharp change. 
 
 Variable Order Methods 
 
 In [2] it was shown that changing formulas "occassionally" 
 
 did not affect stability. Precisely, it was shown that, for any set of 
 
 strongly stable formulas {F. } there exists integers N. . such that if at 
 
 1 ji 
 
 least N. . steps of formula F. are taken after a switch from formula F., 
 ji i J 
 
 then the combined method is stable. This result carries over to the 
 
 A i 
 current situation with a much simpler proof. Suppose S is the matrix 
 
 S corresponding to formula F. at constant stepsize. Let H. be such that 
 
 Is 1 !! < r. < 1. Choose the r., b. and B. such that Is 1 !! < r. < 1 
 1 ' ' 'H. — i 11 i M n"H, - i 
 
 i i 
 
 for all step ratios in the range [b., B.]. Now choose r such that 
 
 max r. < r < 1. (We only consider finite sets of formula in practice. 
 
 i 
 
 The mathematically minded must add requirements such as uniformly 
 
 strongly stable to guarantee that sup r. < 1.) Now consider the max norm 
 
 i 
 
 | # | of S where all steps are taken with method i. 
 ■ i i i m 
 
 I |c n | I < I liT 1 ! I I |q n ! I I | W I I 
 
 i i m i i _ i i x i i i ' m ' 'h. ' ' i' ' 
 
 l 
 
 MhT 1 !! Mh.II r n_m 
 
 — ' ' i ' ' ' ' 1' ' i 
 
 < UIh: 1 !! ||h.|| A) n " m ] r n_m 
 
 — ' ■ i ' ' ' ' i 1 ' r 
 
 Since we have chosen r. < r, there exists an integer N. such that 
 
 i i 
 
 _i r • 
 
 l|H. I I Ih.I (— ) P < 1 for all p > N 
 
10 
 
 Thus, if at least N. steps are taken with formula F., we can bound 
 
 1 1 
 
 Is for a sequence of variable-step, variable-order formulas by 
 
 m 
 
 Kr where K = max(||H. I Ih.II). This can easily seem to be 
 
 1 ' i ' ' ' ' l 1 ' 
 
 i i ^n i i 
 sufficient to show that S is bounded as before. (The proof above 
 
 1 ' m 1 ' 
 
 yields N.. = N. for all j . In practice, a smaller value of N.. can be 
 7 ji i Ji 
 
 obtained by setting 
 
 N. . = 
 ji 
 
 log llHj h, 1 ! 
 
 log r - log r t 
 
 and 
 
 K- maxCllH" 1 !! ||H ||). 
 i,j 
 
 Only the presentation of the proof is complicated, not the principle.) 
 Consequently, a practical step and order control restriction 
 which is guaranteed to be stable for strongly stable methods is as follows: 
 
 1. If the order was changed within too few steps, stay at 
 the current order; otherwise, consider changing the order 
 in the next step. 
 
 2. Compute the "optimal" step and order for the next step 
 using your favorite algorithm. 
 
 3. If the step increase is greater than permitted, use 
 the maximum allowed. 
 
 4. If the step decrease is larger than allowed, change to 
 a one-step method to "restart"; otherwise, take the 
 recommended step. 
 
 (The question of restarting is the subject of a paper in preparation by 
 this writer. It is possible to restart efficiently at high order using 
 a one-step method which is automatically stable.) 
 
11 
 
 BIBLIOGRAPHY 
 
 1. Gear, C. W. and K. W. Tu, "The effect of variable mesh size on the 
 
 stability of multistep methods," SIAM Journal on Numerical 
 Analysis , 11 (5), October 1974, 1025-1043. 
 
 2. Gear, C. W. and D. S. Watanabe, "Stability and convergence of 
 
 variable order multistep methods," SIAM Journal on Numerical 
 Analysis , 11 (5), October 1974, 1044-1058. 
 
 3. Henrici, P., DISCRETE VARIABLE METHODS IN ORDINARY DIFFERENTIAL 
 
 EQUATIONS, John Wiley & Sons, Inc., New York, 1962. 
 
 4. Zlatev, Z., "Stability properties of variable step-size, variable 
 
 formula methods," Mathematics Institute, University of 
 Copenhagen, Preprint Series //38, 1977. 
 
„mAEC-427 U.S. ATOMIC ENERGY COMMISSION 
 
 SfSm UNIVERSITY-TYPE CONTRACTOR'S RECOMMENDATION FOR 
 
 DISPOSITION OF SCIENTIF C AND TECHNICAL DOCUMENT 
 
 I See Instructions on Reverse Side ) 
 
 AEC REPORT NO. 
 
 COO-2383-0048 
 
 STABILITY OF VARIABLE-STEP METHODS 
 FOR ORDINARY DIFFERENTIAL EQUATIONS 
 
 TYPE OF DOCUMENT (Check one): 
 
 [x] a. Scientific and technical report 
 [] 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) 
 
 RECOMMENDED ANNOUNCEMENT AND DISTRIBUTION (Check one): 
 
 Kl a. AEC's normal announcement and distribution procedures may be followed. 
 
 J b. Make available only within AEC and to AEC contractors and other U.S. Government agencies and their contractors. 
 ] c. Make no announcement or distribution. 
 
 REASON FOR RECOMMENDED RESTRICTIONS: 
 
 SUBMITTED BY: NAME AND POSITION (Please print or type) 
 
 C. William Gear 
 Principal Investigator 
 
 Organization 
 
 Department of Computer Science 
 
 University of Illinois at Urbana-Champaign 
 
 Urbana, Illinois 61801 
 
 Signature 
 
 Date 
 
 July 18, 1978 
 
 FOR AEC USE ONLY 
 
 AEC CONTRACT ADMINISTRATOR'S COMMENTS, IF ANY, ON ABOVE ANNOUNCEMENT AND DISTRIBUTION 
 RECOMMENDATION: 
 
 PATENT CLEARANCE: 
 
 Lj a. AEC patent clearance has been granted by responsible AEC patent group. 
 LJ b. Report has been sent to responsible AEC patent group for clearance. 
 LJ c. Patent clearance not required. 
 
BLIOGRAPHIC DATA 
 IEET 
 
 1. Report No. 
 
 UIUCDCS-R-78-931 
 
 2. 
 
 3. Recipient's Accession No. 
 
 Title and Subtitle 
 
 STABILITY OF VARIABLE-STEP METHODS 
 
 5. Report Date 
 
 Julv 18, 1978 
 
 FOR ORDINARY DIFFERENTIAL EQUATIONS 
 
 6. 
 
 Author(s) 
 
 C. William Gear 
 
 8. Performing Organization Rept. 
 
 No -UIUCDCS-R-78-931 
 
 Performing Organization Name and Address 
 
 Department of Computer Science 
 
 University of Illinois at Urbana-Champaign 
 
 Urbana, Illinois 61801 
 
 10. Project/Task/Work Unit No. 
 
 COO-2383-0048 
 
 11. Contract/Grant No. 
 
 ENERGY/EY-76-S-02-2383 
 
 , Sponsoring Organization Name and Address 
 
 Department of Energy 
 9800 South Cass Avenue 
 Argonne, Illinois 60439 
 
 13. Type of Report Si Period 
 Covered 
 
 Research 
 
 14. 
 
 . Supplementary Notes 
 
 . Abstracts 
 
 It is proved that, for any multistep formul 
 
 constant stepsize, there exist constants b • 
 
 of adjacent steps satisfies b < h /h , < B 
 
 — n n-1 — 
 
 of the formula is stable. A practical step- 
 
 a which is strongly stable at 
 < 1 < B such that if the ratio 
 , the variable-step implementation 
 -order control restriction which 
 
 guarantees stability is given. 
 
 Key Words and Document Analysis. 17a. Descriptors 
 
 variable-step 
 
 ordinary differential equations 
 
 stability 
 
 multistep methods 
 
 variable-order 
 
 • Identif iers/Open-Ended Terms 
 
 • COSATI Fie Id /Group 
 tAvT 
 
 Availability Statement 
 
 unlimited 
 
 19. Security Class (This 
 Report) 
 
 UNCLASSIFIED 
 
 20. Security Class (This 
 Page 
 
 UNCLASSIFIED 
 
 21. No. of Pages 
 13 
 
 22. Price 
 
 CM NTIS-15 ( 10-70) 
 
 USCOMM-DC 40329-P7 1 
 
SEP 2 1 t978 
 
\m