Digitized by the Internet Archive 
 in 2013 
 
 http://archive.org/details/rungekuttastarte938gear 
 
si om 
 
 no *93o 
 
 t 7*v*_ 
 
 « 
 
 dec. 
 
 UIUCDCS-R-78-938 
 
 UILU-ENG 78 1733 
 C00-2383-0053 
 
 RUNGE-KUTTA STARTERS FOR MULTISTEP METHODS 
 
 by 
 C. W. Gear 
 
 September 1978 
 
UIUCDCS-R-78-938 
 
 RUNGE-KUTTA STARTERS FOR MULTISTEP METHODS 
 
 by 
 C. W. Gear 
 
 September 1978 
 
 DEPARTMENT OF COMPUTER SCIENCE 
 
 UNIVERSITY OF ILLINOIS AT URB ANA- CHAMPAIGN 
 
 URBANA, ILLINOIS 61801 
 
 Supported in part by the U. S. Department of Energy under contract 
 US ENERGY/EY-76-S-02-2383. 
 
I. INTRODUCTION 
 
 Twenty years ago discussions of the relative merits of multistep 
 versus RK methods made frequent reference to the "difficulty" of starting 
 a multistep method; that is, the cost of obtaining the additional k - 1 
 values needed to start a k-step method when just the initial value of the 
 problem is available. Some early programs for constant step, constant 
 order, multistep methods used several RK steps to start. For example, 
 a four-step Adams formula can be started with three fourth order RK 
 steps to obtain y!, i = 0, 1, 2, 3 and y from y . No error control was 
 provided in such a program, but at that time the error control in most 
 programs consisted of an external adjustment of the step by the user. 
 Nordsieck [12] designed one of the first multistep methods with automatic 
 starting, but the introduction of variable-order methods finally "solved" 
 the problem because these methods can start at whichever order corresponds 
 to a one-step method — usually first or second order. See, for example, 
 Gear [6], Hindmarsh [9], Krogh [11], Shampine and Gordon [13]. (It is 
 interesting to observe that Nordsieck' s scheme, which consisted of 
 integrating forwards and backwards over several steps several times, 
 amounts to a variable-order scheme: because Adam's methods are used and 
 because of the representation used, it can be shown that the order is 
 increased at each step in the Nordsieck starter.) 
 
 Starting at first or second order and relying on the order 
 control to increase the order to an appropriate value is programmatically 
 simple because the mechanism is already present in the program. It also 
 has the advantage of providing local error control, which was not really 
 present in the Nordsieck scheme. However, it is not very efficient for 
 any but low accuracies as very small stepsizes must be taken at low orders 
 
to maintain accuracy. This inefficiency is often offset by the improved 
 efficiency of high order multistep methods for many classes of problems 
 (see Hull et al. [10]). However, there are large classes of problems, 
 arising, for example, in simulation, in which there are frequent 
 discontinuities in first or higher derivatives. Since multistep methods 
 are based on polynomial interpolations over an interval of several steps, 
 they break down when the interval includes a discontinuity in a derivative 
 of degree lower than the polynomial. Hence, a restart is necessary. 
 One-step methods such as RK do not have this difficulty as long as a mesh 
 point falls exactly on each discontinuity, and, therefore, RK methods are 
 usually superior for these "non-smooth" problems. 
 
 This paper presents a RK-like technique for starting (or 
 restarting) any of the currently popular automatic multistep ODE integrators 
 at a high order. (An order four starter will be recommended for general 
 use.) The technique uses fewer function evaluations than are used if a 
 conventional RK method is applied repeatedly, and provides some degree of 
 error control, as well as data for estimating the step size to be used in 
 the first step of the multistep method. 
 
 Section 2 presents the idea in a non-RK framework to establish 
 the existence of methods of the desired type and to derive an upper bound 
 on the number of function evaluations needed in an explicit RK starter 
 (1 + p(p - l)/2 for a p-th order method). Section 3 examines the lower 
 bound on the number of function evaluations in explicit RK starters for 
 p <^ 4, and suggests a particular fourth order method that uses six 
 function evaluations (the lower bound). Section 4 examines implicit RK 
 starters which might be useful for some stiff problems and proves that 
 the minimum number of stages for a p-th order starter is p, and that p-stage 
 methods with desirable stability properties exist. Section 5 considers 
 
programming considerations and error control, while section 6 discusses 
 the detection of discontinuities. 
 
2. STARTING BY EXTRAPOLATION 
 
 If a k-step multistep method is to be started at p-th order, it 
 is desirable to generate k - 1 additional starting values which are p-th 
 order accurate.* Since k is p - 1 or p in most methods in use (Adams 
 or backward differentiation formulas), we will consider the case k = p. 
 If we can generate h q y ,~, , q = 0, 1,..., p, with error no greater than 
 0(h ) , we can easily compute the starting values needed for a code which 
 stores either y (approximations to y(t.)), y! , or h^y^ q /q! (the 
 Nordsieck vector). One way of doing this is by extrapolation as shown 
 below. (It should be emphasized that the following is an existence 
 
 argument. It is not the way to compute the required values!) 
 
 i 
 Let y., i = 0,..., £ be the results obtained from the application 
 
 of the forward Euler method for £ steps to the differential equation 
 
 y' = f(y, t) , y(0) = y_, using step size h = H/&. We will assume that 
 
 f(y, t) has as many derivatives as needed (for example, bounded fourth 
 
 partial derivatives will be needed for a fourth order method). Then it 
 
 is known that 
 
 (2.1) yf=y(ih)+ I h q e (ih) + 0(h P+1 ) 
 
 q=l q 
 
 (see Stetter [14]. Also, there exist sets of numerical differentiation 
 formulas for k + 1 equally spaced points of the form 
 
 (2.2) Z d.,z(ih) = h k z (k) (kh/2) + Z c h S z (s) (kh/2) + 0(h P+1 ) 
 i=0 ±tC s=k+l SlC 
 
 (In fact, c , = when s - k is even.) 
 
 * If restarting occurs a fixed number of times in a fixed interval, 
 errors of order hP in the starter do not cause a total error of more than 
 hP in the global solution, but a practical algorithm may restart more and 
 more frequently as the error tolerance is decreased, so we will try and 
 keep the restarting error to 0(hP + l). 
 
Applying eq. (2.2) to eq. (2.1) with k = p, h = h = H/p we get 
 
 (2.3) 
 
 S d. y? = h P y (p) (H/2) +0(h P+1 ) 
 1-0 ip l p p 
 
 Since the left-hand side of eq. (2.3) can be evaluated, we can compute an 
 order p accurate approximation to H y . From now on, we will omit the 
 term 0(h ) to ease the eye strain. 
 
 Now apply (2.2) to (2.1) with k = p - 1 and h = h 1 and h . 
 
 We get 
 (2.4a) 
 
 P Z d. y?" 1 = h P "^ y (P_i; (H/2) + h P , ej p - i; (H/2) + c„ 1 
 1=0 
 
 p-1 _ p-1 „(p-l) 
 'i p-1 y i p-1 
 
 (P-D 
 p-1 "I 
 
 (P) 
 
 PP P+l 
 
 (2.4b) 
 
 p-1 
 
 Ed. 1 y P = h 15 " 1 y (p_1) (H/2 - h /2) + h P e 1 (p " 1) (H/2 - h /2) 
 
 i= Q X P-1 ! P P P 1 P 
 
 + c h P y (p) (H/2 - h II) 
 PP P P 
 
 Note that h?- 1 y (p - 1) (H/2 - h ID - h^ 1 y (p " 1) (H/2) - \ h P y (p) (H/2) 
 P P P 2 p 
 
 and h P e 1 (p_1) (H/2 - h /2) = h P e n (p_1) (H/2) (neglecting terms of order h P+1 ). 
 pi p pi 
 
 Furthermore, we know h y (H/2) from eq. (2.3) to order p, hence eqs. (2.4) 
 can be rewritten as 
 
 h P-l y (P- 1 ) (H /2) + h P , e n (p_1) (H/2) = known values 
 p-1 p-1 1 
 
 and 
 
 h P-l y (P- 1 )(H/2) + h P e 1 (p_1) (H/2) = known values. 
 P P 1 
 
 These can be solved for H P_1 y (p_1) (H/2) and H P e^ P_1) (H/2) to order h P 
 
 since h ^ h 
 
 P P-1 
 
 The process continues in a similar fashion. Suppose we have 
 already computed H S y (s) (H/2) for s = k + 1, . . . , p, and H^"*" 3 ef s) (H/2) for 
 j=l,...,p-k-l, s=k+l,...,p-j. We now apply (2.2) to (2.1) 
 
with h = h. , h. ,,,... , and h . We get 
 
 K. K.+1 p 
 
 (2.5a) I d. k y k = h k y (k) (H/2) + h£ +1 e< k) (H/2) +. . .+ h£ e ( ^(H/2) 
 i=0 ^ 
 
 (2 - 5b) ,* n d i k yk+1 ■ h u + i * (k)(H/2 - h k+ i /2) + C e i k>(H/2 - \ + i /2) + - 
 
 1=0 
 
 k 
 
 E 
 i=0 
 
 (2.59) E d., y? = h k y (k) (H/2 - (p-k)h /2) + h k+1 e n (k) (H/2 - (p-k)h /2) +. 
 , rt ik i p p pi p 
 
 (k) (k) 
 As before, the values of y and e at H/2 - mh can be shifted to values 
 
 at H/2 by means of Taylors series using the already computed higher 
 
 derivatives of y and e . After this transformation, eqs. (2.5) are a set 
 
 q 
 
 (k) 
 of p + 1 - k linear equations in the P + 1 - k unknowns y (H/2), 
 
 (k) (k) 
 
 e 1 (H/2),..., e , (H/2) whose coefficient matrix is a modification of the 
 1 p-k 
 
 Vandermonde matrix, and is nonsingular since the h 's are all different. 
 Thus, this process yields order p correct values to H 4 y (H/2) for 
 q = 1, 2,..., p. The "cost" of this in function evaluations can be seen to 
 be 1 + p(p-l)/2 because f(y n » 0) is evaluated once and used in all p 
 integrations. An example is presented at the start of the next section. 
 
3. EXPLICIT RUNGE-KUTTA METHODS 
 
 The extrapolation technique of the previous section can be viewed 
 as a Runge-Kutta method. We will first illustrate this for the case p = 3. 
 Let h = H/3 and consider an autonomous system. We have 
 
 r 3 
 
 (3.1) ( 
 
 y x =y + hf(y ) =y + k Q 
 
 yJ-yJ + hfCy?) - y + % + \ 
 
 jl-7 3 2 + mj 3 2 ) =y + k o + k i + k 2 
 
 y l = y o + T f( V = y o + 2 k o 
 
 where k Q = hf (y Q ) 
 
 3 
 where k = hf (y ) 
 
 where k = hf(y_) 
 
 2 23h £/ 2v ■ 3 . , 3 , 
 
 y 2 - yj+T^) = y + 2 k o + I k 3 
 
 k. 
 
 yj - y + 3hf (y Q ) = y Q + 
 
 3k. 
 
 where k = hf (y ) 
 
 We use the difference formula given in (3.2) below where all derivatives 
 
 4 
 are evaluated at t = 3h/2. 0(h) terms are neglected. 
 
 "3 3 3 _ 3 3 (3) 
 
 D 3 - y 3 - 3y 2 + 3 Yl - y Q = h y 
 
 i n 3 3 3 3 ,,2 (2) 3 (2) 
 
 2 = y 3 " y 2 " y l y y e l 
 
 (3.2) 
 
 : n 2 2 2 3h.2 (2) ,3hv3 (2) 
 
 J D 2 = y 2 " 2y l + y = ( T } y + ( T ) e l 
 
 t 
 
 n 3 3 3 . (1) .2 (1) , .3 (1) h 3 (3) 
 
 1 y 2 " y l = y e l e 2 24 Y 
 
 (1) , 9^2 _(1) ^27^3 _(1) ^9^3 „(3) 
 
 1 
 
 D l = y 2 " y = 3hy + 2 h e l + T h e 2 + « Y 
 
 n l 1 iT , (1) Q ,2 (1) ._, 3 (1) , 27 ,3 (3) 
 
 ^- 1 = Y l " y = y e l e 2 24 y 
 
 (Symmetric differences have been used above to reduce arithmetic.) From 
 these we find that 
 
namely 
 
 (3.4) 
 
 'V y (3 > - D ] + 0(hS 
 
 (3.3) I h 2 y <« - | D 3 - | »\ + o(h 4 ) 
 
 9 n 3 4 ^2 J 1 _1 9 „3 
 
 2 D ri D l + 6 D l + lV° (h) 
 
 Fairly clearly, we can view this as a Runge-Kutta process, 
 
 JL 
 
 f k. = hf(y n + Z 3.. k. ), 
 i i ml ij j-1 
 
 i = 0,.. . , q-1 
 
 Lh S y (s) = Z k -_!T. s + 0(h P+1 ), s = 1,..., p 
 where, in the example above, q is four, the matrix of 3 coefficients is 
 
 B = 
 
 
 
 
 
 
 
 
 
 1 
 
 
 
 
 
 
 
 1 
 
 1 
 
 
 
 
 
 3 
 
 
 
 
 
 
 
 and the matrix of y coefficients is 
 
 3/2 
 
 "-3/8 -1/6 1 
 9/4 -2 
 9/8 3/2 1 
 _ -2 -4/3 0_ 
 
 This gives estimates for the derivatives at H/2 = 3h/2. Estimates of the 
 derivatives and/or function values to the same order of accuracy can be 
 obtained at any point within a constant multiple of the interval H. We 
 used the midpoint above to take advantage of symmetry, but for convenience 
 in the discussion below we will consider the problem of computing the 
 
derivatives at the origin. In that case the matrix above becomes 
 
 1 -5/3 
 
 3 
 r ° = 1 
 
 -1/3 
 Obviously, the estimate for hy (0) is simply k . 
 
 In the example above it took A function evaluations for 
 third order. Using the technique of section 2, it would take 7 function 
 evaluations for fourth order and 11 for fifth. It is well known that 
 Runge-Kutta methods of third, fourth, and fifth orders are possible with 
 3, A, and 6 function evaluations, respectively. This naturally raises the 
 question: Can the first p derivatives be found to accuracy 0(h ) with 
 fewer than 1 + p(p-l)/2 function evaluations? Note that this is not the 
 same problem as finding embedded RK methods of several orders such as 
 the RK Fehlberg method. See Bettis [1]. For p <^ 3 the answer is no, 
 but for p = A it is possible in six but no less. The cases p = 1 and 
 p = 2 are trivial. The problem is examined in detail for p = 3 below 
 and for p = A in the appendix. 
 
 3.1. Non-Existence of 3-Stage Method of Third Order 
 
 In the discussion below, all variables are evaluated at 
 t = 0, y = y , unless specified otherwise. Thus we can write 
 hy' = hf = k . For a system of m equations in y , y ,..., y 
 
 (3.5) 
 
 h 2 y" = h 2 I^f 1 
 
 i=l dy X 
 
10 
 
 We will write the right-hand side as h f-,f. Similarly 
 
 (3.6) 
 
 3 m = h 3 ™ ( jh_r l^^-^Kfh 
 
 h y 
 
 i,j=l 8y X 8y J 
 
 3y 9y J 
 
 = h 3 (f" 2 f 2 + f 2 f) (definition of f £ ) 
 
 Writing ot. = £ 3.., the k. in eqs. (3.4) can be expanded to 
 
 k = h f 
 
 ot 
 k = hf + a h 2 f 1 f + -y h 3 f 2 f 2 + 0(h 4 ) 
 
 ct 
 k 2 = hf + a 2 g 2 f 1 f + -y h 3 f 2 f 2 + a x 3 22 h 3 f 2 f + 0(h 4 ) 
 
 Therefore, the second of eqs. (3.4) can be written as 
 
 (3.7) 
 
 1 
 
 1 
 
 1 
 
 
 n i 
 
 
 
 d 
 
 
 
 a l 
 
 a 2 
 
 r = 
 
 
 
 i 
 
 
 
 
 
 a 2 /2 
 
 a 2 / 2 
 
 
 
 
 
 
 l 
 
 
 
 
 
 a i e 22_ 
 
 
 
 
 
 
 l 
 
 where T is the matrix [y. ]. The rows of this equation correspond to 
 
 is 
 
 2 2 
 the terms f, f,f, f 9 f , and f,f> respectively. Clearly, the first column 
 
 T 
 of T is [1, 0, 0] and the first row is such that the sum of all rows is 
 
 [1, 0, 0], Thus, the first row and column can be calculated if values 
 
 can be found to satisfy the remaining six equations in eq. (3.7). 
 
 Therefore, we will drop the first row and column in (3.7) to get 
 
11 
 
 (3.8) 
 
 a 
 
 a. 
 
 a 2 
 
 
 2 
 a 2 
 
 r = 
 
 a i 3 22_ 
 
 
 where T is a 2 x 2 matrix. If this is to have a solution, a S ^ 0. 
 The equation in the (3, 1) position of eq. (3.8) implies Y 91 = 0. 
 Then the equation in the (2, 1) position implies Y, , = so that is 
 impossible to satisfy the equation in the (1, 1) position. Hence, a 
 3-stage third order method does not exist. However, we have demonstrated 
 the existence of a 4-stage method. 
 
 3.2. Fourth Order Methods 
 
 Equations (3.4) must be satisfied as identities in each of the 
 elementary differentials of f of orders up to p. This leads to a large 
 system of nonlinear equations for large p. The number of elementary 
 differentials for p < 5 are given in Table 3.1. 
 
 TABLE 3.1. Elementary Differentials of f 
 
 Order Number 
 
 Representation 
 
 f i f 
 
 2 2 
 
 V • £ i f 
 
 3 2 2 3 
 f 3 f , f 2 f x f , fl f 2 f , f x f 
 
 f 4 f4 ' f 3 f l f3 ' f l f 3 f3 ' f 2 f3 ' f2f l f2 > 
 2 2 2 2 4 4 
 
 f 2 (f 1 f) z , fjfjfjf , qf 2 f , fjf 
 
 The elementary differentials are represented in a prefix operator notation: 
 for example, f. is a ternary operator given by 
 
12 
 
 m m m 
 f~ abc = 
 
 3 f ij k 
 - a be 
 
 • i • i i t 3y.3y.3y, 
 
 i=l j=l k=l 'i J 2 k 
 where m is the dimension of the system and a, b, and c are the three 
 vector operands of f~. For an alternate notation due to Butcher, see 
 Stetter [15, p. Ill ff.]. There are 8 elementary differentials of order 
 up to four. Thus, the second of eqs. (3.4) can be written as 
 
 a r = 
 
 1 
 
 
 
 
 
 
 
 
 
 1 
 
 
 
 
 
 
 
 
 
 1 
 
 
 
 
 
 
 
 1 
 
 
 
 
 
 
 
 
 
 1 
 
 
 
 
 
 
 
 3 
 
 
 
 
 
 
 
 1 
 
 _o o o : 
 
 2 2 3 2 2 
 where the rows correspond to f, f-,f> f^f > f-,f> f-jf > fo^i^ > f-i^of » an( ^ 
 
 ■2~r 
 
 12 
 
 f,f, respectively. For a q-stage method A is an 8 by q matrix whose 
 entries are completely determined by the $.., and T is a q by 4 matrix. 
 As before, the first row and column of all matrices can be discarded 
 
 because the first row of A is [1, 1,..., 1] and the first column is 
 
 T 
 [1, 0,..., 0] . This leaves a system of 21 nonlinear equations in 3(q-l) 
 
 + q(q-l)/2 unknowns. Counting unknowns is of no direct value in 
 
 determining the existence of solutions, although it can be a guide to the 
 
 prospects. Although for q = 5 there are 21 equations in 22 unknowns, it 
 
 is shown in the appendix to the report [8] that no solution exists. 
 
 However, for q = 6, when there are 21 equations in 30 unknowns, there 
 
 exists a 9 parameter family of solutions. This is given in the appendix. 
 
 A particular case of this is 
 
13 
 
 B = 
 
 
 
 
 
 
 
 1 
 
 
 
 
 
 
 
 2 
 
 
 
 
 
 o - 
 
 
 
 
 
 
 
 
 
 
 
 
 
 2 
 
 
 
 /3 
 
 2_ 
 
 , r = 
 
 1 
 
 -5/6 
 
 4/9 
 
 -1/9 
 
 
 
 
 
 
 
 
 
 
 
 1/2 
 
 -4/9 
 
 1/9 
 
 
 
 7/3 
 
 -19/9 
 
 7/9 
 
 
 
 -3 
 
 10/3 
 
 -4/3 
 
 
 
 1 
 
 -11/9 
 
 5/9 
 
 3/4 9/4 
 
 1/2 1 1/2 
 1/12 2 1/4 2/3 
 
 Naturally, one wonders about higher-order techniques. For 
 example, since there are 17 elementary differentials of orders up to 
 5, we are led to a system of 16 by 4 equations (after discarding the 
 first row and column) in (q+8)(q-l)/2 unknowns. Nine is the first value 
 of q for which the number of unknowns exceeds the number of equations, 
 but it is not known if a solution exists for q = 9. Neither does it 
 seem practically important since the increase in the amount of work 
 does not seem to justify the small increase in order. 
 
14 
 
 4. IMPLICIT RUNGE-KUTTA METHODS AND STIFF EQUATIONS 
 
 It is well known that explicit methods do not function well for 
 stiff equations. However, normally immediately after a discontinuity 
 a differential equation is not stiff because the solution is in a transient 
 region where the step must be small to control stability. This can be 
 seen by considering the equation 
 
 (4.1) y 1 = A(y - F(t)) + g(t) where g = F' almost everywhere and X << 0. 
 Its solution is 
 
 (4.2) y = ce At + F(t) 
 
 After a small time, the terms ce ' is small compared to F, so the behavior 
 of F dominates the solution. However, if there is a discontinuity in the 
 right-hand side of (4.1) at some £, and this discontinuity is not consistent 
 with g = F' , then the solution (4.2) will have a discontinuity in c and F, 
 probably making the first term significant for another short time. If this 
 is the case, the explicit methods discussed earlier can be used. However, 
 in some cases the discontinuity may not stimulate the rapidly decaying 
 components sufficiently so that in principle a step size large compared 
 to 1/X can be used. If the explicit methods Of the previous sections are 
 applied to (4.1), we will obtain approximations to the scaled derivatives 
 of F plus polynomials in hX. For example, if the 6 stage method of order 
 4 is used, the approximation to the p-th derivative of y will be 
 
 h P y <P> = h P F <P> + c(h A)P + 0(h 5 ) + cefhX) 5 + ce JhX) 6 
 
 lb pb 
 
 where the 0(h) term depends on derivatives of F only. If c is small, 
 
 we would like the derivatives of F to dominate this expression. However, 
 
 if hX is large, the (hX) and (hX) terms may well dominate. For these 
 
 problems it is worth considering implicit RK methods. It is known 
 
15 
 
 (Butcher [2]) that q-stage implicit RK methods can achieve order 2q. 
 It is also known (Ehle [4]) that such methods are A-stable. In fact, 
 we are not concerned with A-stability or variants, such as stiff- 
 stability, here because only a single RK step is to be used. What we 
 would like is approximations to the first p derivatives of y to accuracy 
 h with the property that they are not large for hA in the negative 
 half plane. However, note that we would like these approximations to be 
 nonzero as hA ■*■ °° so that if there is, in fact, a significant ce 
 component present in the solution, it can be detected and the stepsize 
 used in the starter reduced. 
 
 It is known (Gear [7]) that a q-stage implicit RK scheme 
 applied to y' = Ay leads to the calculation of terms of the form 
 
 (4.3) R(hA)y= ^Ml y 
 
 where the degree of the polynomials P and Q cannot exceed q. Since we 
 
 wish to get an order p approximation to h y , we must be able to compute 
 
 a term of the form (hA) P + 0(hA) P . This implies that q _> p. In fact, 
 
 we will demonstrate that it is possible to get all the desired terms with 
 
 q = p, in which case we will be calculating the rationals 
 
 2 p 
 
 y + a y +.. .+ a . y F 
 
 R x (y) = l ^— 
 
 1 + a- y +. . .+ a y 
 1 P 
 
 2 3 p 
 
 y + ay +...+ a y K 
 
 R 2 (y) = 1 ^— 
 
 (4.4) 1 + a n y +. . .+ a y P 
 
 1 P 
 
 R (y) = 
 
 1 + a, y +.. .+ a y P 
 I P 
 
16 
 
 Note that for these approximations to be bounded when U is in the negative 
 
 half plane, and to be nonzero when y ■+ °°, all of the a must be nonzero 
 
 P 
 
 and all roots of Q(u) = 1 + a u +. . .+ a u must be in the positive half 
 
 plane. (It is probably sufficient that a ^ plus the root condition so 
 
 that at least the approximation to h y is nonzero when the stiff 
 
 components are significant. In fact, a should not be too large, so that 
 
 a large value of h y can be used as a warning to reduce the step and 
 
 try again.) 
 
 The existence of p-stage order p methods for starting can be 
 
 demonstrated without reference to elementary differentials — fortunately! 
 
 Let us assume that £. are the exact values of hf(y(t. + ot.h)). Consider 
 
 i l 
 
 P 
 k. = hf(y(t ) + Z $ k.) i = 1,..., p. 
 1 u j= i ^ J 
 
 If the a. are all different, the interpolation coefficients 3.. can be 
 chosen so that 
 
 y(t ) + I S k 
 
 -,=1 J -J 
 
 is an approximation to y(t^ + a.h) with an error of 0(h). Then, k. is 
 
 Ox i 
 
 ^ p+1 
 
 an approximation to k. with error of 0(h ). Therefore, solution of 
 
 P 
 (4.5) k. = hf(y(t ) + Z 3.. k.) 1=1,..., p. 
 
 will lead to approximations to k. with errors of 0(h ). With these 
 
 approximations to hy' at p different points we can obtain approximations 
 
 s (s) P+1 
 
 toh y, s=l,...,p with errors of 0(h ). (In fact, determination 
 
 of the interpolation coefficients 3.. leaves one degree of freedom in each 
 set, which, together with the free choice of the a.'s leaves 2q degrees 
 
17 
 
 of freedom. The $'s can be chosen so that y + Z 3.. k. has an 
 
 ij J 
 
 0(h ) error in its approximation to y(t n + ot.h). In that case, 
 
 s (s) 
 the errors in the estimates of the h y will be proportional to 
 
 terms in h y + 0(h ) and not involve separate elementary 
 
 differentials of order p + 1.) 
 
 It remains to show that the coefficients can be chosen so that 
 
 the root condition on Q(y) is satisfied, that is, so that Q(y) has zeros 
 
 in the positive half plane. It was shown in Gear [7] that Q(y) = 
 
 det(By - I) where B is the matrix of beta coefficients. If the betas 
 
 are completely determined by the alphas (as happens if the suggestion 
 
 from the last paragraph is followed), then Q(y) is completely determined 
 
 by the alphas.* If the alphas are chosen as Gaussian points corresponding 
 
 to the Butcher 2p-th order implicit RK, then it is known that Q(y) is of 
 
 degree p with all roots in the positive half plane (Ehle [4]). 
 
 * The zeros of Q are the inverses of the eigenvalues of 
 
 diag[a 1? a ,..., a ]A diag[l, 1/2, 1/3,..., l/p]A 
 12 p 
 
 -1 
 
 where A = 
 
 P" 1 
 1 a, a!f 
 
 P" 1 
 1 a 2 a£ 
 
 , p-1 
 
 1 a a 
 
 P P 
 
18 
 
 5. PROGRAMMING CONSIDERATIONS 
 
 Both the implicit and explicit techniques can be used either to 
 generate the Nordsieck vector, a set of function values, or derivative 
 values. However, it apparently lacks an error estimator, so what stepsize 
 should be used for the RK step, and what stepsize for subsequent multistep 
 methods? We propose to copy the technique followed in many of the better 
 codes available today, and use the estimate of the error in a method of 
 one order lower than that used! (This is often called local extrapolation: 
 it can be viewed as forming an error estimate, and using it not only for 
 step control but to correct the solution.) The RK Fehlberg methods [5] 
 are in this category if the higher order answer is accepted, as are 
 predictor-corrector methods if the corrector is one order higher than the 
 predictor and the predictor-corrector difference is used to estimate the error, 
 
 Since an estimate of h y is available in the p-th order method, 
 this can be used to select the stepsize for the first multistep step, and to 
 decide if the RK starter stepsize is reasonable. Assume, for the moment, 
 that we have a rough approximation to a step size. An RK starter step can 
 be taken and h y calculated. Suppose e is an error control parameter. 
 Three cases can occur. 
 
 (i) ||h P y (p) || « E 
 (ii) ||h P y (p) || » £ 
 
 (iii) Neither much too big nor much too small. 
 In the first case, the estimates of h y will be badly tarnished with 
 roundoff errors, so it is not possible to scale. Therefore, the RK starter 
 should be re-executed with a step about [e/ | |h P y P | |] larger. In the 
 second case, there is a danger of the estimates being badly tarnished by 
 truncation error, so the RK starter should be repeated with a smaller h. 
 
19 
 
 In the case of explicit methods, the ratio [e/ | |h y |] may be about 
 right; at least it should serve to bring h y within the large tolerance 
 allowed. In the case of implicit methods for stiff equations, an excessively 
 large step will cause |h y | to tend to l/|a |, so it may be difficult 
 to determine the amount of step reduction. However, if implicit methods 
 are used, some estimate of the Jacobian must be available, and this can be 
 used to get a crude upper bound for eigenvalues (e.g., the max norm), so 
 that a step that is certainly small enough to avoid stiffness problems 
 altogether can be tried the second time around if the first attempt fails. 
 In the third case, the initial stepsize for the multistep method can be 
 calculated from 
 
 h ° " M c|| h P/P»||' 
 
 p l I 
 
 where C is the error coefficient of the multistep method to be used. Then 
 P 
 
 the appropriate starting values can be calculated with Taylor's series. 
 
 What should be the criteria for "too large" and "too small" in the 
 
 acceptance tests. "Too small" must be based on roundoff considerations. As 
 
 s (s) 
 long as the roundoff errors in the estimates of h y cannot contribute an 
 
 error to the multistep starting values of size greater than £ (or something 
 
 a little smaller for safety), the RK starter can be accepted. 
 
 "Too big" is a harder question to answer. Assuming that 
 
 |y| = 1 (or that relative errors are used), my intuition is to require 
 
 that |h P y | not exceed e P P . This is based on the idea that if 
 
 (s) - iS , , . , , p+1 (p+1) ~ , , 
 
 y = A , then this bound causes h y = £, so that the errors in the 
 
 RK starter are of about the right size. 
 
20 
 
 6. DETECTING DISCONTINUTIIES 
 
 If an automatic step control integrator is given no information 
 about discontinuities, the step control mechanism can be very inefficient 
 in the neighborhood of a discontinuity. This occurs because an attempted 
 step that straddles a discontinuity yields a very large error estimate, 
 usually causing the code to reject the step and reduce the stepsize 
 sharply so that the step not only no longer straddles the discontinuity 
 but often so that it falls far short of the discontinuity. Consequently, 
 several small steps are taken until the step is once again increased to 
 a reasonable size so that an attempted step again straddles the 
 discontinuity. This process is repeated with a large cost in function 
 evaluations. Codes have used a number of techniques to reduce the problem. 
 DIFSUB [6], an early code, reduces the order to one and restarts when this 
 type of difficulty is first encountered, and then quits if it occurs three 
 times in succession — a simple but not too useful solution. Some codes 
 resort to binary division of the interval when the error tests are 
 inconsistent. (The error has a known asymptotic behavior if no 
 discontinuities are present; it is possible to check for gross deviation 
 from such behavior as an indication of discontinuities.) Binary division 
 may be as good as any strategy if no further information about the 
 discontinuity is available. However, if we know exactly when each 
 discontinuity occurs, it is far more efficient to integrate exactly to the 
 discontinuity and then restart. 
 
 Discontinuities can arise from two sources; the driving 
 (t-dependent) terms and the dependent variables. An example of a driving 
 term discontinuity is smog simulation. The sun's energy input is a time- 
 dependent term that has two discontinuities each 24 hours as the sun rises 
 and sets. Dependent variable discontinuities frequently arise in simulation. 
 
21 
 
 For example, a relay leads to equations of the form 
 
 if |lR| > I critical then switch **■ on 
 
 if |lR| < I min then switch ■*■ off 
 
 if switch = on then V = 
 else 1=0 
 where IR is the current in the relay coil, and V and I are the voltages across 
 and the current through the relay contracts. This means that as IR (which 
 is a function of the dependent variables) increases to the point it passes 
 I critical, the relay turn on current, the coefficients in the differential 
 equations change discontinuously. A similar problem arises in some 
 mathematical models that, although infinitely dif f erentiable , exhibit such 
 sudden changes they are indistinguishable from discontinuities. For 
 example, a diode is sometimes modeled by an exponential function such as 
 
 t 40V IN 
 
 I = c(e - 1) 
 where c is the reverse leakage current. The constant 40, or some large 
 number, causes a behavior similar to that shown in Figure 1. 
 
 -*■ V 
 
 Figure 1. Analytic function which behaves like discontinuity. 
 In such cases either the model should be changed to a discontinuous model 
 (e.g., a piecewise linear function) or detection schemes which indicate 
 when the nature of the model changes radically should be added. 
 
 The basic problem is to detect when a discontinuity occurs in 
 time so that when a step is about to straddle it, the stepsize can be 
 
22 
 
 reduced. This is trivial for time-dependent discontinuities, so let us 
 
 consider the problem of functions of dependent variables. Let us suppose 
 
 that a discontinuity occurs when an expression e(y) crosses zero. Since 
 
 y is a function t, we can evaluate e(y(t)) at each integration step. We 
 
 could use a set of past values e ., i = 0,..., k from a multistep method 
 
 n-i 
 
 to estimate the time at which the next zero crossing will occur and 
 reduce the step of that time is within the next interval. However, that 
 is time-consuming because of an inverse interpolation in each interval 
 and error prone because extrapolation gives larger errors than interpolation. 
 Therefore, we would like to be able to take the next step and then interpolate, 
 This is possible only if there is no discontinuity in the step. The 
 technique employed by many people (see Carver [3]) is to delay changing 
 the "switch" until the completion of a step. Thus, if e(y) is an 
 expression such that a switch is to change when e(y) changes from negative 
 to positive, e(y) is evaluated only at the end of each step. If e(y) 
 remains negative, the integration proceeds normally. However, if the value 
 of e(y) is found to be positive at the end of a step, it indicates a 
 discontinuity somewhere in that step. Inverse interpolation can then be 
 employed to determine the t value at which e(y) is zero. Interpolation 
 can then be employed to calculate the dependent variables, although it is 
 probably better to integrate from the last value of t, and possibly iterate 
 this procedure one time. 
 
23 
 
 Acknowledgment 
 
 I would particularly like to acknowledge some discussions with 
 F. Cellier of ETH, Zurich, on this problem. 
 
 Bibliography 
 
 [1] Bettis, D. G., "Efficient Embedded Runge-Kutta Methods," in Numerical 
 Treatment of Differential Equations , Lecture Notes in Mathematics #631, 
 eds. Burlirsch, R., Grigorieff, R. D. t and Schroder, J. Springer-Verlag, 
 Berlin, 1976. 
 
 12] Butcher J. C., "Implicit Runge-Kutta Processes," Math . Comp . 18 , pp 
 50-64, 1964. 
 
 [3] Carver M. B., "Efficient Handling of Discontinuities and Time Delays 
 in Ordinary Differential Equations," Proceedings of the Simulation '77 
 Symposium, Montreux, Switzerland, ed. Hamza, M., Acta Press, Zurich. 
 
 14] Ehle B. L., "A-Stable Methods and Pade Approximations to the 
 Exponential," SIAM Journ . Math . Anal . 4^, pp 671-680, 1973. 
 
 [5] Fehlberg E., "Klassiche Runge-Kutta Formeln vierter und niedrigerer 
 Ordnung mit Schrittweiten-Kontrolle und inre Anwendung auf 
 Warmeleitungsprobleme," Computing 6_, pp 61-71, 1970. 
 
 16] Gear C. W. , "DIFSUB for the Solution of Ordinary Differential 
 Equations," CACM 14, pp 185-190, 1971. 
 
 [7] Gear C. W. , "Rational Approximations by Implicit Runge-Kutta Schemes," 
 BIT 10, pp 20-22, 1970. 
 
 18] Gear C. W. , Numerical Initial Value Problems for Ordinary Differential 
 Equations , Prentice-Hall, Englewood Cliffs, New Jersey, 1971. 
 
 (9] Hindmarsh A. C, "GEAR: Ordinary Differential Equation Solver," Rept. 
 UCID-30001, rev. 3, Lawrence Livermore Laboratory, CA, 1974. 
 
 [10] Hull T. E., Enright W. 11. , Fellen B. M. , and Sedgwick A. E. , 
 "Comparing Numerical Methods for Ordinary Differential Equations," SIAM 
 Journ . Num . Anal . £, pp 603-637, 1972. 
 
 Ill] Krogh F. T., "VODQ/SVDQ/DVDQ - Variable Order Integrators for the 
 Numerical Solution of Ordinary Differential Equations," Section 314 
 Subroutine Write-up, Jet Propulsion Laboratory, Passadena, CA, 1969* 
 
 112] Nordsieck A., "On the Numerical Integration of Ordinary Differential 
 Equations" Math . Comp . 16, pp 22-49, 1962. 
 
24 
 
 [13] Shampine L. P., and Gordon M. K. , Computer Solution of Ordinary 
 Differential Equations : the Initial Value Problem , W. H. Freeman, San 
 Francisco, 1975. 
 
 [14] Stetter, H* J., "Aymptotic Expansions for the Error of Discretization 
 Algorithms for Nonlinear Functional Equations," Num . Math . 7_, pp 18-31, 
 1965. 
 
 [15] Stetter, H. J., Analysis of Discretization Methods for Ordinary 
 Differential Equations , Springer Tracts in Natural Philosophy #23, 
 Springer-Verlag, New York, 1973. 
 
25 
 
 APPENDIX 
 
 A. Non-Existence of 5-Stage Method of Fourth Order 
 
 For a 5-stage, fourth order method to be possible we must 
 be able to solve 
 
 (A.l) 
 
 where 
 
 A = 
 
 1 
 
 2 
 a i 
 
 
 
 a 
 
 22 a i 
 
 a. 
 
 Ar = 
 
 1 
 
 
 
 
 
 
 
 2 
 
 
 
 
 
 1 
 
 
 
 
 
 
 
 6 
 
 
 
 
 
 3 
 
 
 
 
 
 2 
 
 
 
 
 
 1 
 
 a. 
 
 (B 32 a x + B 33 a ) 
 
 a. 
 
 a 
 
 a, 
 
 (3 42 a l + 3 43 a 2 + 3 44 a 3 } 
 
 a, 
 
 a A A 34 
 
 \l a l + B « 4 + B 44 a 3 
 
 B « A 32 + \k A 33 
 and T is a 4 x 3 matrix. 
 
 Clearly, the columns of V are independent because the columns of 
 the right-han 1 side of eq. (A.l) are independent. Let B be the 4x4 matrix 
 
 
 
 a 2 A 32 
 
 a 3 A 33 
 
 
 
 a i A 32 
 
 B 32 a i + 3 33 a 2 
 
 
 
 
 
 3 33 A 32 
 
26 
 
 consisting of the last four rows of A. From (A.l) the first two columns 
 
 of T are independent eigenvectors of B for the eigenvalue zero. Therefore, 
 
 B has rank no greater than 2. If a = 0, the method degenerates to a 
 
 4-stage method, so we can assume a ^ 0. Hence the first row of B is 
 
 independent of the other three because of the zeros in the remaining 
 
 places in column 1. Hence the last three rows must have rank at most 1. 
 
 Since A = 0, this is possible only if A = A = 0. Since A = a A 
 
 72 52 62 62 1 32 
 
 and a ^ 0, A = 0, so that 3 ?? = 0. Hence A = 0. Again, for the 
 
 last three rows to have rank 1, A_„ = A.„ = 0. Note that A_. = 
 
 53 63 74 
 
 3,~ A„ + 3,, A„_ = 3,/ A„„ because we have shown that A„„ = 0. Clearly, 
 43 32 44 33 44 33 32 
 
 A . f (otherwise last row of A is zero and it is impossible to satisfy 
 eq. (A.l)). Hence A ~ ^ 0. Therefore, since = A ~ = a A , a~ = 0. 
 The A matrix is now 
 
 A = 
 
 a 
 
 ex. 
 
 
 
 
 
 a 
 
 
 
 33 
 
 34 
 
 1 
 
 3 
 
 a 2 
 
 
 
 3 
 
 a 4 
 
 
 
 
 
 
 
 A 54 
 
 
 
 
 
 
 
 A 64 
 
 
 
 
 
 
 
 K, 
 
 74 
 
 where A_., k, . and A_, are nonzero (so that eq. (A.l) can be solved). 
 54' 64 74 
 
 Therefore V = T =0. Consideration of the first, second and fourth 
 ' 41 42 
 
 rows of (A.l) as equations for T and T i = 1, 2, 3 leads to a 
 singular system without a solution, so a 5-stage, fourth order method 
 does not exist. 
 
27 
 
 B. A 9-Parameter, 6-Stage, Fourth Order Method 
 Again, we must solve 
 
 (B.l) 
 
 Ar = 
 
 1 
 
 
 
 
 
 
 
 2 
 
 
 
 
 
 1 
 
 
 
 
 
 
 
 6 
 
 
 
 
 
 3 
 
 
 
 
 
 2 
 
 
 
 
 
 1 
 
 where 
 
 A = 
 
 a 
 
 a. 
 
 a 
 
 a l 
 
 
 
 
 
 
 a. 
 
 a 
 
 a. 
 
 a P 
 2 2 
 
 a i P 2 
 
 a. 
 
 a 
 
 a. 
 
 °3 P 3 
 
 a, 
 
 a. 
 
 a 4 P A 
 
 a r 
 
 a r 
 
 a 5 P 5 
 
 with 
 
28 
 
 r 
 
 p„ = 
 
 6 22 a i 
 
 (B.2) 
 
 P~ = 
 
 P 4 
 
 P. = 
 
 e 32 a ± + 3 33 a 2 
 
 8 42 a x + 6 43 a 2 + 3 44 a 3 
 
 3 52 a x + 3 53 a 2 + 3 54 a 3 + 8 55 c* 4 
 
 Qo = 
 
 B 22 c^ = a ± P 2 
 
 (B.3) 
 
 Q 3 = 3 32 a* + 8 33 a\ 
 
 Q 4 = 3 42 a i + 3 43 a 2 + 3 44 a 3 
 
 1^5 
 
 8 52 a 1 + 3 53 a 2 + 3 54 c* 3 + 3 55 « 4 
 
 R 2 = 
 
 (B.4) 
 
 R 3 - 8 33 P 2 
 
 R 4 = e 43 P 2 + 3 44 P 3 
 
 R 5 = S3 P 2 + 3 54 P 3 + 3 55 P 4 
 
 The first six parameters we will pick are a., i = 1,..., 5, and P~. For 
 the moment consider P_ as another parameter. These should be chosen so 
 that the first 2x2 and 3x3 minors of A are nonsingular. (Some 
 additional nonsingularity conditions will be imposed later.) Note that 
 these determine 3 ?9 trivially from eq. (B.2.1). We first want to show 
 that these determine $ ~ and & ~ anc * hence Q and R_ from (B.3. 2) and 
 (B.4. 2). We will also show that they allow R. to be determined in terms 
 of P. and Q., i = 4, 5. As before, let B be the last four rows of A. It 
 has rank no greater than 3 because it is a 4 x 5 matrix with two vectors, 
 
 F and T , such that 
 
 Br. = o 
 
29 
 
 T 
 Therefore, there exists a row vector c = (c, , c , c„ , c.) that 
 
 — 12 3 4 
 
 (B.5) c T B = 
 
 From the first column of B we have immediately that c = 0. From the 
 second column we find that a c + a c = 0. Without loss of generality, 
 take c 9 = a 1 , c„ = -a_. From the third column we obtain 
 
 (B.6) a 3 P 3 c 2 + Q 3 c 3 + R 3 c 4 = 
 
 T 
 If r is the third column of T, we have from (B.l) that Bl\ = [6, 3, 2, 1] . 
 
 Post multiplying (B.5) by this column vector we have 
 
 c T Br 3 = c T [6, 3, 2, 1] T = 
 
 or 
 
 3c_ + 2c» + c. = . 
 
 2 3 4 
 
 Thus, c,= -2c - 3c 9 = 2a - 3a . Inserting this in (B.6) we get 
 
 a 2 Q 3 -(2a 2 - Sa^ = a ± a 3 P 3 
 
 If we substitute (B.3.2) and (B.4.2) in this we get 
 
 (B.7) a^ a 2 3 32 + (a 2 - (2a 2 - Sa^P^e^ = ^ a 3 P 3 
 
 This, together with (B. 2. 2), forms a pair of equations for 3^ 9 and 3„«. 
 
 Consideration of the fourth and fifth elements of eq. (B.5) yields 
 
 the relation 
 
 a. P. c„ + Q. c + R. c. =0 
 ix2 i 3 l 4 
 
 Using the values found for c . we get 
 
 3 
 
 ou Q. - ol a, P, 
 
 (B.8) R. =-^-i k-^- 1 
 
 x 2a - 3a 
 
 Note that the last row of A is now a linear combination of the fifth and 
 sixth rows consistent with the right-hand-side of (B.l). Hence we can 
 now ignore the last row of A. 
 
30 
 
 We now break A into a 3 x 5 matrix C of its first three rows, 
 and the fourth through sixth rows, D. Then, eq. (B.l) can be written, 
 in part, as 
 
 10 
 
 (b.9) cr = 2 
 
 _0 1 0_ 
 
 We have assumed that the first 3x3 principal minor of A is nonsingular 
 so the solution of (B.9) can be expressed as 
 
 (B.10) 
 
 r = 
 
 A ll 
 
 A 22 
 
 
 
 A 22 
 
 A 22 
 
 
 
 A 33 
 
 A 32 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 y il 
 
 y i2 _ 
 
 y 21 
 
 y 22 
 
 y 31 
 
 u 32 
 
 \l 
 
 
 
 
 
 y 52 _ 
 
 where jj. = {\i. .}, i = 1, 2 are null right vectors of C and E is any 2x3 
 matrix. Note that since the first three column of C are 
 
 a i a 2 a 2 
 
 a, a 2 a 2 
 
 o p 2 p 3 
 
 the values of A., are rational polynomials in P„, linear in the numerator 
 and denominator. 
 
 The remaining steps consist of ensuring that we can choose E in 
 (B.10) so that 
 
31 
 
 (B.ll) 
 
 Dr = 
 
 
 
 
 
 6~ 
 
 
 
 
 
 3 
 
 
 
 
 
 2_ 
 
 thus solving (B.l). 
 
 Let A. = {A. .} i = 1, 2 be column vectors and let e.. be the 
 -i l 13 
 
 elements of E. From (B.ll) we observe that 
 
 D Al + e H °% + e 21 % = ° 
 
 D -2 + 6 12 D % + e 22 ^2 = ° 
 Hence, D\ and DA must be in the space spanned by (Dy , Dy ). Since the 
 third column of (B.ll) tells us that 
 
 e 13 PUj + e 23 D M2 = [6, 3, 2] T 
 
 we know that [6, 3, 2] is in this same space. If DA and DA are independent, 
 
 T 
 which we now assume, then they must also span this space. Hence [6, 3, 2] 
 
 is in span [DA , DA ] , or 
 
 (B.12) 
 
 det 
 
 DA 1 , DA 2 , 3 
 
 = 
 
 This condition is used to determine P_ in terms of the six parameters a. 
 
 3 i 
 
 and P ? . (In fact, a, and a^ do not appear in the value.) 
 
 If we examine the structure of jj and y_ , we see that they are 
 linear in P, and P,., respectively, but otherwise depend only on the parameters 
 ot. and P ? . Furthermore, y . , and y depend only on the parameters, and not 
 on P, and P-. Because of this Dy is linearly dependent on P, and Q, , while 
 Dy is linearly dependent on P, and Q . From the earlier argument that DA 
 and D_A„ are in the space spanned by Djj and Dy_ 9 we see that Dy and Dy must 
 be in the space spanned by DA. and T)X . Hence 
 
32 
 
 (B.13) det [DA , DA 2 , Dy. ] =0 
 
 and 
 
 (B.14) det [D,^, DA 2 , D]J 2 ] = 
 
 Since DA . have been determined in terms of the parameters a. and P„, 
 i l 2 
 
 (B*9+i) i = 4, 5, is a linear relation between P. and Q.. Thus, we can 
 select P, and P- as parameters, determining Q, and Q^, and hence R, and 
 R,. from (B.8). (In some cases, vanishing coefficients fix P, and/or P_, 
 leaving Q, and/or Q- as free parameters.) Now we observe that eqs. (B.i.3) 
 i = 2, 3, 4 are a set of three nonsingular equations which determine 3,. 
 j = 2, 3, 4 uniquely. Finally, note that eqs. (B.i.4) i = 2, 3, 4 are a set 
 of three equations for $t-. j = 2, 3, 4, 5. 3c,- can be chosen as a parameter, 
 determining the remaining betas. Now the matrix A is known, so the 
 determination of T from (B.l) is straightforward. 
 
 Each of the steps given above is linear with one exception- 
 solution of eq. (B.12) for P~. We will now show that this will always 
 have a real solution, and thus demonstrate the existence of solutions 
 almost everywhere in terms of the parameters a., P , P , P , and 3,-c- 
 (The existence of one solution indicates that the linear equations are 
 nonsingular almost everywhere. Such a solution is given in the paper.) 
 
 There seems to be no simple approach to demonstrating that (B.12) 
 in fact reduces to a quadratic that can be factored over the reals except 
 to plod through part of its derivation. We will do this for the benefit 
 of the determined. 
 
 If (B.7) and (B.2.2) are solved for 3 32 and 3^ 3 , and the result 
 
 is used to calculate Q-. from (B.3.2) we get 
 
33 
 
 Q 3 = ot 1 
 
 a 2 a 3 (a i~ a 2 ) + (2a 2 _3CX l )P 2 
 
 a 2 (a i" a 2 ) + ( 2a 2 " 3a i )P 2 
 
 Computing the X's in (B.10) we find that 
 
 h~ [ 
 
 2 2 2 
 a 3 P 2 - a 2 P 3 a i P 3 
 
 A ' A 
 
 2 „ 
 
 ' a i P 2 T 
 -^ , 0, 0] 1 
 
 (B.15) 
 
 A 2 = 
 
 a 2 a 3 (a 2 -a 3 ) 
 + 2a 2 P 3 -2a 3 P 2 
 
 a i a 3^3 _a i^ a l a 2^ a l~ a 2^ 
 
 " 2a i P 3 
 
 + 2ai P 2 
 
 , 0, 
 
 where A = P_ a a (a -a ) - P~ a a (a -a ) . From now on we will work with 
 A_A and AA ? since this scaling has no effect on eq. (B.12) as long as 
 A ± 0. Note that from (B.15) AA and AA 2 are dependent (parallel) if 
 
 P„ a„(a -a ) = P~ a„(a_-a ). (This can be verified by routine manipulation.) 
 
 T 
 Hence, det [ADA , ADA , [6, 3, 2] ] = has a root satisfying this 
 
 condition — which happens to be A = so should be ignored. By computation 
 
 we find that 
 
 ADA = 
 
 a^(a r a 3 )P 2 - ^(a^^ 
 
 a 1 P 2 P 3 (a 2 -a 3 ) 
 
 a i a 2 ^ a i" a 2^ ( a 2 ~ a 3^ P 2 P 3 
 a 2 (a i~ a 2 ) + ( 2a 2 ~ 3a i )P 2 
 
 and 
 
34 
 
 ADA 2 = 
 
 2 2 2 
 
 a^a^o^Cc^-a ) + 0^(0^-0^) + a (a -a ) 
 
 + 2a 1 [a 2 (a 1 -a 2 )P 3 - a^a^a^P^ 
 
 a 1 a 2 a 3 [P 3 (a 1 -a 2 ) - P^-a^ 
 
 + 2a x P 2 P 3 (a 3 -a 2 ) 
 
 2 2 a 2 a 3 (a -a ) + (2a -3a )P. 
 a a (a -a )P + a a (a -a )P — = 
 
 1 J J L Z l Z l l J a 2 ( V a 2 ) + < 2 V 3a l )P 2 
 
 2a 1 a 2 (a 1 -a 2 )(a 2 -a 3 )P 2 P 3 
 
 2 
 a (a -a ) + (2a 2 -3a ;[ )P 2 
 
 Note that these are both linear in P~. Hence det [ADA , ADA. , [6, 3, 2] ] 
 is quadratic in P~. We have already seen that P„a~(a -a ) - P„a ? (a 9 ~a ) = 
 is a zero of this that is of no interest because A = 0. Hence there is 
 exactly one other real zero for P~. 
 
Form AEC-427 
 
 (6/68) 
 
 AECM 3201 
 
 U.S. ATOMIC ENERGY COMMISSION 
 
 UNIVERSITY-TYPE CONTRACTOR'S RECOMMENDATION FOR 
 
 DISPOSITION OF SCIENTIFC AND TECHNICAL DOCUMENT 
 
 I See Instructions on Reverse Side ) 
 
 1. AEC REPORT NO. 
 
 COO-2383-0053 
 
 2. TITLE 
 
 RUNGE-KUTTA STARTERS FOR MULTISTEP METHODS 
 
 3. TYPE OF DOCUMENT (Check one): 
 
 [x| a. Scientific and technical report 
 
 I I 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): 
 
 PCI a. AEC's normal announcement and distribution procedures may be followed. 
 1] b. Make available only within AEC and to AEC contractors and other U.S. Government agencies and their contractors. 
 ~2 c. Make no announcement or distribution. 
 
 5. REASON FOR RECOMMENDED RESTRICTIONS: 
 
 6. SUBMITTED BY: NAME AND POSITION (Please print or type) 
 
 C. W. Gear 
 
 Professor and Principal Investigator 
 
 Organization 
 
 Department of Computer Science 
 University of Illinois 
 Urbana, Illinois 61801 
 
 Signature 
 
 ^Kjji^tj^* 
 
 Date 
 
 September 1978 
 
 FOR AEC USE ONLY 
 
 7. AEC CONTRACT ADMINISTRATOR'S COMMENTS, IF ANY, ON ABOVE ANNOUNCEMENT AND DISTRIBUTION 
 RECOMMENDATION: 
 
 8. PATENT CLEARANCE: 
 
 CD a. AEC patent clearance has been granted by responsible AEC patent group. 
 I) b. Report has been sent to responsible AEC patent group for clearance. 
 I I c. Patent clearance not required. 
 

 BIBLIOGRAPHIC DATA 
 SHEET 
 
 1. Report No. 
 
 UIUCDCS-R-78-938 
 
 4. Title and Subtitle 
 
 RUNGE-KUTTA STARTERS FOR MULTISTEP METHODS 
 
 3. Recipient's Accession No. 
 
 5. Report Date 
 
 September 1978 
 
 6. 
 
 7. Author(s) 
 
 C. W. Gear 
 
 8. Performing Organization Rept. 
 
 No -UIUCDCS-R-78-938 
 
 9. Performing Organization Name and Address 
 
 Department of Computer Science 
 University of Illinois 
 Urbana, Illinois 61801 
 
 10. Project/Task/Work Unit No. 
 
 11. Contract /Grant No. 
 
 US ENERGY/EY-76-S-02-23*3 
 
 12. Sponsoring Organization Name and Address 
 
 U. S. Department of Energy 
 Chicago Operations Office 
 9800 South Cass Avenue 
 Argnnnp, Tllinnia fiOA^Q 
 
 13. Type of Report & Period 
 Covered 
 
 14. 
 
 15. Supplementary Notes 
 
 16. Abstracts 
 
 Runge-Kutta-like formulas are presented which enable a multistep method to start 
 or restart at a high order after just one RK step. These greatly improve the 
 efficiency of multistep methods in one situation in which they were previously 
 outperformed by RK methods — problems with frequent discontinuities or sudden 
 large increases in derivatives, either of which cause an automatic program to 
 reduce order and stepsize suddenly. 
 
 17. Key Words and Document Analysis. 17a. Descriptors 
 
 Runge-Kutta 
 
 Multistep 
 
 Integration 
 
 Ordinary Differential Equations 
 
 7b. Identifiers/Open-Ended Terms 
 7c. COSATI Field/Group 
 
 
 
 8. Availability Statement 
 
 unlimited 
 
 19. Security Class (This 
 Report) 
 
 UNCLASSIFIED 
 
 21. No. of Pages 
 
 36 
 
 20. Security Class (This 
 Page 
 
 UNCLASSIFIED 
 
 22. Price 
 
 °"M NTIS-35 ( 10-70) 
 
 USCOMM-OC 40329-P71 
 
rug ia».