■■nfi um ' 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, and underlining of books 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 JUN 1 8 1936 SEP 23i99 7 JUN 3 1997 When renewing by phone, write new due date below previous due date. L162 Digitized by the Internet Archive in 2013 http://archive.org/details/numericalintegra894fatu 7 £dg ' U ^U&CDCS-R-77-894 COO-2383-0044 UILU-ENG 77 1751 ^7 NUMERICAL INTEGRATORS FOR STIFF AND HIGHLY OSCILLATORY DIFFERENTIAL EQUATIONS by Simeon Ola Fatunla October 1977 UIUCDCS-R-77-894 NUMERICAL INTEGRATORS FOR STIFF AND HIGHLY OSCILLATORY DIFFERENTIAL EQUATIONS by Simeon Ola Fatunla* October 1977 DEPARTMENT OF COMPUTER SCIENCE UNIVERSITY OF ILLINOIS AT URB ANA- CHAMPAIGN URBANA, ILLINOIS 61801 *Author on leave from the Department of Mathematics, University of Benin, Benin-City, Nigeria. Work supported in part by the U. S. Energy Research and Development Administration under contract US ERDA/EY-76-S-02-2383. ACKNOWLEDGMENT The author is indebted to Professor C. W. Gear for his invaluable suggestions, criticism and advice. He has gone out of his way to ensure that I was comfortable during my short visit at the Department of Computer Science. -1- ABSTRACT Some L-stable fourth order explicit one step numerical integration formulas which require no matrix inversion are proposed to cope effectively with systems of ordinary differential equations with large Lipschitz constants (including those having highly oscillatory solutions). The implicit integration procedure proposed in Fatunla [10] is further developed to handle a larger class of stiff systems as well as those with highly oscillatory solutions. The same pair of nonlinear equations as in [10] is solved for the stiffness/oscillatory parameters. However, the nonlinear systems are transformed into linear forms and an efficient computational procedure is developed to obtain these parameters. The new schemes compare favorably with the backward differentiation formula (DIFSUB) of Gear [12, 13] and the blended linear multistep methods of Skeel and Kong [22], and the symmetric multistep methods of Lambert and Watson [16]. -2- 1. INTRODUCTION The development of numerical integration formulas for stiff as well as highly oscillatory systems of differential equations has attracted considerable attention in the past decade. The reason for this cannot be farfetched, realizing that the mathematical models of physical situations in kinetic chemical reactions, process control and electrical circuit theory often generate systems of ordinary differential equations whose Jacobians have at least one eigenvalue with a very negative real part or very large imaginary part. Both situations are respectively described as stiff and highly oscillatory. Consider the following model problems (1.1) y' = A(y-c(x)) + c'(x), y(a) = y Q , where y(x) £ R , X a complex constant with re X << and c(x) has the property |c'(x)| ~ in the finite interval a < x < b; (1.2) !- e J, z' = z» z(a) = z Q , ' -co -e 2 with z(x) E R and co >> 0, e a positive constant such that e ~ 0. Problem (1.1) has theoretical solution (1.3) y(x) = c(x) + y Q e Ax whose component c(x) is slowly varying in the specified interval while the second component decays rapidly in the transient phase. The analytic solution to problem (1.2) is given by (1.4) / a cos cox + 6 sin cox'. z(x) = e CX 1+8 cos cox - a sin cox, where a and 6 are the arbitrary constants of integration. -3- The transitory phase for problem (1.1) is of the order of -1/A while that of problem (1.2) is the entire interval of integration with o)/2tt complete oscillations per unit interval. Almost invariably, most conventional numerical integration solvers cannot effectively cope with problems (1.1) and (1.2) as they lack adequate stability characteristics. Any attempt to impose the stability properties will in effect constrain the integration mesh-size to be intolerably small. This may ultimately have adverse effects on the accuracy due to propagation of roundoff errors. Besides, the computing time and cost may be too excessive. Existing algorithms developed for problems of the type (1.1) can be classified into the following categories: (i) generalized Runge-Kutta schemes (ii) implicit Runge-Kutta schemes (iii) trapezoidal rule with extrapolation (iv) multiderivative multistep formulas (v) backward differentiation multistep formulas. As of now the most widely used numerical integration code for stiff systems is the class (v) schemes particularly Gear's DIFSUB [12, 13]. The method is efficient and reliable provided the eigenvalues of the Jacobian are not close to the imaginary axis where the higher order schemes exhibit poor stability properties (as evidenced from example 2). Dahlquist [5] established that neither an explicit linear multistep scheme of any order nor an implicit multistep method of order greater than two can be A-stable. He proved that the trapezoidal rule (which of course is not 1 3 (3) L-stable as y /y , -*■ 1 as Ah -* -°°) has the smallest error of (±)to h y. ,. J t t-1 1^ V.x; Other schemes which behave better than DIFSUB when the eigenvalues are close -4- to the imaginary axis includes the second derivative multistep method Enright [6]. As regards problems of the type (1.2), earlier efforts include Gautschi'tl [11] nonlinear multistep schemes which produce exact solution to algebraic (or trigonometric) polynomials up to certain degrees. The main drawback of this scheme is that it requires an a priori knowledge of the period of the systems under consideration. Other numerical integration solvers to oscillatory systems include Amdursky and Ziv [1, 2, 3], Snider and Flemming [23], Miranker, et al. [18], Miranker and Wahba [19], Miranker and Veldhuizen [21], and Fatunla [8, 9]. Unfortunately, none of these existing routines has been properly (adequately) put to test for various kinds of oscillatory problems. In this paper, we propose some two-point numerical integration formulas which effectively cope with systems of ODEs whose characteristics are identical to problems (1.1) and (1.2). We shall consider initial value problems (1.5) y' = f(x,y), y(0) = y Q with y(x) £ R in the finite interval S = [0,x f ] C R where x = Nh for some positive integer N > 0. It is assumed that y(x) is sufficiently dif ferentiable. We adopt the vector notation: y = ( y' y, . . . , y) , f=(f, f,..., f). The numerical estimates y to the theoretical solution y(x ) at the points x = th, t = 0(1)N are to be generated. On every subinterval [x , x +h], the theoretical solution y(x) is approximated by either the interpolating function fix -fi„x (1.6) F(x) = (I-e )A - (I-e l )B + C, A, B, nnd C being m-tuples with real entries, I is the identity matrix whilst r u and U are diagonal (stiffness/oscillatory) matrices or -5- a x n*x (1.7) F(x) = (I-e )R + (I-e )R* + S where R, S are m-tuples with complex entries and (*) denotes complex conjugate. The choice of interpolation formula is determined by equation (3.12). The following definitions are worthwhile: Definition 1 . A one-step numerical integration scheme is considered L-stable if apart from being A-stable, when it is applied to the scalar initial value problem (1.8) y' = Ay, y(0) = n (X being a complex constant with negative real part) , the resultant numerical solution is given by (1.9) y fc+1 = u(Xh)y t with the characteristic equation y(Xh) having the property: (1.10) lim |y(Xh)| = 0. re(Xh)-x» Definition 2. A numerical integration scheme is said to be exponentially fitted at a complex value X = X if when it is applied to the initial value problem (1.8) with exact initial condition, the characteristic equation y (Xh) satisfies the relation X h (1.11) y(X Q h) = e U . Liniger and Willoughby [18] and Jackson and Kenue [14] have respectively discussed A-stable one- and two-step numerical integration methods which are exponentially fitted at infinity. Both schemes ensure exponential fitting by a suitable choice of a free parameter. This approach was further extended to construct a stiffly stable k-step method of order k+2 in Enright [6]. -6- We shall construct from both equations (1.6) and (1.7) explicit one-step numerical integration formulas of fixed order four which possess adequate stability and convergent characteristics to cope with both stiff and highly oscillatory systems of ordinary differential equations. For linear systems the interpolating functions are global in the sense that the stiffness/oscillatory matrices have constant entries which are determined by solving a set of linear equations at the first step of the integration procedure. The implicit scheme proposed in Fatunla [10] is further developed using the interpolating function (1.7). The need to solve nonlinear equations for the stiffness/oscillatory parameters is eliminated. -7- 2. DEVELOPMENT OF THE INTEGRATION FORMULAS Let y . . denote the numerical estimate of the theoretical solution t+j y(x ,.) and x = x , , and adopt the notation f ,. ■ f(x_.., y_. .) and set 7 t+j t+j r t+j t+j J t+3 (2.1) y t+ . ■ F( V] ), 3=0,1 and (2.2) f^ k) = i™ y k = 0, 1. The imposition of. these constraints on the interpolating function (1.6) gives the integration formula (2.3) y t+1 - y t + « t + sf <» , where (2.4) r = Oy-^a) , and (2.5) s = y +o ; Y and a are diagonal matrices with entries x ft h (2.6) i Y =| r^ , i = l(l)m, 1 J2 1 ( 1 fi 1 + 1 fl 2 ) and (2.7) ^ = -^ : t-^- , i = l(l)m. 1 fi 2 ( 1 fi 1 + 1 fi 2 ) In the event that a component of the stiffness/oscillatory matrix (ft say) does vanish, by L'hopital rule, the corresponding component of Y in (2.4) is (2.8) 3 -h_ J ft "2 and the resultant integration formula is given by -8- (2.9) j y = j y + h j f + (-7^ + j a) j f^ 1) . We now consider the case of the complex interpolating function (1.7). Let the parameters be expressed as (2.10) tt = A + iu and * 0, = A - iu By imposing constraints (2.1) and (2.2) on (1.7), we obtained the integration formula (2.3) with r and s given by ,~ , , s ,, >. e {(A -u )sin(hu)-2Au cos(hu)} + 2Au \Z.A.L) r^A,u; - * r , u(A +u ) \ v» /0 10N ,, v e {A sin(hu)-u cos(hu)} + u u(A Z +u ) Near the origin, equations (2.11) and (2.12) reduce to (2.13) r(A,u) = h and (2.14) s(A,u) = ^- , and thus the integration formula (2.3) reduces to the second order Taylor series. -9- 3. EVALUATION OF STIFFNESS PARAMETERS By using the Taylors expansion of y .. = y(x +h) about x = x « 1 h -Oh t and the Maclaurin's series of e , e in the integration formula (2.3) 1 2 it is observed that the coefficients of h , h , and h vanish identically. With the view to obtain numerical estimates for the stiffness/oscillatory 3 4 parameters, we simply allow the coefficients of h and h to vanish thus yielding the following pair of nonlinear equations (3.1) (n 2 -fi 1 )f^ 1) - n 1 n 2 f t - -f£ 2) , and (3.2) -(flJ-fl^+G*)^ + ^ 2 (fi 2 -fi 1 )f^ 1) = -f^ 3) . By adopting (3.1) as a definition of £Lft f , equation (3.2) becomes (3.3) (fi 2 -fi 1 )f^ 2) - fi 1 " 2 ^ 1) = ~f^ 3) • If the set of equations (3.1) and (3.3) were to be meaningful, it is desirable that i f (2) i,(l) r r t t (3.4) det i # 0, i = l(l)m. \ I Let (3.5) X D = X ^ 2 - l Q , i = l(l)m; and (3.6) i E = 1 fi 1 1 fi 2 , i = l(l)m Equations (3.1) and (3.3) can now be expressed as m pairs of linear equations (3.7) -10- ^Vf (2) - Vf (1) =-V 3) , t t t ' for i = l(l)m. l D i f (D . V ft = -i f <2) , These m pairs of equations can be readily solved for D and E by the Cramer's rule to give i f i f (3) _ i f (l)i f (2) (3 ' 8) lD ■ i f (l)i f (l) . I i f (2) - 1=1(1 > m t t t t and i f (l)i f (3) _ i f (2)i f (2) < 3 ' 9 > ±E ■ i f (i)i f (i) _ i f fc i f (2) > f"»-- t t t t The numerical values obtained for D and E can be substituted in equations (3.5) and (3.6) to generate the oscillatory/stiffness parameters as (3.10) V = \[-\ ♦ /(V+^E)], and (3.ii) x si = x n + X D. The complex interpolation formula (1.7) is adopted if in (3.10), the following relationship holds (3.12) V < 4 X E -11- 4. STABILITY CONSIDERATIONS We now apply the integration formula (2.9) to the scalar test equation (4.1) y' = Ay + g where X is a complex constant with negative real part and g is a real constant. The numerical solution is given as (4.2) y t+1 = pU,fi 2 ,h)y t + gh where the characteristic equation p(A,fl„,h) is given by 2 "S h 2 (4.3) p(X,fl 2 ,h) = 1 + Ah + X (SI h+e -l)/« 2 Equation (3.11) gives the following relationship: (4.4) tt 2 = X D = D If we use equation (4.1) in (3.8), we readily obtain (4.5) D = -X = ft 2 This in equation (4.3) gives (4.6) P(X,-X,h) = e Ah . We further consider the application of the integration formula (2.3) to test problem (4.1) for the case when the stiffness/oscillatory parameter is imaginary. Thus, r(X,u) and s(X,u) used in the two integration formulas (2.3) to the test problem (4.1) as specified by ,, ,v fy s sin(hu) (4.7) r(X,u) = J — '- , and /, Q \ /-, n l-cos(hu) (4.8) s(X,u) = 2 . u -12- Let (4.9) X = iz, i 2 = -1. The resultant integration formula is given by (4.10) y t+1 = q(X,u,h)y t + gh, where the characteristic equation can be obtained as (4.11) q(X,u,h) = 1 + (l-cos(hu))X 2 + sin(hu)X < u The application of equation (4.1) and (4.9) in (3.8) yields (4.12) D - -iz = -X. ■k = n . This in (4.11) gives (4.13) q(u,u,h) = 1 - (l-cos(hu))+i sin(hu) ihu = e Xh = e Equations (4.6) and (4.13) together with definitions 1 and 2 thus establish the L-stability and exponential fitting of the proposed integration formulas. -13- 5. LOCAL TRUNCATION ERROR We now associate with the integration formula (2.3) the operator V[y(x),h] specified as (5.1) V(y(x),h) = y(x+h) - y(x) + (Sl^-Qf )f (x,y) - ( Y -to )f (1) (x,y) for an arbitrary function y(x) G c (S). The local truncation error T at x = x i is hence given as V[y(x ),h] where y(x ) denotes the solution to the initial value problem (1.5). By using the Taylors expansion of V[y(x),h] about x = x with the localizing assumption that there is no previous error (i.e., y = y(x )), the truncation error for the integration formula (2.3) with constraints (3.1) and (3.3) can be derived at (5.2) T t+1 = yj-{(fi 1 +^ 2 )" 1 [(n i +fi 2 )f^ 4) - ^(f t (1) +fi 2 f t ) + ^ 2 (f^ 1) -fi 1 f t )}+ 0(h 6 ). The corresponding truncation formula for the integration formula (2.9) is given by (5.3) T t+1 =^ (f^-^a)) +0(h 6 }) while the truncation error when r(A,u) and s(X,u) are specified by equations (2.11) and (2.12) is (5.4) T t+1 = 3y [f^ 4) + (A 2 +u 2 )(3A 2 -u 2 )f t + 4X (u 2 -A 2 )f^ 1) ] + 0(h 6 ). From equations (5.3) and (5.4) we deduce that all the proposed explicit integration formulas are of fixed order four. -14- 6. EXTENSION OF NONLINEAR SCHEME We recall the implicit numerical integration formula proposed in Fatunla [10] : (f, -n ,J S+1 ] ■- „t S ] TT T [s] n-1 r„[s] f [s] , ft - , (6.D y t+1 - y t+1 - [i- J t+1 l [y t+1 - f t+ r G t ] ' s " °' l ' 2 '" [si where J. +1 denotes the Jacobian specified as (6 2) J [s] = — (x v [s] ) f [s] = f(x v [s] ) and (6.3) G t = y t - (Y-HT)f t . The components of 9, y, a are obtained as ln h . - ln .Ji X fi (1-e ) + fl.. (1-e 2 ) (6.4) S- — : ~ , i = l(l)m, fi ft (e -e ) -*Q h ln h (6.5) \=— (1 " e , } , i = l(l)m, - ft„h fl n h ^ 1 (e -e ) and X ft h _1 ^2 h (6.6) ^ = — P^ : *- , i = l(l)m. -Vh Vh ft~(e -e ) The corresponding truncation error for the integration formula (6.1) was obtained as (6.7) T t+1 = j^Q {9f t 4) + t ( fi 2" fi l )f t 3> + 5 (" 1 -^ 2 +fi 2 )f t 2) + [5n^«2 " 9n n (n^-n n -k^]^ + °( h7 )- ■ -15- In the case when the oscillatory/stiffness parameters have complex values, the components of and (y+<$) in equation (6.1) are replaced by re o\ i Q _ e [u cos(hu) - X sin(hu)] + u _ ,>.» ( } ^7X~2~~7~r~r — ■ " () e (X +u ) sxn(hu) ,, nx i, . ~ N ue - [X sin(hu) + u cos(hu)] (6.9) (Y+6) = 2 — 2 (X +u ) sin(hu) with the truncation error (6.7) replaced by (6.10) T t+1 = -jjQ {9f^ 4) - lOXf < 3) + 5(3X 2 -u 2 )f^ 2) + 36(X 2 -u 2 )f^ 1) + 4(X 2 +u 2 )(8X 2 -u 2 )f t } + 0(h ? ) . -16- NUMERICAL EXAMPLES Example 1 We first consider the linear system: r -.1 -49.9 (7.1) y' = I -50 L 70 -120J y(0) =! l in the interval £ x _< 15. The eigenvalues of the Jacobian of (7.1) are ^ = -0.1, ^ ? = -50, ^„ = -120. The theoretical solution is given as 1 , x -O.lx , -50x y(x) = e + e 2 , , -50x y(x) = e 3 , s -50x , -120x y(x) = e + e The new explicit scheme performs better than the DIFSUB [13] and the blended DIFSUB [22]. Details of the numerical results are given in Table 7.1. TABLE 7.1. Numerical Results for Example 1. ., . , Max _ Function Back LU Accurate Method _ , Steps _ -. _ - _ _. : .. Order Calls Solves Decomp Digits DIFSUB BLENDED DIFSUB EXPLICIT SCHEME 6 353 1024 969 18 11 234 595 1056 22 75 75 9.3 10.4 12.5 -17- Example 2 We further consider the linear problem: (7.2) ~-10 100 on -100 -10 -4 y. = -1 -.5 -.1- y.(0) =1, i = 1(1)6 in the range <_ t < 20. The eigenvalues of the Jacobian are A _ = -10 ± lOOi, ■•-» * A- = -4, A, = -1, A,. = -.5, and A, = -.1. This problem is particularly troublesome for both DIFSUB and the blended DIBSUB as can be seen from Table 7.2, TABLE 7.2. Numerical Results for Example 2. G Max Func Back Order Eval Solves LU Accurate Decomp Digits DIFSUB io" 2 (4) (1001) (3002) (2959) (7) (1.1) "too much work" ; integration abandone .d at x = 8.3 blended DIFSUB io" 10 (12) (1001) (2917) (5580) (21) (10.3) "too much work"; integration abandone d at t = 5.1 Explicit Scheme unspec- ified 4 200 200 - 14.2 -18- Example 3 We now consider the initial value problem of Liniger and Willoughby [18]: (7.3a) y' = '-2000 iooo\ /i jo) y + | , y(0) = \ ; \°/ in the interval < x < 5. This problem is contained in the test batch recommended by Bjurel, et al. [4], and was also considered in Lambert [15]. In addition, it is solved with the blended DIFSUB [22]. The eigenvalue of the Jacob ian to (7.3a) are A = -2000.5 and A„ = -0.5 thus yielding a stiffness ratio 4001. The stiffness matrices which have constant entries were obtained as -0.499875 -0.499875 2000.5 i 2000.5 (7.3b) n - fi 2 = / Table (7.3) gives the numerical results for the explicit scheme, Lambert's scheme and the blended DIFSUB. The global relative error for the blended DIFSUB is 0.743400 D-5 while the global relative error for the explicit scheme is 0.5746777037 E-05. The computational cost and function evaluation on IBM 360/75 for the blended DIFSUB are respectively $1.84 and 107 while that of the explicit scheme is $0.77 and 10, respectively. The global errors were computed as (7.3c) (7.3d) m e = max E 0* •H C r- CM CO m v3- CO o r-~ 00 CO •U O m as -d- OS rH m m CO CO r^ a) -h H 00 CM CM rH r^ CM vO a\ rH }-l -M o 3 CM CO in vO r^ r- 00 00 00 OS a; rH ,£ o H CO rH 4J O r*« OS CM vo rH CU d >* 00 CO rH ,0 i—i rH 00 CM CN rH r- CM vO as rH g LO ll CO rH CM CO m vO r^ r-~ 00 00 00 as ►J •— ' JC T3 VO m m CT\ CM rH as m rH O CO r-» rH co vO vO -tf rH a> Pn CM rH M CM CM CO m vO r^- r^ 00 00 00 as PQ O "— ' 4J •h lh vO CO r» rH m m 00 •H O O CO r-^ rH CO vO vO ov CTi 0) rH X o H CO rH u O rH r- m CTi ,0 r— i O ov vO rH m 00 rH CO ~* m 6 m II CO rH vO VO r^ 00 00 00 o> as as as (J •— ' 43 CM CM rC CUvO 00 o CM O CM CO CO vO o- rH T3 .— ■ rH 1 co m vO OS vO co O CM CM o> cu pa X O o vO co m VO 00 CO CM r- 00 13 !=> CO rH rH o\ vO rH m 00 rH CO as as rH tH CO PQ Q > 4J •U •H m 00 m m CM co - 2.4 ± 2.8i -+ -0.052 ± 8.8i ■> -2.0 -*• 9.5i ■> -5.9 ■> 4.5i ■* -2.0 and -12 -v 0.050 and -15 + 1.1 and -3.4. The numerical computation was effected with both the explicit and the implicit schemes using uniform mesh sizes h = 0.2, 0.15, 0.1, 0.05, 0.025 and 0.0125. The same problem was solved with the blended -3 DIFSUB using an initial mesh size h = 10 as suggested in Enright, et al. [6]. From Table 7.4 the new schemes compare favorably with the blended DIFSUB. -21- TABLE 7.4. Numerical Results f or Example 4. EXPLICIT No. of Fn. Eval. h yd) y(2) 5 0.2000 1.8716065 -0.14358810 7 0.1500 1.8711045 -0.14574713 10 0.1000 1.8705973 -0.14610294 20 0.0500 1.8694380 -0.14823599 40 0.0250 1.8694389 -0.14823587 80 0.0125 1.8694388 -0.14823588 Blended DIFSUB No. of Fn. Eval. h =10"3 e =io -6 y(D 1.86944 y(2) -0.148236 IMPLICIT No. of Fn. Eval. h yd) y(2) - 0.2000 * * - 0.1500 * * 42 0.1000 1.8693953 -0.14824187 81 0.0500 1.8694357 -0.14823631 161 0.0250 1.8694387 -0.14823589 321 0.0125 1.8694389 -0.14823587 * No convergence -22- Example 5 (7.5a) /-icf 5 100 \ y' = -5 y(0) = < X < lOlT. 1 \-100 -10 This is the model problem (1.2) with £ = 10 and w = 100 whose theoretical solution is _ in ~5 sin lOOx (7.5b) y(x) = e iU X \ cos lOOx The numerical results in Table (7.5) are obtained with a uniform mesh size h = tt/20 and only one evaluation of the stiffness/oscillatory parameters was obtained and given by -10" 5 -1000i o\ ?'-10 5 +1000i o\ °1 = 1 -5 I' "2 = I -5 \0 -10 -10001/ \0 -10 +1000i TABLE 7.5. Numerical Results for Example 5. x = 10 , uniform h = tt/20 X in 12 X t 10 x T t+1 12 2 10 ±Z x T 7T 0.16918 0.00949 2tt 0.02447 0.02233 3tt 0.29342 0.03504 4tt 0.61126 0.04402 5tt 0.92906 0.05704 6tt 1.2468 0.07060 7T7 1.56451 0.08478 8tt 0.09729 0.09991 9tt 1.29056 1.06429 10tt 1.60815 0.12150 t+1 one evaluation of oscillatory/stiffness parameters -23- Example 6 (Stiefel and Bettis [24]) We finally consider the nearly periodic initial value problem which was earlier studied by Stiefel-Bettis [24] and Lambert-Watson [16]: (7.6a) y M + y' = 0.001e ix , y(l) = 1, y'(0) = 0.99951, y e R 1 whose theoretical solution is y(x) = u(x) + iv(x) , u, v G R (7.6b) <( u(x) = cos x + 0.0005x sin x j v(x) = sin x - 0.0005 cos x Equations (7.6b) represent motion on a perturbation of a circular orbit in the complex plane in which the point y(x) spirals slowly outwards such that its distance from the origin at any time x is given as (7.6c) t(x) = /(u 2 (x) + v 2 (x)) The initial value problem (7.6a) can be expressed as |V = 2 y Vo) = i V = -S + 0.001 cos x , 2 y(0) = (7.6d) { | V = y , y(0) = o V = " 3 y + 0.001 sin x , 4 y(0) = 0.9995 The system (7.6d) was solved with the explicit formulas in the range <_ x _< 40tt which corresponds to 20 orbits of the point y(x). The integration was performed using uniform mesh sized h = tt/4, tt/5, tt/6, tt/9 , and tt/12. Two sets of numerical results were generated — in the first set, the oscillatory/stiffness parameters are obtained once at the first step of integration while in the second set, the oscillatory parameters are evaluated at every step of integration. -24- The same problem was solved with the symmetric multistep method of Lambert and Watson [16] as well as the StBrmer-Cowell five-step multistep formula [24] (both of order 6): (7.6e) y t+5 - 2^ + y^ - JL (18f t+5 + 209^ + kt^ + l« t+2 , - 6f t + l +f t >- The exact distance from the origin x(x) is given by (7.6c) -jo 3 2 and the approximate distance is t = /( y + y ) at x = 40tt. All the solutions generated by the new scheme spiral outward in agreement with the theoretical solution as well as Lambert's scheme whilst the first three values generated by Stormer-Cowell scheme spiral inward. From the Tables (7.6a, b and c) which respectively show x, |t(x)-t| and |y(x)-y| = /[ ( y(x)- y) + ( y(x)- y) ] , we see that despite the fact that the new scheme is of lower order, yet it is more accurate than both the symmetric multistep method as well as the five order Stormer-Cowell multistep scheme. -25- TABLE 7.6a. x f = 40ir , t (x ) = 1.001972 h Stbrmer-Cowell Symmetric Explicit one evaluation repeated eval. of parameters of parameters tt/4 0.965645 1.003067 1.002311 1.001972 tt/5 0.993734 1.002217 1.002205 1.001972 tt/6 0.999596 1.002047 1.002140 1.001972 tt/9 1.001829 1.001978 1.002050 1.001972 tt/12 1.001953 1.001973 1.002016 1.001972 TABLE 7.6b, x f = 40rr, x (x f ) = 1.001972 h 10 x|t(x)-t| 10 9 x|t(x)-t| Stormer-Cowell Symmetric Explicit one evaluation repeated eval. of parameters of parameters tt/4 36 327 1 095 339 204 tt/5 8 238 245 233 66 tt/6 2 376 75 167 26 tt/9 143 6 78 3 tt/12 19 1 44 TABLE 7.6c. x f = 40tt, u(x ) = 1, v(x f ) = 0.062832 h 10 6 x y(x J-yJ 10 9 x y(x f )-y | Stormer-Cowell Symmetric^ Expli* one evaluation of parameters :it repeated eval. of parameters TT/4 48 014 31 272 389 384 tt/5 13 136 7 300 252 159 tt/6 4 494 2 303 176 77 tt/9 405 188 79 15 tt/12 73 33 45 5 -26- 8. CONCLUDING REMARKS The proposed explicit scheme (2.3) and (2.9) is considered to be more efficient and accurate than the DIFSUB and the blended DIFSUB for linear stiff systems of ordinary differential equations. It is equally efficient for highly oscillatory systems as it is capable of admitting fairly large mesh size and still maintains high degree of accuracy. The major drawback is the need to generate higher order derivative but automatic generation of higher order derivatives is practicable for an extensive range of problems. -27- LIST OF REFERENCES [1] Amdursky, V. and A. Ziv (1974). On numerical treatment of stiff, highly oscillatory systems, IBM Technical Report 015 , IBM Israel Scientific Center. [2] Amdursky, V. and A. Ziv (1975). On the numerical solution of stiff linear systems of the oscillatory type, IBM Technical Report 032 , IBM Israel Scientific Center. [3] Amdursky, V. and A. Ziv (1976). The numerical treatment of linear highly oscillatory ODE systems by reduction to nonoscillatory types, IBM Technical Report 039 , IBM Israel Scientific Center. [4] Bjurel, G. , Dahlquist, G. , Lindberg, B., Linden, S. and L. Oden (1970). Survey of stiff ordinary differential equations, Computer Science Report NA 70.11 , Royal Institute of Technology, Stockholm, Sweden. [5] Dahlquist, G. (1963). A special stability problem for linear multistep methods, BIT 3, pp. 27-43. [6] Enright, W. H. (1974). Second derivative multistep methods for stiff ordinary differential equations, SIAM J. Numer. Anal. 11 (2), pp. 321-331. [7] Enright, W. H. , Hull, T. E. and B. Lindberg (1975). Comparing numerical methods for stiff systems of ODEs, BIT 15 , pp. 10-48. [8] Fatunla, S. 0. (1976). A new algorithm for numerical solutions of ODEs. Journ. Computer and Mathematics with Applications 2 (3/4), pp. 247-253. [9] Fatunla, S. 0. (1977). A variable order one step scheme for numerical solution of ODEs, Journ. Computer and Mathematics with Applications (to appear) . [10] Fatunla, S. 0., An implicity two point numerical integration formula for linear and nonlinear stiff systems of ordinary differential equations, J. Mathematics of Computation (to appear). [11] Gautschi, W. (1961). Numerical integration of ODEs based on trigonometric polynomial, Num. Math. 3, pp. 381-397. [12] Gear, C. W. (1968). The automatic integration of stiff ordinary differential equations, Proceedings IFIP Congress _1, North- Holland, Amsterdam, pp. 187-194. [13] Gear, C. W. (1971). Algorithm 407: DIFSUB for solution of ordinary differential equations, Comm. of the ACM 14, pp. 185-190. -28- [14] Jackson, L. W. and S. K. Kenue (1974). A fourth order exponentially fitted method, SIAM J. Numer. Anal. 11 (5), pp. 965-978. [15] Lambert, J. D. (1973). Nonlinear methods for stiff systems of ordinary differential equations, Proceedings Conference on Numerical Solution of Ordinary Differential Equations 363 , University of Dundee, pp. 75-88. [16] Lambert, J. D. and I. A. Watson (1976). Symmetric multistep methods for periodic initial value problems, J. Inst. Math. Appl. 18 , pp. 189-202. [17] Lindberg, B. (1971). On smoothing and extrapolation for the trapezoidal rule, BIT 11 , pp. 29-52. [18] Liniger, W. and R. A. Willoughby (1969). Efficient numerical integration of stiff systems of ordinary differential equations, IBM Research Report RC 1970 , IBM Yorktown Heights, New York. [19] Miranker, W. L. and G. Wahba (1976). An averaging method for the stiff highly oscillatory problems, Math. Comp. 30 , pp. 383-399. [20] Miranker, W. L. , van Veldhuizen, M. and G. Wahba (1976). Two methods for the stiff highly oscillatory problem, P roceedings of nu meric a l an alysis conference held in Du blin , Ed. J. Miller (to appear) . [21] Miranker, W. L. and M. van Veldhuizen (1977). The method of envelopes, IBM Research Report RC 6391 (#27537) , Mathematics Division, IBM Yorktown Heights, New York. [22] Skeel, R. D. and A. K. Kong (1976). Blended linear multistep methods, UIUCDCS-R- 76-800 , Department of Computer Science, University of Illinois, Urbana, Illinois. [23] Snider, A. D. and G. L. Flemming (1974). Approximation by Aliasing with applications to "certain stiff differential equtions," Math. Comp. 28, pp. 465-473. [24] Stiefel, E. and D. G. Bettel (1969). Stabilization of Cowell's method, Numer. Math. 13, pp. 154-175. mAEC-427 U.S. ATOMIC ENERGY COMMISSION i 6 / 6 ?' UNIVERSITY-TYPE CONTRACTOR'S RECOMMENDATION FOR DISPOSITION OF SCIENTIFIC AND TECHNICAL DOCUMENT ( See Instructions on Revnm Side ) AEC REPORT NO. COO-2383-0044 2. TITLE NUMERICAL INTEGRATORS FOR STIFF AND HIGHLY OSCILLATORY DIFFERENTIAL EQUATIONS TYPE OF DOCUMENT (Check one): [x] a. Scientific and technical report l~l 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): PH a. AEC's normal announcement and distribution procedures may be followed. (~l b. Make available only within AEC and to AEC contractors and other U.S. Government agencies and their contractors. |~1 c. Make no announcement or distribution. REASON FOR RECOMMENDED RESTRICTIONS: 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^TL 61801 Signature ^-JA- ^ C ^ f Date October 1977 FOR AEC USE ONLY CONTRACT ADMINISTRATOR'S COMMENTS, IF ANY, ON ABOVE ANNOUNCEMENT AND DISTRIBUTION .RECOMMENDATION: PATENT CLEARANCE: □ 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. 3GRAPHIC DATA 1. Report No. UIUCDCS-R-77-894 3. Recipient's Accession No. 5. Report Date October 1977 and Subtitle NUMERICAL INTEGRATORS FOR STIFF AND HIGHLY OSCILLATORY DIFFERENTIAL EQUATIONS 6. ior(s) Simeon Ola Fatunla 8- Performing Organization Rept. No 'UIUCDCS-R-77-894 orming Organization Name and Address Department of Computer Science University of Illinois Urbana, IL 61801 10. Project/Task/Work Unit No. 11. Contract /Grant No. US ERDA/EY-76-S-02-2383 nsoring Organization Name and Address US Energy Research and Development Administration 9800 South Cass Avenue Argonna, IL 60439 13. Type of Report & Period Covered Technical 14. plementary Notes Some L-stable fourth order explicit one step numerical integration formulas which require no matrix inversion are proposed to cope effectively with systems of ordinary differential equations with large Lipschitz constants (including those having highly oscillatory solutions). The implicit integration procedure proposed in Fatunla [10] is further developed to handle a larger class of stiff systems as well as those with highly oscillatory solutions. The same pair of nonlinear equations as in [10] is solved for the stiffness/oscillatory parameters. However, the nonlinear systems are transformed into linear forms and an efficient computational procedure is developed to obtain these parameters. The new schemes compare favorably with the backward differentiation formula (DIFSUB) of Gear [12, 13] and the blended linear multistep methods of Skeel and Kong [22], and the symmetric multistep methods of Lambert and Watson [16]. y Words and Document Analysis. 17a. Descriptors ordinary differential equations, L-stable large Lipschitz constants stiffness parameters oscillatory parameters exponential fitting entifiers/Open-Ended Terms 3SATI Field/Group liability Statement unlimited TIS-3S 110-70) 19. Security Class (This Report) UNCLASSIFIED 20. Security Class (This Page UNCLASSIFIED 21- No. of Pages 28 22. Price USCOMM-DC 40329-P7 1 MH 2 5 1978 %, I MBti