, 143 min read

Text: A Stiffly Stable Integration Process Using Cyclic Composite Methods

JOEL MARVIN TENDLER

B.E., The Cooper Union, 1964

DISSERTATION

Submitted in partial fulfillment of the requirements for i the degree of Doctor of Philosophy in the Graduate’ School of Syracuse University, May, 1973 , approved eden, £Luker= Date 2¢ Fobra eu 1472

A STIFFLY STABLE INTEGRATION PROCESS USING CYCLIC COMPOSITE METHODS

by

JO MARVIN TENDLER

B.E., The Cooper Union, 1964

ABSTRACT OF DISSERTATION ' Submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy in the Graduate School i of Syracuse University, May, 1973 ) approved Lodoru ALI pate 2( February /%73

ABSTRACT

Let y = fly, t) with ALY = Yo possess a solu ion 4 y(t) for t > t,. Set t = t, + nh and let y denote the approximate solution of ye) defined by the composite | linear multistep method 1 | L I (25,5 Yneak-143 h8; 5 Yngex-143 - 9 j=-k+1 i i=l, ..., 2, I with Yn = fly, t;) and n=1, 2, ... . Asymptotic stability properties of the method as reflected in the behavior of i ypasn-+= are examined; in particular, Ala)-stability and stiff stability are treated. It is shown that, when the I method is used to numerically solve y = qy, q constant, 4 [ with a fixed integration step-size h, the characteristic polynomial, a polynomial in two variables, { and A = gh, , [ may be examined to determine the asymptotic behavior of the given method. The Lambda and Zeta Loci determined by the characteristic equation are introduced and their properties examined. Specifically, the Lambda Locus is the A-plane image of the {-plane unit circle under the mappings . [ Ay = £02), i=l, ..., 2, and determines the stability regions for the method. Similarly, the image of the A-plane imaginary axis in the g-plane under the mappings ty = 95) ° . j=1, ..., k, determines the Zeta Locus, a useful tool in

checking if a method is A-stable. Note: 1X; = £, (5) for each i is a zero of the characteristic polynomial viewed as a polynomial in A; similarly, t5 = 93) for each § (3 a zero of the characteristic polynomial viewed as a polynomial in {. A modified Zeta Locus is also discussed and it is , shown that it can be used to check if a method is stiffly stable. A computer program to plot both the Lambda and Zeta ) Loci for a given method is described. Using an interactive l version of this program new and efficient cyclic composite linear multistep methods, that is, methods given as

I 05,5 ¥ngek-145 ~ 0 I B55 Ynpex-145 = 27 j==k+i j=1

i=l, ..., 2%, . of orders 3 through 7 are derived. These new methods are stiffly stable and exhibit better stability properties than ‘ the backward differentiation formulas of corresponding order. [ A new variable step, variable order integration algorithm incorporating these new methods is presented. This new f algorithm is ideally suited for the numerical integration of stiff systems of first order differential equations where- . in the Jacobian matrix fy has complex eigenvalues. Numerical results to support this claim are presented, including a comparison of the new algorithm with Gear's widely accepted ’ computer program DIFSUB and with a third algorithm based on

the backward differentiation formulas. For cone of the test problems solved an improvement of better than an order of i magnitude in computer program execution time was note | using

the new algorithm over results obtained using Gear's algorithm DIFSUB. ad :

[ iv

’ ¢ ~ A STIFFLY STABLE INTEGRATION PROCESS USING CYCLIC COMPOSITE METHODS

JOEL MARVIN TENDLER B.E., The Cooper Union, 1964 DISSERTATION Submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy in the Graduate School of Syracuse University, May, 1973 “ | Approved Led NITE 8 1 Date __2 Febriorsx 1472

AJ

£ \

[1] ) To the two Philips in my life— one of whom has affected me i as I pray I can the other }

ii |

PREFACE With the advent of larger and faster computers new problems have been opened to solution. More complex Sys- tems, described by systems of ordinary differential egqua- i tions, are now being solved than ever before. This places } a bigger burden on numerical integration subroutines as they must be capable of efficiently and accurately supply- g i ing answers. In the case of stiff systems, that is, sys- tems with both fast and slow transients, the numerical i integration process can be both time consuming and expen- sive. This work was undertaken to help alleviate this prcblem. A new integration algorithm is presented that 1 does not suffer from many of the shortcomings of presently available algorithms. In the very least, the techniques 1 developed to arrive at this new algorithm should provide I others with the insight and tools for pursuing this field. ‘ I would like to express my sincerest. appreciation to my advisor, Professor Theodore A. pickart, for his guidance and supervision during the course of this work. i His insight and knowledge of this research area has been invaluable I am also grateful to Mr. William B. Rubin, . 1 a fellow student, for the many enlightening conversations we had. Support provided by the Syracuse University Research Corporation in the form of unlimited free computer time and f iid

. L :

’ (

secretarial services has proven invaluable. I would like to H especially thank Mr. James D. Rodems for providing me with the opportunity of pursuing this work while still remaining t as an employee, and Mrs. Christine Peters and Miss Janet H sinclair for typing and overseeing this manuscript. The financial support provided by the National Science Founda- tion, to myself in the form of a traineeship and to my advisor under Grant GK-23010, and by Syracuse University is also gratefully acknowledged. 1 Finally, and perhaps foremost, I am indebted to my wife Debbie for her understanding and encouragement during i my studies.

I i.

ff iv

. pe

CONTENTS

A te. LIST OF ILLUSTRATIONS . . « + « + = * = * = ° s 0 09 no IE ¥ LIST OF TABLES . + « « + + = = * Cc eee oo ® ob ANE i CHAPTER I CHSOODUCEEON + o « co co oc cco vo soo 13 CHAPTER II PRELIMINARIES . . + « + + = + + = = = = & ¢ 6 1 2.1 Problem Statement . . . EEE EE 6 2.2 Linear Multistep Methods . . . + «+ » « = = = 8 | PEE eso soc cocsance B i 2.4 A(a)-Stability . « + « « « « » oo = © = = = ¢ 13 2.5 Runge-Kutta MethodS . « « « s + » = = = » = 14 2.6 Composite Linear Multistep Methods . . . + - 17 2.7 The Characteristic Equation . . « + + o oo » 21 CHAPTER III STABILITY BY GRAPHICAL METHODS 5 a aw [] 3.1 Introduction . . . + + + os + oc e000 26 3.2 The Lambda Plot . . + « «+ « = « = = = = = * ¢ 29 i BC od 3.84 The Seta PIO , o + + os oo oo oo 0 0&0 0s 43 3.5 A Plotting Program . . . « « + = + = * = = * 51 CHAPTER IV A NEW INTEGRATION PROCESS , . + « « « » + oo = 67 " 4.1 Free Parameters . . » » + +» + + + + = * = °° 67 i | 4.2 Some New Cyclic Composite Linear , Multistep Methods . . + « + + + = = = = = = °¢ 69

wv Fe

CHAPTER V A NEW VARIABLE ORDER/VARIABLE STEP-SIZE ’ 1 INTEGRATION ALGORITHM cs so » hE WSS. 85 5.1 Introduction . . + + c+ c+ "C7 REE 5.2 Computer Program organization and Logic . + « 85 5.2.1 Solution of the Implicit Equations . . 1] 5.2.2 Error Control . . » + + + * °° vu & 9 5.2.3 Step-Size and order Selection Logic .102 5.2.4 Start-Up procedure . . + + + + + + * .111 ] 5.3 Sample Problems and Comparisons . . « + «+ = L112 CHAPTER VI CONCLUSIONS AND FUTURE WORK . . +. « « =» + = .130 APPENDIX A COMPT'TER PROGRAM LISTING FOR BASIC PLOTTING PROGRAM as « oo a oo oh 1] APPENDIX B COMPUTER PROGRAM LISTING FOR INTERACTIVE VERSION OF SUBROUTINE ACCEPT . . . «+ + = = .146 I APPENDIX C FLOWCHART AND COMPUTER PROGRAM LISTING 4 Cp resaau EI I CC EEE RR ER RR 3 EE EE RR ER cate { vi

& [*] .

. » a

LIST OF ILLUSTRATIONS

FIGURE 3.1 pefinition of the Contour Batt 20 Salat . 4 FIGURE 3.2 Plotting Program Organization =. + + * * ° eo . Kd FIGURE 3.3 Main program—Initialization procedure - - + 99 | FIGURE 3.4 Main pProgram—Start-up procedure . + + * * * 58 . FIGURE 3.5 Main Program—Locus computation + + + + °° 64 FIGURE 4.1 Lambda Locus for Order 1 and 2 Methods - . B80 FIGURE 4.2 Comparison of Lambda Loci Between New Methods and BDF for orders 3 through 6... 81 ps i FIGURE 4.3 Lambda Locus for New Order 7 Method +... » + B83 FIGURE 5.1 organizational Flowchart for Subroutine I sew rrr Tr EERE ER FIGURE 5.2 Lambda Loci and Ray pefined by (5.49) ep me wr REE Ray t FIGURE B.1l organizational Chart for Interactive ACCEPT. 147 FIGURE C.1 Detailed Flowchart for DIFJMI . = « «= =~ 162 .

¥ o a

LIST OF TABLES

TABLE 4.1 Linearly Independent Multistep Formulae i OE D3, coco Do oo oc cco ssc soo UB TABLE 4.2 New Cyclic Composite Multistep Methods . . . 76 TABLE 4.3 Comparison of Widlund Wedge Angles . . . - - 78 TABLE 4.4 Comparison of y-values . . + + + + + + + + * 78 i TABLE 5.1 piscretization ro Constants for Methods Used in DIFJMT . . . . + «+ « « + = = 103 TABLE 5.2 Maximum Allowable Step-Size Changes . . . . . 108 . TABLE 5.3 EE EE RI TABLE 5.4 EE i TABLE 5.5 III I IO ER ar TABLE 5.6 se ss sso soe BF TABLE 5.7 2 = - «>> scn-10 TABLE 5.8 Conditions Used for the Numerical Solutions . 120

I ¢

viii 3

« i ~~ CHAPTER I : INTRODUCTION We are concerned with the numerical integration of a % stiff system of first order ordinary differential equations* f ye) = £ly(e), ©), (1.12) given the initial conditions i y(t) = Yor (1.1b) By stiff we mean a system with widely separated eigenvalues I of the Jacobian matrix J = fy of (1.1) evaluated along a i particular solution y(t) = g(t) . Note: stiff systems are ' encountered in many disciplines, €.g., reactor calculations, circuit analysis, chemical kinetics, etc. JOR —— *1+ is assumed that (1.1) possesses a unique solution y(t) for t 2 t,. Such will be the case if fly), t) satisfies i a Lipshitz condition, i.e., there exists a constant L > 0 such that i Ely, &) - £(z, ell <x lly - zl for all t > 0 and y, ze "of I AJ

¢ TN -

¥ 2 Conventional numerical integration processes when B applied to stiff systems limit the integration step-size, h, to a value determined by the largest eigenvalue or the 1 Jacobian matrix, even after the transient effects due to the larger eigenvalues have decayed. This unnecessarily ) ! causes the numerical integration process to become both time consuming and expensive. Thus there is a need for numerical integration processes which perform well when ’ | used to solve stiff systems. A detailed discussion of the history and state-of- the-art as of 1970 of the methods used for the numerical ] solution of stiff systems can be found in Bjurel et al & {BJU70). A more recent bibliography covering the field can [ be found in [ENR72]. Almost all of the previous work has been limited to classical multistep and Runge-Kutta type processes. of all [ the schemes proposed the most widely accepted process is a variable order, variable step algorithm due to Gear [GEA71a]. Gear's algorithm is based on the backward differentiation [ multistep formulas described in [HEN62]. Recent modifica- tions to Gear's algorithm by Brayton, Gustavson, and Hachtel R [BRA72a) and by Krogh [KRO72] have attempted to improve on some of the shortcomings in Gear's method, but their pro- cesses are based on the same formulas used by Gear. The performance of all three algorithms deteriorate when they '

. Vv

¢ i ~~ FA (

are used to integrate a stiff system wherein the eigenvalues of the Jacobian matrix are complex. Recently Sloate and Bickart [SLO71b] introduc 4 the concept of a composite linear multistep method wherein a 3 block of ¢ future points are simultaneously computed using § @ifferent equations. Sloate {S1071a] has reported an al- i gorithm suitable for the solution of stiff systems that is unaffected by the presence of complex eigenvalues in the ’ i Jacobian matrix. However, the algorithm consisted of solv- ing for two future points and retaining only the first of these points at each stage of the process. [| Since then Bickart, Burgess, and Sloate [BIC71], Bickart and Picel [BIC72], and Watts [WAT71] have reported i on composite linear multistep methods wherein all & future points are retained with the additional desirable property of requiring only one backpoint making the processes self- I starting and simplifying step~-size changes. These processes : are applicable to the numerical solution of stiff systems 1 and, to varying degrees, are not as sensitive to complex i eigenvalues in the Jacobian matrix as are the processes based on the backward differentiation formulas. However, [ they suffer the disadvantage of requiring the simultaneous ’ 1 solution of & coupled equations. One of the goals of the research effort reported here- |i in was to obtain a new integration process applicable to » pi .

p i ~~

[| “ 4 stiff systems that is efficient to use and that is not "too" I sensitive to complex eigenvalues in the Jacobian matrix. As composite linear multistep methods satisfy two of the th ‘ee ¥ requirements they were considered. In order to make the process efficient attention was further directed toward the class of cyclic composite linear multistep methods where- i in each of the & component formulas are solved cyclically for one of the % future points. The concept of a cyclic composite linear multistep method is not new. It was first formally i considered by Donelson and Hansen [DON71] and later by Chesler and Pierce [CHE71], and by Dyer, Pierce, Haney, and Chesler I [DYE72]. However, their interest was in orbit computation and the processes they considered are not suitable to stiff systems. 1] A search of the available literature soon revealed ’ that tools for characterizing the general class of composite linear multistep methods that would be required in obtaining the new integration processes were noticably missing. Hence a second goal of the research effort reported herzin was to obtain new techniques to be used as tools in characterizing numerical integration processes, specifically composite lin- ear multistep methods. Chapter II begins with a discussion of classical multistep numerical integration processes and procedures used to characterize them and concludes with a description B

rs 3 *

'e

of composite linear multistep methods. Scme historical i perspectives are also presented. Furthermore, since it is a key ingredient of subsequent results reported, t e B concept of the characteristic equation of a multistep process is introduced. Chapter III describes a program used to obtain stability boundaries for multistep processes i from the characteristic equation. The program described in Chapter III was used to obtain a new composite linear | multistep method; that result is reported in Chapter IV. 2 variable order, variable step-size algorithm using this new method is described in Chapter V. Results obtained with this new algorithm from some sample problems are also re- ported. Chapter VI contains conclusions about the new J t method and an outline for future research.

CHAPTER II

{ PRELIMINARIES 2.1 Problem Statement Consider the system of first order differential i equations y(t) = fly(t), t), ylty) = Yor (2.1) It is desired to obtain a sequence {y,} such that y ap- proximates y(t) - gle, + nh) . Historically predictor-corrector methods and ex- plicit Runge-Kutta type processes have been used to solve (2.1) numerically. However, when (2.1) is stiff, classi- cal methods often encounter serious problems. For stiff systems one wants to increase the step-size, once the transient effects due to smaller time constants have de- ) A cayed, in order to decrease computer program execution { time. However, this is in conflict with the small step- size required to keep the numerical processes stable. Hence one is forced to use a small step-size resulting K in unnecessarily long and expensive program execution times.

o Fay .

i ’

The phenomenon of a numerical integration process i becoming unstable even when used to integrate a stable differential equation can be exhibited by considering the forward Euler formula I Yn+l =¥5 os hy,» (2.2) I when used to integrate the differential equation y=ay. y(0) =v. (2.3) gq real and negative. The differential equation (2.3) is stable; however, when (2.2) is used to numerically inte- grate it, the solution becomes

yp = (1 + ah)" yo (2.4) which for h > =-2/q is unstable. This follows from the fact that, for h > -2/q, |¥p41! > ly,! , hence, y += as { n+w, Thus, the forward Euler method is unstable under these circumstances, even though the differential equa- tion itself is stable. Our aim is to present a new integration method which avoids the pitfall exemplified by the forward Euler method. }

¢ i - » !

2.2 Linear Multistep Methods (LMM)*

i { Linear multistep methods may be written in the form

I (ayy, = hBy¥n,s) = 0 (2.5) j=-k+1 The LMM as given by (2.5) may be used ~o0 numeri- cally integrate (2.1) by setting Y: = fly .ty)- As the differential equation requires an initial condition so does (2.5). In fact (2.5) requires Yore+= gal as start- ing values. As information at k different points in time is required it is called a k-step method. [ Next consider

Lly(t)), hl = I logyleg,y) - hBy¥ (tn) 1} (2.6) j=-k+1 Under the assumption that y(t) in (2.6) possesses a con- vergent Taylor series at t = the (2.6) can be written as

N | Liye), ml = I oeply®ien, 2.7

i ———

*For notational simplicity only the scalar form of (2.1) will be considered. The results extend, in a straight- forward manner, to the vector case.

« i“ ~~

where

1.4 1 i-1 eg= I qr¥ey-gmmrd fy * j=-k+1 The basis now exists for pefinition 2.1: The linear multistep method (2.5) is said to be of order p if C; = 0 for i=0, «ces P but Cp+1 # 0, where Cs is given by (2.8). The order of a method is related to the local discretiza- tion error as follows: If

  1. the multistep method is of order p;

. |

‘ Ww ~~ A 5 {

y(t) = gq y(t) (2.9) with any fixed positive h where Re{g} < 0. Substituting (2.9) into (2.5) yields

. = qhB, . = 0. (2.10 j==k+1 The solution of (2.10) is given as

yo= I 2% (2.11)

where, with gh = ), the t's are the roots* of

I (ay - 28g) F143 =o, (2.12) j==k+1 ’ JEN—— #The form of the solution as given by (2.11) assumes that the zeros of (2.12) are all distinct. In general the solution is given as [HEN64]

j-1,n ' ( EE) ag gg (2.11%) i=1 j=1 where the g's are the s distinct roots of (2.12) with multiplicity m; respectively, and the a;4's are functions ’ of the initial conditions. As the primary concern of the f present discussion is with A-stability it suffices to as-= Pume that there are k distinct zeros of (2.12).

|! -

and the a's are functions of the initial conditions. (Equation (2.12) is called the characteristic equation of the method.*) i As a result of (2.11), if an LMM is to be A-stable, i then Izy! <1, i=l, ..., k, for all values of ) such that Re{)} < 0. i The forward Euler method given by (2.2) is not ) A-stable. Its characteristic equation is g~-(1L+21)=0, (2.13a) or | 4 | zg =1 + gh. (2.13b) It is readily verified that for Reig} < 0, zg] > 1 for 11+] > 1, (2.142) or h>- ela) : (2.14b)

——— 1 #*The characteristic equation is discussed in detail in ’ | Section 2.7. |

¥ sd

‘ » ~~

H 12 On the other hand, the backward Euler method, given i as Yos1 ™ Yn * Dna (2.15) is A-stable. Its characteristic equation 1s gb -2) -1=0, (2.16a) or

23 — & (2.16b)

For Re{g} < 0 it is seen that |g] < 1. [ The backward Euler method is implicit* and hence | more difficult to solve as ¥, ., appears on both sides of the equation. Dahlquist (DAH63] has proven that no ex- ! plicit* LMM can be A-stable. He further showed that if (2.5) is to be A-stable the order of the integration formula cannot exceed 2.

As high order integration processes that do not restrict the step-size due to stability considerations are desirable, researchers have looked for ways to

i dA *A linear multistep method is called implicit if in (2.5) » By #0. If LY = 0 it is said to be explicit.

pe i v r

¢ Ww ~~

EY ie © ;

BN ; : (

pahlquist's concept of A-stability is seen to corre~ spond to A(n/2)-stability. Widlund proved that an explicit method cannot Lb i A(0)-stable and that for all ac[0,n/2) there exist Ala)~ i stable methods of third and fourth order. By allowing a to be less then n/2, methods as high g i as 12-th order have been reported [JAI70]). Indeed, Gear [GEA71a] and Brayton, Gustavson, and Hachtel [BRA72a] both have operational general purpose programs specifically de- signed for stiff systems. Both programs are similar and incorporate and make use of multistep methods up to sixth i order, namely the backward differentiation formulas [HEN62], each of which is A(a)-stable for some ac (0,7/2]. The effect of decreasing a is to restrict the eigen- values of the Jacobian matrix to a smaller region of the ’ complex domain wherein the method remains asymptotically stable for arbitrary, fixed h. ! | As the order of the integration process may exceed 2 for A(a)-stable methods, a<w/2, researchers have focussed their attention on obtaining higher order methods that are A(a)-stable, with a as large as possible, for the solution of stiff systems.

1 2.5 Runge-Kutta Methods (RKM) In addition to linear multistep methods Runge-Kutta ’

’ oo

v4 £ ~ . - .
1 15 type formulas have been used to numerically integrate (2.1). |] A g-stage REM may be written as [BUT64) ki Ypsr = Yn +B I Byky (2.192) i » where kyo ky» pp kq satisfy the equations

ky = fly, +b I Ajgkyr tq * czh). (2.19b) i=1 i

with i = 1, «..r 9» and )

c= 1 Ay for i = 1, «ees 9- (2.19¢)

Li All RKM require only one initial condition and hence are ) i self-starting. However, for every step taken, a g-stage Runge-Kutta process requires gq function evaluations* in B contrast to an LMM which only requires one. i If in (2.19b) Ay = 0 for j > i, the process is said to be explicit. Otherwise it is said to be jwplicit. The commonly used fourth order Runge-Kutta method is a 4-stage explicit process. It is given as EEE *Evaluating £(y(t)., t) is termed a function evaluation. If } this is a complicated function, as it oftentimes is, it may be a costly operation.

PS A |

Ls oF ——

y ~ ~~ 1 § 16 aid * + hk, + 2k, + 2k 4 + kde (2.20a) where Fa fly, Lp k, = £( + nk t+ i 2 Fy ¢ Phyo & + PV ¥ 1 1 (2.20b) ky = £ly, + Zhky, ty + FT) | li i | ky = £ly, + hky, t, +h). In the notation of (2.19) ’

1 2b, = by = by = 2by = 7 i oly = jel 22

Cy=Ry = 3 Ay,5= 0 j=2,3, 4, (2.21) 1 = oo " i Cy=R3,= 3 Ayy=0.3 1, 3, 4, [ Cq=Ag,3= ti Bg g=0 3=1. 2, 4 Though it has the desirable property of being explicit it is not A(a)-stable for any ac (0,n/2] and is hence unsuited for the solution of stiff systems. Butcher [BUT64] has explored the general class of g-stage RKM. He has shown that an implicit g-stage process can have order 2q. Ehle [EHL68] has shown that the root of ’ . the resulting characteristic equation is the Poa Padé ap- proximation to the exponential and hence is A-stable. { Implicit RKM, hence, can have arbitrarily high order

Ld ao

and still be A-stable at the cost of gq function evaluations per time step for an order 2q method. This is in contrast ‘ to implicit LMM which require only one function evalua, on { per time point.* Though implicit Runge-7utta methods are applicable to the solution of stiff systems of ordinary differential equa- | tions they were not considered in this research effort due to the high number of function evaluations required per time point. b 2.6 composite Linear Multistep Methods (CLMM) In 1971 the concept of a composite linear multistep method was introduced by Sloate and Bickart [SLO71b]. It I consists of ¢ linear multistep formulas used to simultan- I eously solve for & future time points.’ They may be written

¢ q *Both implicit LMM and implicit RKM may have to be solved by an iterative method resulting in a higher number of function

« 4 TN § 18

¥ I (03,5 Yngey = PBi5 Yngey) = 07 Ju-k+l [| with £ = 1, ..., t. (2.2 As with an LMM one may consider : : g 3 - h v ) Ly ly (tg ,) hl = I fog, g vty.) 8; y¥(tn, 40k jm-k+1 1 with i = 1, ..., 2. (2.23) i Under the assumption that y(t) in (2.23) possesses a con- vergent Taylor series at t = toe’ (2.23) may be rewritten

i, (3) Ly ly (tp) hl Icy, 407y 7 (Ey) . I 3=0 I with i =1, ..., %, (2.24) ‘ : where

1.3 _ 1 j=1, cy I Gree. wor Pia’ n=-k+1 with i = 1, ..., ¢ (2.25) The order of a CLMM may now be defined, as in Def- inition 2.4, which follows. pefinition 2.4: The i-th equation of the composite .

3 : 0 3 pe Send |

[ oF

linear multistep method (2.22) is said to be of order Py if Ego 00x]=9 Less P but Cipnn? 0 8 Cc is given by (2.25). let p= min {(p;)}. Ten 1.3 ledet © the composite linear multistep method (2.22) is said to be of order p. As with linear multistep methods the order of a g CLMM is related to its local discretization error as follows: If i) the composite linear multistep method is of order pi; : p+l = of k ii) Yne+s = y(t 44) + O(n )e 3 k+l, «ees 0;

. iii) y(t) is (p+l) times continuously differenti- able then y = gle.) + omP*l) for i=l 2 , ne#i neti p evoe Be The definitions of A-stability and Ala) -stability are also readily extended to CLMM with the words “com- $ posite linear multistep method” replacing the words #linear multistep method” in Definitions 2.2 and 2.3. Sloate and Bickart (SL071b] have offered the follow- ing conjecture: Conjecture: The order of an A-stable composite linear multistep method cannot exceed 2%. For ¢ = 1 Dahlquist [DAH63] has proven the conjecture. i Sloate [SLO7la] has derived a fourth order A-stable method ’

| |

¢ of Pan i 20 with Lt = 2, the limit anticipated by the conjecture. For only one backpoint, k = 1, the conjecture has heen proven [SLO71b), though the maximal order of 2¢ cannot be chieved H for ¢ > 1.* Though at first glance the introduction of a composite system seems to complicate the problem, as L equations must be solved simultaneously, the real benefit is derived from the hope that A-stable processes of at least order 2%, the limit hypothesized in the conjecture, can be found. i High order composite linear multistep methods have been reported. Reference has already been made to Sloate's algorithm [sLo71la) which is a 4-th order method and consists of two equations.’ A-stable methods of order £ + 1 which retain all future points and with only one backpoint® have ] been independently reported by Bickart, Burgess, and Sloate [(BIC71]) and Watts [WAT71). The former reference also des- cribes results obtained from a program utilizing these for- mulas for L = 1, ..., 6. Watts [WAT71] has proven that T— *For & = 1, the trapezoidal formula which has only one back- point is A-stable, of order 2. For k=1and £ > 1, only (2 + 1) degrees of freedom exist in choosing the parameters for each of the ! equations. Hence the maximal order . achievable is % + 1. tSloate's algorithm only maintains one of the two future time points solved for in each stage of the process. §The advantage of only one backpoint, k = 1, manifests itself in the ease with which the step-size may be changed.

gis,

« of ~— . p :

these methods are not A-stable for ¢ > 8 but has shown that 3 it is possible to achieve A-stability with order equal to h the number of equations, for all ¢, again with only oc back~ B point. ’ Composite linear mul .istep methods that achieve high { order, with the desirable properties of A-stability and only : one function evaluation per time point retained, reported so far have done so at the expense of increased computation 1 since f£ equations must be solved simultaneously at each stage of the process. As one of the primary objectives of this research effor:. was to derive processes that were ef- ficient, cyclic composite linear multistep methods were considered. (A composite linear multistep method is said } to be cyclic if in (2.22) 5,5 = Bij =0 for j >» i.) At I first it was hoped that high order A-stable cyclic methods ’ could be found. This goal was not achieved. However, the new integration methods described in Chapters IV and V are ' A(a)~stable and hence applicable to the solution of stiff systems. Moreover, they exhibit better stability proper- ties than the backward differentiation formulas in Gear's widely accepted algorithm. 1 2.7 The Characteristic Equation Consider the composite linear multistep method (2.22). The set of 1% scalar equations can be expressed ’

9 A!

¢ of ~—— i 22 in matrix notation as Y (a, Y, * Aate- * °° . A_x¥n-x! } -hiz, ™ B_1%n-1 +e +B lnk’ " 0 E I (2.26) where i [Fomsines Rens) Ya-3 = FN EO : § ¥ (n-3) 241 Y(n-3)2+1) : O1,-ja4L +o “Loa I Ag = : : ; (2.27) I By, -ja4s ove O2,-3241

By,-jass +++ P1,-jos1 B_j s : : | , By,-je4p +++ Br,-je+l for j= 0, .... K, K is the integer part of (k + & = 1)/%, and %.,m = Bi, m = 0 whenever m < - k + 1. Now consider using the CLMM as expressed by (2.26) I ’

: al p, \

P’ of § 23
to integrate the linear differentiaj equation y(t) = g y(t), or equivalently § .=q¥ .,3=0, cop) K. With ' = ah, ~n-j -n-J (2.26) becomes fi * BROKE, & ove ® (Ag = AB_)¥, x = 9 (2.28) which is a vector difference of equation of the form

I Ry (NY, 5 = 0, (2.29)

where Ry (1) H Ay - AB. Define

plg,A) = det { [ Ry 7} 2 (2.30) j=-K Then for a given step-size, ji.e., fixed A, (2.29) has a | solution, governed by the zeros gt, (2) of plL,A), given as '

{ , = 1 Cie (A), (2.31)

where the vectors Civ i=1, ..., KL are functions of the

k my

a. re ~~

initial conditions.* (Note: Of the KL zeros of p({,A).,

ER a . #The formulation as given by {2.31) can be proven by an extension of the scalar case as given in Section 2... As in the scalar case a simplification has been made in (2.31). In general one rust consider the invariant poly- nomials of the polynomial matrix [GAN60]

R(z,A) = Ry ¢*. (2.32) /

i | The invariant polynomials, 1, (z, 2), 0 1, (5.0) , are defined as - D (5,2) _ Cpl , i ye = CT with § = 1, ..., 2, (2.33)

where D4 (2,2) is the greatest common divisor of all g minors of order j of the determinant of R(gz,)) and

o-1 .n,, ' = 1k Ee ET, (2.31) i=1 j=1 o=1 where the g's are the s distinct zeros of (2.30) of multiplicity m, 5 in the j-th invariant polynomial

150.0 and the C,; 3 o's are functions of the initial ~dedr conditions. (Note: It is possible that m, =0 2 2 ro ne ri, 1. Moreover, m; i, < m 50) For the scalar case, i.e., '=1, (2.31') reduces to (2.11).

4 @h

‘ of —— / ] : [ ¢ TT \ + pi | KL - k of them are identically zero for all values of A.) H Equation (2.30) is called the characteristic polynomial of the composite linear multistep method (2.26). § Now if a method is to be A-stable, then Ay 4 as H nsw, Therefore by (2.31) a wothod will be A-stable if and F only if legal <1 for i=1, ...s Kt, for all yalues of g | A» with Re{)} < 0. similarly, if a method is to be Ala)- stable for some ac (0,m/2), then lz; | <1, for i = 1, | eevee Kb, for all values of A ¥ 0 such that arg(-i)| < @. | i In the next chapter the basis for necessary and suf- ficient conditions for graphically determining the stability | boundaries of a composite linear multistep method is derived, i.e., determine the regions of the A-plane for which all | [ solutions of a CLMM, when used to integrate a linear dif- i ferential equation, tend to zero as n * %. The character- » istic equation (2.30) is examined for those values of A I for which |g; (M] < 1.

’

3 - CT i

a. or ~——

i CHAPTER III H i STABILITY BY GRAPHICAL METHODS i i 3.1 Introduction H Consider the linear differenti»l equation : yi) = q y(t), (3.1) i and the composite linear multistep method > H

ag, Er A_g¥n-x’ a h{B YX, ame B_x¥n-x! a i (3.2) ¥ i Asymptotic stability for a CLMM will be as defined in

Definition 3.1 which follows. b] I pefinition 3.1: The composite linear multistep H I method (3.2) is said to be asymptotically stable

for Ay =q h, if all solutions of (3.2) tend to zero H [ as n + » when it is used to integrate the linear dif- | ferential equation (3.1).

If (3.1) is substituted into (3.2), it has been l shown in Section 2.7 that the solution is governed by the | zeros, gz; (0), of the characteristic polynomial of (3.2),

namely 0 | ll ple, A) =det { I (a; - TR B (3.3) i=-K

. “3 a iE a

‘ rg — / : (

' pefinition (3.1) is equivalent to saying the cLMM (3.2) is H asymptotically stable at the point A if all zeros gy ag) of S plz.) are well defined, that is, assume & unique v lue, and less than unity in modulus. § The characteristic polynomial can be written as the product of three polynomials as

pg.) = 2(2) P(E.) ADV, (3.4)

where z(g) and A(X) are, respectively, polynomials in § and ' A only, and pz. ) is a polynomial in both I and )\ that can- not be further factored into the form (3.4). If z(g) is a non-constant polynomial and its zeros, zy, are such that for some i EN > 1, then the method is H not asymptotically stable for any A. On the other hand, if for all i Iz < 1, then the stability of the method is not affected by z(g). For those values of A such that A(X) = 0 : the zeros Tg; of p(g,)) are not well defined (equivalently, i jll-defined) .* However, in a deleted neighborhood of each

l #Note: The zeros of A()) are also zeros of det {A -1B }; at the other zeros of det{a -2B }, the zeros gy are . unbounded (equivalently, undefined) . Thus, the zeros TM of ple, either are not well defined or are unde- fined at those values of A for which det{a -1B,} vanishes. Now, it is for just these values of )A that (3.2) does not J I have a unique solution when used to solve (3.1). ’

oe ; |

« of — / ; ( i 28 of these points, (FAL is well defined. (These zeros are ¥ discussed in greater detail in Section 3.2.) Considering the complex A-plane one may collec. ali : points A for which the CLMM is asymptotically stable into a set, S(1), i.e., s(A) = {A: CLMM is asymptotically stable}, (3.5a) 4 or, equivalently, sa) = {): 6; is well defined and i lg 0] <1, 4=1, ooo KAD. (3.5) : If S()) contains the open left half plane then the method I is A-stable; i.e., with L, = {A: Re{d} <0} if L,C s(}) ' then the CLMM (3.2) is A-stable. Note: In the complement of the set S(1), denoted as g(a), the composite linear ' multistep method is not asymptotically stable, by defini- [ tion of S(A), and, hence, there must be either at least one zero of the characteristic polynomial such that ty) is well defined and |g;(A)| > 1 or the t;(0)'s are not well

defined, corresponding to a zero of A(X). Let B()A) be defined as

( B(x) = 50 Ns) - (a: AN = 0}, (3.62) po

‘ a 1 a :

‘ 5s ~~

: where 513) is the closure of S(i). B()) may also be defined H - § B(A) = {A: for at least one index Jj ty) § is well defied and fey] = 1}. (3.6b) In the next section it is proven that B()) can be described by continuous curves, Yi (0). in the extended complex A-plane. These curves, vy (A), together separate the A-plane into ’ | regions* wherein the CLMM is either asymptotically stable or not asymptotically stable.’ i 3.2 The Lambda Plot The characteristic polynomial (3.3% may be written I as a polynomial in two variables as | 3 L

ple, = I I ade (3.7) §=0 i=0 ‘ J—— *The term components is used by odeh and Linigar [ODE71]. f tExcept possibly at a finite set of points where the (FAY are ill-defined, i.e., at those values of ) for which A(X) = 0. SThe characteristic polynomial has a factor Krk, i.e., 4 in (3.4) z(Z) has at least a factor ¢**k since this 1 factor in no way affects stability considerations, as it { only adds at most a multiple root at § = 0 for all values of A, it can be dropped from consideration. Therefore, in all that follows, p(Z,}) will be treated as though this ’ | factor has been removed. The resulting polynomial is hence i of degree k < K2 in L. |

= Bn .

p oy ~~ i 30 or alternately as

ple.) = I pytont, (3.8) so where

pyle) = I agyed, wien i=o0, oot (3.9)

Similarly, one can express pl(g,2) as ;

plz, 2) = [ aye, (3.10)

I where

I ayn) = a ts with 3 = 0, ..., k. (3.11)

It shall be assumed that (3.7) is irreducible.* If it is not, then p(Z,A) can be written as the product of irreducible polynomials and the discussion to follow anplies [ to each of the irreducible polynomials that is not a func- tion of A or r alone. Since (3.8) is a polynomial in ¢ and A, it may be EE . *A polynomial in two variable is irreducible if it cannot [ be expressed as the product of two polynomials neither of which is constant. By assuming that plt,)) is irreducible {A: A{()) = 0} will be empty, and need not be considered. These points will be reconsidered later.

1 une

p 4 ~~ / : (

written as* plz,2) = py (2) A-£, (2) [x-£, (2) Jeee [-£,(5)], (2.12) where each £, (0) is a wel. defined algebraic function of §, . analytic except at, at most, a finite number of isolated singularities. Denoting these singular points as c¢, they are either the zeros of the resultant of p(%,)) and ’ alp(g,A)1/3¢ or the zeros of p,(%). At the singular points Ch such that py (cp) ¥§0, Ahlfors [AHL66] has shown that each £,(0) is bounded, that is, these singular points are ordinary algebraic singular- . ities (equivalently, branch points). Furthermore, at those points c such that py (cy) = 0 there exists a least rational 1 number at non-negative and finite, such that £, (8) (gc) remains bounded, that is, these singular points are algebraic poles or, in the special case corresponding to m=0, are or- ! dinary algebraic singularities. Now, consider pted®,n) =o, (3.13) J— #The material in this and the following paragraph is excerpted ’ from Chapter 8 of Ahlfors, Complex Variables [AHL66]. The rational number m is bounded from above by the multi- plicity of the zero c, of P,(Z). .

: & : 3 8 hb

a. 5s — i 32 where j = /-I and 6 ¢ [0,2%). From the previous development, i (3.13) defines 1 loci NG

i 30 | / A; (8) - £; le ), for 6 ¢ (0,27) and isl, «cco 2, (3.14) in the extended complex A-plane. The basis now exists for pefinition 3.2: The set of points , i yO) = {az A =2,(8), 8c [0,2m), i=l, .... t}, (3.15) for a given composite linear multistep method, is termed the Lambda Locus, or the Lambda Plot, of the method. The Lambda Locus, Y(1), is a set of closed curves which divide the extended complex plane into a finite col- [ lection of open connected regions. Let N(}) be one such i 3 open connected region and 3N()) be its boundary, i.e., aN(A) = N(X) - N(A). (3.16) f Note: aN(A) € y()A) but N(A)Ny(2) = ¢, ¢ being the empty set. The stage has now been set for Theorem 3.3: For all points A contained within the open region N(1), one, but not both, of the following ’ statements is true: } ’ i) The composite linear multistep method is asymp- totically stable.

he k

a. 5g —~ / g (

ii) The composite linear multistep method is not asymptotically stable. NX The consequence of Theorem 3.3 is that yv()) s parates / the complex A-plane into regions wherein the method is either ! asymptotically stable or rc asymptotically stable. Proof of Theorem 3.3: Assume ‘there exists a point ); € N(}) such that A e SA), i.e., A is a stable point of the CLMM. ’ | Next, assume that there exists another point A, e N()) such that Ay ¢ S()). Now, consider the two equations p(t.) =0 4 and p(z,2;) = 0. oy a development similar to (3.12), we have,* with i = 1, 3 L plz,2y) - qy (3) [z=g, (Ay) ] [z-g, (24) Jree [Z-gy (3;)] = 0. (3.17) { Two cases need to be considered as follows: ! Case 1. Neither A nor 1, is a singular point of g; (0), §=1, ..., k. Note: under these | conditions a, (33) # 0, i=1,2. case 2. Either, or both, Ay or Ay is a singular point of 950) for some J. ’ Case 1: By construction, for at least one index Jj, lgg Op) | <1 .

*See (3.10) for a definition of gq (A). '

_ or ~~ BE g and lg5x01 > 1. If EAP = 1 then 1, ¢ v(}) but NO OAYO) = 6. \ Hence EAP >1. As 950) is analytic except at a finite number of points there must exist a path from LY to 1,. f say L, completely contained within N()) and along which { 950) is analytic. As CAEP > 1> |gs(,)| there must exist J & - a point A; € L such that EATS =1, i.e., A; € (2). But N(A) )Y(2) = ¢ which leads to a contradiction that i; € ¢. ; Hence, A and Ay cannot both be regular points of N()) with Ay € S()) and Ay g£ sh). Case 2: As before, ror at least one index ji, EAT < 1 and lg; 021 > 1. Again it can be assumed that EAP # 1 for : otherwise A, € vy(1), which contradicts Y(AODNN) =o. i Assume Ay is a singular point of 950) but ay (2) #0, ] i.e., ), is an ordinary algebraic singularity of 940). ’ § Now, consider the function 650) defined as Gyn = lay 1s (3.17a) where ( lim Gy (2) = xd, G5). (3.17b) The function G5 (0) is well defined and continuous in a neighborhood of Aye Hence, for each ¢ > 0 there exists a f v

p Eg . 35 é > 0 such that [65023 - G51) <e€ (3.18) for all i, for which |i; = Al <6. Therefore, as Gy (25) > 1, ) an ¢ can be chosen such that G5 (23) > 1 for A, € N(A). At ) A= dg 9; (0) is analytic and |g503) | > 1. Therefore, by a similar argument to that in Case 1, a contradiction results. Next, suppose a, (2,) =0, i.e., A, is an algebraic pole 1 of 950) .* thus g95(2,) + ® as A + Aye In some sufficiently small deleted neighborhood of i, 952) is analytic and lag] > 3. ] Again a contradiction results. S Now, let A be a singular point of 950). Note: A cannot be an algebraic pole, by construction; hence gy (34) #0 ] and A is an ordinary algebraic singularity, As before, ' consider the function G5 (N) defined as

  1. = lay 1. (3.19a) where JRU—— ] *The zeros of 9 (N) correspond to those values of A for which the zeros of the characteristic polynomial, ZC, (a). | are undefined, i.e., gg) = det {A_-)B J. Further, by assumption, none of these points are ill-defined—A(}) = i : constant. ’

~ §

“__ 4 ~~ 36 { 4 lim 2 G50) = A*dy Gy). (3.19b) G50) is well defined and continuous in a neighboxhc 4d cf / “ / Ape Hence, for each ¢ > 0 there exists a § > 0 such that 16523) - G50) <e (3.20) for all A, for which |i; - Al < 6. Therefore, as Gy(2y) <1, an ¢ can be chosen such that 6523) <1 for Ay € N()). At [ | A= Ay 95) is analytic and lay ag) | < 1. Once again a con- tradiction results. Q.E.D. } From this theorem follows 1] Corollary 3.4: If N()) is an open connected region in the extended complex A-plane determined by the Lambda Locus and if there exists a point A € N(}) such { that the composite linear multistep method is asymp- , totically stable at A then the CLMM is stable for all A € N(A).

Corollary 3.5: If in any open connected region, H(A), of the extended complex A-plane determined by . the Lambda Locus, there exists a point Ao such that n<k of the values pg.) = 0 are greater than one in ) modulus then this condition is true for all A € N(}A). ’

Moreover, if N,N) and N, (0) are two such regions

Pp Eg —— Ve E (

ky) such that an, (NaN, (0) # ¢ then as one crosses from B N, (0) to N, (2) at least one root of the character- istic equation has gone from being less than unity i in modulus to greater than unity, or vice versa. I As previously noted, y()) divides the extended com- > plex A-plane into a collectiqn of open connected regions. d At most +l such regions can be formed, each of which must 1 be tested to determine if the composite linear multistep method is asymptotically stable or not in that region. ] Moreover, as the set of points in the extended A-plane de- fined by (3.13) is a mapping of a closed contour, its image I in the extended A-plane must also be a set of closed con- ' I tours [AHL66]. Hence, if one of the f,(e3®) of (3.14) does not form a closed contour in the extended A-plane as 6 goes . from 0 to 27, it must join at least one other segment to ( complete the contour. ‘ It has been assumed that the characteristic polynom- jal is irreducible. If it is not then the results of this section apply to each of its irreducible factors, say py (2.2), Eel pn (E/N). Associated with each irreducible factor will be a stability region, s;(). The stability region of the CLMM, s()), is then given as SA) = 8; (A N+=- NS (A). (3.21) .

’ ’ 5 1 CI 5

38 3 Excluded from the discussion have been polynomials which are functions of \ alone and of ¢ alone; i.e., it has been assumed that in (3.4) both A()) = constant and Z(g) * constant. If A(A) is a non-constant polynomial the zeros (FALY) of the characteristic polyuomial p(g,A) are ill-defined . J at the zeros of A(A) and {A: A(X) = 0) cannot belong to the ) stability region S(}) of the method. If Z(g) is a non-con- stant polynomial and its zeros, z;, are such that for some i, lz] > 1, then the method is not asymptotically stable for l any A. On the other hand, if for all i, lz <1, then the stability of the method is not effected by Z(%). An example of a factorization of this type is given . i in (3.22) [WAT71]. I Yne1 = Yn » ny, ii neo! (3.22) Yn+2 =¥, * iy, ud Yne1! The characteristic polynomial is given as Ey [4 -1-a(z+1)) plz.) aoe{( 8 z-1-1 )}- (3.23a)

p(g,A) = ¢l(z=1) = (Z+1)A1(1%A). (3.23b) In the notation of (3.4), '

Ek hk E

_ E38 ~~ / : : (

z(g) = ¢» | pg.) = (g=1) = (g+1)X, (3.24) 5

i AD) = 14. | The zero of Z(g) at ¢ = 0 douse not affect stability and can be ignored. The stability region for the irreducible polynomial i TRY) consists of the open left-half plane as can be seen | from

I Plea) = - (54D) O-55P - (3.25) I i The Lambda Locus for pg.) consists of the imaginary axis " | which closes at infinity. Moreover, plz.) has an algebraic i pole at A = +41. The zero of A(}) at A = -1 indicates the zero | of plz.) is ill-defined there; hence, the point A = -1 must ' be excluded from the stability region determined by a con- I I sideration of plz.) alone. The Lambda Locus for the irre- ! | ducible polynomial A(X) consists of a single point, * namely red. 3.3 Stiff Stability For a composite linear multistep method to be useful | for the solution of stiff systems of ordinary differential Sie | *Note: The Lambda Locus corresponding to A(}) consists of ’ 1 isolated points ccrresponding to the zeros of A(}). | A! Cv ; 3 A —\

_ “ —— Ve : . (

i equations one desires asymptotic stability in a "large" i region of the left half of the complex A-plane. The concept - of stiff stability includes this attribute. Stiff + tability J may be defined as in Definition 3.6 [BIC72], which follows. / | pefinition 3.6: (et Ss, be an open connected region of the extended complex A-plane such that i i) {r: Re{)} <v = o}cs,. ii) A=0c¢ 5, Then the composite linear multistep method (3.2) is i said to be stiffly stable if all solutions of (3.2) tend to zero as n+= when the method is used to inte- i grate the differential equation i y(t) = q y(t), with any fixed positive h, such that I gh = 1¢€ Sy- I If a ‘method is stiffly stable and if {A: A is real and negative}C Ss, then it is A(a)~-stable for some a € (o,n/2]. Moreover, a method is A-stable if it is stiffly stable for all A € Sgr i.e., the method is stiffly stable with vy = 0 in pefinition 3.6. a The development of Section 3.2 permits one to deter- : mine the stability boundaries of a CcLMM and hence determine if it is stiffly stable or A-stable. The method may be sum- marized as

Cw ; 8 —

p F's —— — \ <\

Procedure 3.7: Given a composite linear multistep method and its characteristic polynomial p(C,)2), the h mapping of the unit circle in the g-plane into the J A-plane defined by (3.13) is applied, dividing the extended complex A-plane into disjoint regions. (At most +1 such regions exist.) If a region wherein ‘ the CLMM could be stiffly stable exists, an interior point is chosen, say Le The k zeros gy (0) of Plz,2)) are determined and, if lz; 1 < 1 for ] i=1, ..., k, then the method is stiffly stable. If, in addition, Ly is contained within this region, the i method is A-stable. ] In applying the mapping (3.13), one must be careful 5 to exclude the points in the A-plane for which the zeros . (FALY of plz,)) are ill-defined from the region Sy" Procedure 3.7 has been implemented in a computer pro-= t gram for determining the Lambda Locus of a method. Note: Definition 3.6 requires A = 0 to be in 5, and implies A = = is in 5, Consequently as an aid in determining if a CLMM is stiffly stable the program solves for the zeros of py) and Po (2), corresponding to {g;(=)} and {z; (01, respectively. . —— *These points are the zeros of A()). Each of these points, if any, must be determined and excluded from S,-

p Fe ~ 42 For A = 0, p(Z,)) must have a zero at { = 1;* thus the point A=0f£ 5, he legal, i=l, ..., k, is a continuous function N i of A, if the method is stiffly stable {; (0) and I;(=) can- / not have modulus greater than unity.’ If for some i this p condition is not met then the method cannot be stiffly stable. Moreover, if lgg=)] < 1 and a region does exist

#If the CLMM is to be at least of order zero, then each row ’ i of the matrix

i a= J] A i=-K must sum to zero; hence,

det { I ac) = p(z,0) = p (2) 1x y must have a zero at { = 1. . tan even stronger statement can be made as follows: The modulus of all roots of py (2) = 0 and Po (%) = 0 must be bounded by 1 and the roots of modulus 1 must be "simple" [HEN62]. ! If py (2) or py (8), say Po (2), has a zero of modulus 1 with multiplicity greater than 1 the minimal polynomial of R(Z,A), X(%,)), must be examined at A= 0, If x(g,0) = Po (2) the method cannot be stiffly stable. Howevel, if x(z,0) # [ Po (2) the invariant polynomials AL §=1, vees 2, a8 l defined by (2.33), must be examined at A = 0. If for all i the zeros of 1502.0) have multiplicity 1 then the method [ may be stiffly stable. Otherwise, the method cannot be : stiffly stable. Note: x(Z,) = 1, (5.0) and Lig (BoM 25080 for j=1, ..., %-1.

l x ~~

such that the method might be stiffly stable (that is, a i connected region exists, say S, such that A = 0 ¢ § and A = ¢ S), then the method is stiffly stable and t. ? Zeros A Hi FAY of pl(¢,)) for an additional interior point of the region - i need not be checked. This follows directly from Corollary 3.4 and the fact that in this instance A = « is an interior point. 3.4 The Zeta Locus Since the characteristic polynomial of a composite i linear multistep method is a polynomial* in ¢ and A, it may be written as plz.) = gy (W) (5-9, (115-95 (M1 += + [5-5 (M11 + (3.26) I where each g; (MN) is a well defined algebraic function of 1A, analytic except at, at most, a finite number of isolated

singularities. Denoting these singular points as ane they are either the zeros of the resultant of p(%,)) and 2[p(z,A) 1/3) or the zeros of q (A). Let C, "xe the contour in the A-plane shown in Figure 3.1. As the point R on the ————— #*As was in the discussion with the Lambda Locus, it is assumed that the characteristic polynomial is irreducible. If it is not, then p(%,) can be written as the product of irreducible polynomials and the discussion applies to each of the irre- : ducible polynomials that is not a function of A or § alone. The special case wherein one of the irreducible factors may ] be a function of A or i alone is considered later.

a F's yy

i 44 [] —

7 AT JO —— /| \V/ A) / 3 A i | / JriEmann SPHERE iy 4 —— TN] 7’ COMPLEX I c A-PLANE 1

I PROJECTION OF Cy ON : I THE RIEMANN SPHERE Figure 3.1: pefinition of the Contour Cy

k Cr

_ 3 ~~ i 45 Riemann sphere of Figure 3.1 moves towards the point T, <5 i completely encloses the open left half of the complex A- plane. i Now, consider I plg,A) = 0, 2 € Cy. (3.27) I As in Section 3.2, (3.27) defines k loci, | 5; (0) = g; (2) with i=1, ..., k and A € Cy (3.28) in the complex A-plane. The basis now exists for Definition 3.8: The set of points reg) = {gx ¢=3%;, 0), A eC, i=l, ..., k}, (3.29) : I for a given composite linear multistep method, is

termed the Zeta Locus, or the Zeta Plot, of the ' I method. The Zeta Locus is a useful tool for checking if a method is A-stable. Theorem 3.9: Let U, = {g: |g] <1}. If the singular i points of the characteristic polynomial of a composite linear multistep method, denoted as dpe are such that a, F4 I, then the composite linear multistep method is ° A-stable if and only if its Zeta Locus is contained within T,.

: Be

_ rs —

46 Proof: The proof is a direct consequence of the Maximum Modulus Theorem of complex variables [AHL66]. Under the assumptions of the theorem the g; 's are analytic for a i Ae I: hence, each assumes its maximum on the boundary. If for € C, |ggtA)| <1 then |g (0 <1 for all A ¢ L, and by ) definition the method is A-stable. On the other hand, if a method is A-stable, then for i ¢ Lys FAY < 1. Therefore, I leg] <1 for A € Cy, and the Zeta Locus must be contained within o.. Q.E.D. ; The Zeta Locus, -y itself, does not give necessary and sufficient conditions for checking if a method is A- ) stable, since the location of the singular points must be i determined and checked independently. Now, let an denote those singular points that are in the left half of the com- : plex A-plane, a situation in contradiction with the primary I hypothesis of the previous theorem, and let Ly = L, = az}. | ! As each dan is an isolated point, the g; (M) are analytic [ everywhere for A € 3 and, therefore, each assumes its maxi- mum on the boundary. A proof similar to that for the last theorem easily establishes Procedure 3.10: Given a composite linear multistep method and its characteristic polynomial p(¢,A), the mapping of C, in the A-plane into the g-plane defined . by (3.27) is applied, thereby specifying the Zeta Locus T(z) for the method. If rev, then the |

= 2 CT 2 2 NX

singular points of pl(g,)) in the A-plane, a. are determined. The k zeros gga) of plz.ay) for each a, € L, are determined and if all are contained i. U, the method is A-stable. s 4_ such that gq, (d = } - At those singular poinus 4d. 8 I mn’ 0 there exists at least one index i such that g, (A) contains an algebraic pole at i = doe i.e., g; becomes unbourded as A+d_.* Hence

| Theorem 3.11: A necessary condition for a composite linear multistep method to be A-stable is that the zeros of qa), as defined by (3.11), be in the open right half of the complex A-plane. Note: Singular points may exist in L, and the method could 2 still be A-stable. However, they can only be branch points, that is, they cannot be algebraic poles, namely, zeros of : gy (N). Moreover, if the singular points are not checked,’ erroneous information may be inferred from the Zeta Locus.’ J — *Though g;(N) becomes unbounded, there exists a positive integer n such that (=a) "g; remains bounded as +d. tpetermining the singular points that are not zeros of gq, (V) is a non-trivial task in practice, though possible. Sas an example, the Zeta locus, alone, for the backward difference composite linear multistep method proposed by G Bickart and Picel [BIC72] would indicate that the 6-th ’ order method is A-stable. The Lambda Locus does, however, indicate that it is only stiffly stable.

/ : (

The Zeta Plot, as defined by (3.28), cannot yield any information as to whether or not a method is stiffly stabie | A if it is not A-stable. However, a slight modification of the above formulation establishes | Corollary 3.12: Let uv, = {g: gl] < 1) and let Ly be an ) | arbitrary open connected region of the extended com- plex i-plane. 1f, for a given composite linear multi- | step method: } i) All singular points of its characteristic equation, denoted as ane are such that for | I a ify a, c Li. gg ta) € Up RE se k. =i 1 ii) The k loci, FAY = g; (N) with i=1, ..., k I and ) € aL} are contained within U,. i Z Then the composite linear multistep method is asymp- . totically stable for all A ¢ Li: In particular, if Ly i [ is chosen such that Spr of Definition 3.6, is con- ' | tained within Ly and conditions i) and ii) are met then the method is stiffly stable. | In practice the Lambda Locus has been found more use- | ful in determining if a method is stiffly stable, since Corollary 3.12 requires an a priori knowledge of the set ( Sys which is unknown at the outset. It has been assumed that the characteristic polynomial [ [ is irreducible. If it is not then the results of this section ' [ apply to each of its jrreducible factors. Each irreducible

Gi 3

p ES ~~ / i? (

49 factor and its associated Zeta Locus must be investigated as i discussed in this section. Excluded from the discussion have been polynomia’ 3 ! which are functions of A alone and ¢{ alone; i.e., it has been assumed that in (3.4) buch AQ) = constant and 2(g) = constant. If A(X) is a non-constant polynomial and if 3 {xz AN) = 0}CL, then the Zeta Locus will yield erroneous results. Note: At those points {a: AQ 0} the zeros 1 zg, (0) of the characteristic polynomial p(g,)) are ill-de- i fined. An example of a method exhibiting this type of SEN factorization was given in (3.22). As the Zeta Locus for I that method consists of the isolated point { = 0 and the unit circle in the g-plane one is tempted to assume that the : i method (3.22) is A-stable. However, as A()) has a zero at i A=-1¢lL, the method cannot be A-stable. ' If z(g) is a non-constant polynomial then the Zeta I Locus corresponding to the z(g) consists of isolated points, ! namely, the zeros, 2;, of z(gz). As long as Iz] # 1 for all i the Zeta Locus may be interpreted correctly. However, if for some i, 1241 = 1, a casual examination of the Zeta Locus may lead to erroneous conclusions. An example of this type of factorization is given in (3.30) .*

[ #This example was offered by a colleague, W. Rubin of the Department of Electrical and Computer Engineering at [ Syracuse University. ’

CT. :

‘_ ox —— — .. 4

-12y,,+24y,-12y,,, = h{Sy,;* 3, "%%n:1*¥ne2 } (3.30) By," 127 Hn, = h{(=3¥, 1-59, * 541% 042) The characteristic polynorial is -12(g+1) =A (=9¢+5) 24-2 (z+) plg,2) = det |}, (3.312) 8-a(7g-3) 4(g-3)-2(g-5),) ’ or plea) = ~(z-1) [(3-3an?) - (+302). (3.31b) ] In the notation of (3.4), \ z(g) = -(z-1), [ el a 2 2 plz, 2) = (3-332) - (3+3+27)C, (3.32) and A(X) = 1. "The Zeta Locus associated with plz,) for A c, consists of the unit circle in the g-plane while the Zeta Locus corres- ponding to zZ(z) consists of a single point, namely ¢ = 1. As the Zeta Locus for the method (3.30) is the unit circle and the isolated point { = 1 one is tempted to assume that it is an A-stable method.* However, as (3.31) has a zero for all

JE ——— *Purthermore, the singularities of the method are in the ' right half of the complex A-plane. |

‘ « ~ i 51 A for ¢ = 1 the method cannot be A-stable. t 3.5 A Plotting Program One of the primary orjectives of the work reported ¥ herein was to obtain new and precise methods by which to characterize the stability properties of composite linear i multistep methods. To that end the relationship of the i Lambda and Zeta Loci to those stability properties has been established. To use those relationships in the search for [ | "properly" stable CLMM's a computer program was written® for plotting both the Lambda Locus and the Zeta Locus of a i CLMM. A block diagram showing the five segments cf the pro- i gram is given in Figure 3.2. After program execution is begun control is passed I immediately to the subroutine ACCEPT, which returns the characteristic polynomial in a doubly dimensioned array [34 real numbers, 1cuar. 5 In particular, the characteristic Jo ais *The program was written for execution on 2a XDS SIGMA-5 ' computer. The program is in non-standard FORTRAN making full use of the power of XDS EXTENDED FORTRAN-IV. As an example, many arrays are indexed from zero rather than one as in standard FORTRAN-IV. Conversion to standard FORTRAN is readily made. Listings of all programs cited can be . found in Appendix A. fail computation is in double precision. On the XDS SIGMA-5 computer this corresponds to 56 binary bits of mantissa, or approximately 16.8 decimal digits of accuracy. { Sa11 capitalized words correspond to actual program variable { names.

! SE Ek : ¥ 4

: 52

Es5pg%° .

® < i 3 = 280 28; I] $2. ; $2 4 = = $3: gli: ji: - & Ba 8.3 i LEP dgiz :

Exif 1H g 148 zt a iad o

» : |

£1 8 2x .

I 3 ob & EEE J “ea on

tH : bi : : -~

i. of $2885 1883 2883

cand 3 §iic3 HA FE EME 58%; st gx® ' FERRE

“ g

$ hi.

— \ . \ I

53 polynomial is NL NZ ple.) = § I zcaarcr,nchl. (3.33) 1=0 J=0 For a CLMM consisting of & eguu.ions and k back points,* NL = ¢ and NZ = Kk. i Though only one program is shown in Figure 3.2, there are in reality many programs. This is a consequence of the l fact that numerous ACCEPT subroutines have been written. i i One version reads a CLMM from cards and computes the deter- ' minant as given by (3.3). Another version reads the charac- ] teristic polynomial directly from cards. Each returns the characteristic polynomial, for the moment assumed to be of I maximum degree NL in A and of maximum degree WZ in L. In addition an eight character alpha-numeric identification is ' supplied. This identification appears on all plots and printed output. The type of locus to be plotted, i.e., ' Lambda Locus, Zeta Locus or both, is also specified. Control is then transferred to subroutine INZUT which prints the characteristic polynomial on a line printer and determines if the method is potentially stable by examining the zeros of the characteristic polynomial at A = 0 and A = =, . —_— : This assumes that the factor ¢¥"K nas been removed from the characteristic polynomial and that all future points are re- . tained. If all future points are not retained then NL < L and NZ < k. The program can accommodate this variation of a CLMM.

2 E EB

— A H 54 With the characteristic polynomial written 2s in (3.8), i these zeros are, respectively the zeros of Po (8) and pyle) .* These zeros are computed and printed on the printer. It is possible for a characteristic polynomial of a CLMM consisting of L equations to be of degree n < L in A. An ] example of this phenomenon is found in the CLMM consisting of the forward Euler formula and the backward Euier formula I applied cyclically: [ Yns1 = ¥n + hy, (3.34) I Ynez = Yne1 + Pns2 [ The characteristic polynomial is given by [ t -+1))) pl(%,}) = det ’ (3.35a) i -t t(1-2) ' . or j plg,A) = Zl(=1+2) ~- (142) A). (3.35b) i Note: p(%,) is linear in A even though the method consists ——————— h *as discussed in Section 3.3 these zeros are useful in deter- mining if a method is stiffly stable. p

p: ——

« 3 ~~ 55 of two equations.* Subroutine INPUT checks for the possibility that n < 1 in the process of determining the zeros at h = =. If the - on~- dition does exist, it decrements NL oy one and uses Py-y (%) [ to determine stability at h = =», The procass 1s repeated as is necessary. { Control then returns to MAIN where the points com- prising the loci, as defined by (3.13) for a Lambda Plot or by (3.27) for a Zeta Plot, are computed. y The MAIN program segment consists of three sections: initialization, start-up, and loci computation sections. The initialization section, shown in flowchart form i in Figure 3.3.7 starts by reformatting the characteristic J SEER *In fact the method (3.34) is A-stable as can be seen by setting p(Z,A) = 0 and solving for ¢{. Thus 1+) !

: (The constant zero at { = 0 in no way affects stability and —- be ignored.) It is readily verified that lz| < 1 for Rre{g} < 0; in fact the characteristic equation is the same as that for the second order trapezoidal rule which is well f known as being A-stable. However, this method is just of order 1. But of more interest is the fact that one of the equations is explicit. Hence, an implicit equation need be solved only every other time. The possibility of inter- mixing implicit and explicit equations to form a CLMM is an 3 interesting concept worthy of further research. Tan statement numbers refer to actual statement numbers in the program, a listing of which can be found in Appendix A. r Ad

3 es.

#* ( » =

56 3] : 53) 3 —%-i ® $1: | [4 I | rE | g

, | e t\d IT = pel ‘ A EA 2

CE FFE a 1 J '

sie 45% 3: [382 AVE - : HER Hi FF ff 8 Bd

: y BE CR

4 3 a [e ' in -

5 id 3 | . -

By i . 1 [1] we Pp i - 5 38% " | LEE Se | i §2:i3 (it " { gS EF [22% AJ

g 5,

— \ x T i 57 polynomial from the array ICHAR to a second doubly dimen- i sioned array CHAR. The array CHAR defines a polynomial [34

degree NROOTS whose j-th row* is the coefficient of ¢ [37] ] for a Zeta [Lambda] Plot. Thus the entries of each row determine a polynomial of acgree NDEGREE in A [g] for a Zeta [Lambda] Plot. NROOTS and NDEGREE are set to NL and NZ, © | respectively, for a Lambda Plot and to NZ and NL, respective-

2 ly, for a Zeta Plot. The parameter VAR, a complex variable, 1] l is initialized. VAR assumes the values e3® for a Lambda

§ ] Plot and jw (A = 0 + ju with ¢ = 0) for a Zeta Plot.

F puring the start-up procedure, shown in flow chart

p | | form in Figure 3.4, the starting values for each of the | CER 4 3 NROOTS loci are computed. 5 i : EL =] For a Lambda Plot 6 is initially set to -m. AS only ! g | half the locus need be computed’ ¢ is varied from -w to 0. ’ |

o ———————————

a i *Indexing starts at zero.

g The locus for © varying from 0-7 is the complex conjugate of the locus corresponding to ® varying from 0+-7 and hence

‘ 3 ) — / : [ B ® A 58 i Tan a pom ~ i 1 Epa) TX iy TT aos oF Tae Pesvaenit . Sen sasmsaves sen ¥ Eo A I Petes mu Rew TN TE are ® Aue Reeth Pewed Manes meer sem ReeviD) 2 TE Tuy Root (AReOTS), | 4 NANT MER eerhe eT ree" 1 oF Pu vusmtas Ts i | | Co a— » fs [fin » ram TG RS - RE 4 Swney | a a i Ra fr. 0 ms | I Te ror? l Na {= (r= 3 = LS 4 3 2D ne “ *) )

WE x I— 5 ele SL 1 Bs | - | 3 .| 1

i a: | ro = 3 | Pp

Jo——— i ) [ S = ee Zo Y E ey. ©) Voie Genes Lee—y 4 = Gar> : | Seane | A C5 © ~~ har a Nv Tor ©

0" a ds 3 ar wary Suny, SaOnn y Lei Figure 3.4: Main Program——Start-up Procedure

EL v

At each new plot point 6 is incremented by w»/INCR.* For a Zeta Plot w is varied from 0 to Ar’ (chosen as 100%) in a ( geometric manner as follows: The start-up value is alv ws ! chosen as w = 0. After the initial point for each locus is | computed w is set to AI (chosen as 0.1). At each new plot point, ©oe1 is chosen as EE ———— As ay" 35 —~— 3 L
4 = 2 iJ ord p(t. 2) | II agp ‘l = p(z, 2). §=0 i=0 Hence, if for a given value of {, say Sor Ao is a root of the characteristic equation, i.e., PZ s2;) = 0, then p(T 2p) = Be This implies that LN is a root of plz.) = 0. - As "3° = 3%) the claim is proven. *INCR is chosen as follows: A singly dimensioned complex . array ROOTS (dimensioned at 1500) is used to store all points to be plotted. Because jt is desirable to compute a maximum number of plot points, INCR is chosen as the largest integer satisfying (INCR+1l) NROOTS < 1500. Each ‘ f locus, i.e., each £00) or 95 (0) of (3.12) or (3.26) for a Lambda or Zeta Plot, respectively, is evaluated at INCR+1 1 values resulting in a plot comprising (2xINCR+1) NROOTS points. ‘ Tin an argument similar to that for a Lambda Plot only half the locus need be computed as the locus for w varying from 0+-= is the complex conjugate of the locus corresponding to w varying from 0+=.

{ Ste has been found that, in general, a value of w= 100 is close enough to ® so as not to adversely effect the appear- { ance of the plot.

o & V 3 bs

—

“a4 " DELTA>® , (3.36a) where DELTA = (AP/AT) (M/INCR), (3. %b)

Varying wv as indicated by 17.36), rather than linearly, has been found to yield better results in the smoothness of the 1 plots. 3 For the complex value VAR, 2a polynomial is computed ’ H which is a function of only one variable—i if it is a Lambda 1 Plot; ¢ if it is a Zeta plot. This polynomial, of degree NROOTS, is formed in the complex array POLY, given as i NDEGREE p pory(3) = [I [CHAR(I,)] [VAR], a 1-0 i with J=0, ..., NROOTS. (3.37) J ] This procedure is repeated at each plot point—value of VAR. J The zeros of this polynomial, whose coefficients are complex numbers, are initially computed by the rootfinder COMROOTS using a Newton-Rhapson jteration to extract one root [ at a time.* If the rootfinder indicates success each root liad . *Al1 zeros are computed with equal accuracy as follows: ( Assume that a root of the polynomial p(x) has been found, { say X,;- A Newton-Rhapson iteration is then performed on . j the reduced polynomial p(x) /(x-%x,) to locate a second zero. . |

£ CC ® :

1 \ ( {

This value is then used ..= an initial guess for a Newton- } Rhapson iteration on the initial polynomial, p(x). result- . ing in a more accurate value, say X,. The process is re- peated for the third zero, this time starting with the re- duced polynomial (p(x) / (x-%,)1/ (x-%;) to obtain an initial guess for a Newton-Rhapson iteration on the original poly- , nomial. This process is repeated until all zeros have been extracted. By returning to the original polynomial at each stage, the error incurred in determining one zero (due to round-off, etc.) does not effect the accuracy to which a second zero is determined. Logic to detect multiple zeros is also built int. COMROOTS. Note: 1f two zeros of the L polynomial, say xy and Xyr are close to each other, but not equal, the procedure outlined may assume that a multiple > zero exists. Such a condition can occur if the intial guess -

for the Newton-Rhapson iteration on the original polynomial for the zerc Xx, is "closer" to Xx, than to x,. In this case, the procedure will erroneously yield x, as the zero. This, however, does not present any major problem. Multiple zeros are checked for in the locus computation segment and handled as a special case. This procedure is described later. *Convergence is determined in the following manner: The (n+l)-st guess at the root is given as *n41 = *n i 6x0 where x is the best guess at the end of the n-th stage and ox, = -p(x,)/p' (x), p(x) being the polynomial whose zeros are desired. The process is terminated when [16x] [/g(x 4) < EPS (set at 10710), where : x11. [xl] >1 { g(x) = y , 1, xl <n and v 11x]] = |Reix}| + |zm{x}].

A!
¢ ®

(¥]

Success at this stage—namely, all roots found-—results in being transferred to the last section, locus computation, of the program. If COMROOTS failed to converge at the starting value of VAR and the highest order coeffic ent is non-zero (cor- responding to IER = 3 in the flow chart) an attempt is made B to rescale the characteristic equation. All coefficients of the array CHAR are divided by the modulus of the largest i coefficient in CHAR. The polynomial (3.37) is then recomputed B and the process repeated. If it fails a second time, the program aborts. [|] If COMROOTS failed to converge and the lighest order coefficient, POLY (NROOTS), was zero (corresponding to IER = 4 I in the flow chart), it is assumed that the point is an alge- I braic pole—a zero of py (2) or gy (N) for a Lambda Plot or a 3 zeta Plot, respectively. In this case, a message is printed and the point is omitted from the plot. The next plot point is chosen as the starting point and the process repeated. If the same condition occurs a second time* the program aborts.

{ *As the poles are isolated singularities the condition will re-occur only due to numerical noise in the evaluation of E POLY (NROOTS) or the unlikely possibility that a second pole bas ng found at the particular spacing of the plot point pacing.

1 -~

63 It has been found beneficial to know the extent of i each branch of the locus.* Hence, the starting values, as well as the mid-point values’ are printed. In the last part of the FaIN segment, shown in flow B chart form in Figure 3.5, tu: variable VAR is set to its next value, as previously discussed, and the polynomial whose coefficients are given by (3.37) is re-computed. For each branch of the locus, the zeros of the polynomial at the previous setting of VAR are used as an initial guess to the zeros of the polynomial at the current setting of VAR. How=- ever, since it is possible that VAR may assume the value of a branch point,’ a check is made to sec if two (or more) of the previous zeros coincided. 1f they did, the zeros of

the polynomial at the new setting of VAR are all recomputed | I using COMROOTS.** In addition, in an array ILOCUS, the ’ J —— *This also aids in determining if A()) or 2Z(Z) in (3.4) are ' non-constant polynomials. tas only half of the locus is computed, the last point com- puted is the mid-point of each branch. { Sif an algebraic pole is encountered in this phase of the plotting, the program aborts. ( fThis condition can also occur if two (or more) irreducible . { polynomials assume the same value for a given setting of : VAR. The program handles this situation in a like manner. ##1f this procedure were not followed then REFCROOT would fail : to distinguish one branch from another and one locus, at least, would disappear from the plot. Hence all zeros are recomputed anew. ’

i. v

KL. HE _

t : . ig \ Lal Gli o off

ROH] [3 Pee Hd ol AURA H LD E §| 5 Ni | fl a : 3 | do fib 3 1! 19 | [keys | : 28% 4 h | | : if AN Naly, PRE Bi Dif +O ii Pa ima Va i AR we : SEE ; Ee | : | A \ i 0 : to iq LE— 5 = ih = Mel — p | yet | OB mn 4 5 pp | ; ; itl bd : | iL : Pi ,358Y -

= / $35] | B

ih ) 8 I fl Bt © ii hi 31) | 135] i" r=

— \ \

65 current plot point is stored, so that the plotting routine will have the information needed to lift the pen at this point. The points comprising zach locus are stored in the i complex array ROOTS. Th: .J-th locus (J=1, ..., NROOTS) is determined by the points stored at ROOTS (J+ (1COM-1) NRCOTS) , ! B for ICOM=1l, ..., INCR+l. — computed the subroutine PLOT- | TER plots these points. For each locus the pen is dropped H at ROOTS(J) and a continuous line is drawn through the points | i ROOTS (J+ (ICOM-1) NROOTS) , for ICOM=1, .... INCR+1l. The pen is | then lifted, moved to ROOTS (J+ (ICOM-1) NROOTS) , dropped, and a J i continuous line drawn through the points ROOTS (J+ (1COM-1) NROOTS) , \ for ICOM=INCR+l, ..., 1, completing the locus for the J-th root. This process is repeated for each locus—J=l, «.ss NROOTS—except at those plot points as specified by the array ILOCUS. At these points, for each J the pen is lifted, moved to the next point, dropped, and the process is con- tinued. 1 In addition, subroutine PLOTTER automz.cically chooses i i axes scalings that best suit the data. For a Zeta Locus, the unit circle is plotted as a reference. If the Zeta Plot is i completely contained within the unit circle, axes limits of | +1 are chosen. If it is not, then axes limits of #1.5 are } chosen. If the Lambda Plot indicates the method is A-stable the minimum negative real axis limit is set to zero. y

. Ca ’ . gs |

a = ~

| y For Lambda Plots subroutine PLOTTER also computes the | angle ao for A(a)-stability in accordance with Definition 2.3. The value y for the region 5 of Definition 3.6 '3 also com~ puted and printed. The sutomatic scaling can be overridden | and overplotting—on: Lambda Plot on top of a previous one—

accomplished. | These features have been used extensively in conjunc- tion with an interactive version of subroutine ACCEPT used to i t vary certain "free parameters” in order to "optimize" the P | Lambda Plot. The composite linear multistep methods pre- P sented in Chapter IV were derived by using these properties. OTS) + ] The plotting program listed in Appendix A is limited to a characteristic polynomial of maximum degree 10 in both I A and ¢. This limitation is arbitrary.

he (

4 CHAPTER IV A NEW INTEGRATION PROCESS [§ ’ 4.1 Free Parameters A composite linear multistep method may be written as

) (85,5%ne+3 = hB; , 5¥ne+;’ -0, Jl with i=1, .... 2, (4.1) or, in matrix notation, as ay, x A Yn AES A_g¥n-x’ -niBgl, + Byf,y + oot Bln) = 2 (4D . . . . where Aye B_y¢ ¥n-4° ¥h-3 and K are as defined in (2.27). Each of the % equations (4.1) has 2k + 2¢ coefficients. However, as they must constitute 2 independent equations, only 2k + 2% coefficients are arbitrary. This may be seen as follows: If (4.2) is to be at least of order zero, a solution for must exist when it is used to numerically integrate . { the differential equation y = 0, i.e., Moly + Ay¥py + 0 + Agdng "0 (4.3

¢ v

68 must have a solution for ¥ , or, equivalently, ag must exist. Pre-multiplying (4.2) by apt results in i x, . A Yn An. A_g¥n-x’ i Bl, + Bgtpy + 00 * Bade) = 00 0 & ge - = ak a - he identity . where Ay - Ag Ayr B_4 H A, By» and A, z I, the identity
) matrix.* Hence each of the equations (4.1) only has 2x + 2 arbritrary coefficients. JER —— t *A CLMM when represented as in (4.4) is said to be in canon- 1 jcal form, a term suggested by W. Rubin (Electrical and Computer Engineering Department, Syracuse University). Any given CLMM can be put ‘nto its canonical representation by pre-mult iplying the former by a. The characteristic equa- tion of (4.4), x(%. A), is given as, “ 3 xg, A) = det{ J (Ay ~ xB_je 7) 3=0 1 i whereas the characteristic equation of (4.2), pl(g, A), is given as 4

I plg, A) =det{ J (Ay ~ a_i")

As 3 exists, the latter may be written as

plz, A) = det{agl § (Aj - B_0e5 3) © | |

{ or { plz, A) = det{a;} x(t, N). Hence, within a constant factor, the characteristic equation, ° { and therefore, the stability properties of a method is the same whether or not it is represented in canonical form.

3 is.

i 69 If one requires a composite linear muitistep method i to be of a given order p, then p+l coefficients, per equa- tion, are predetermined. Any remaining coefficients ry then be used to improve "other properties” of the method. i These free parameters, 2k + & - p - 1 per equation, can be . used to i) improve the stability characteristics of the method ; | ii) make it efficient as implemented in an integra- tion algorithm; and d iii) lower the discretization error introduced at each stage of the process. | In the next section a numerical integration process consisting of 7 composite linear multistep methods, of order | 1 through 7, is presented. Free parameters were adjusted to make each method stiffly stable and efficient. ) 4.2 Some New Cyclic Composite Linear Multistep Methods As one of the primary objectives of this research effort was to derive processes that were efficient to use, the class of composite linear multistep methods investigated was cyclic, i.e., in (4.1) o; 4 = 8; 4 = 0 for j > i, or, equivalently, in (4.2) A, and B, were limited to being lower triangular relative to the skew diagonal.* Moreover, as ’ ET Sp — AJ *As pointed out in Section 2.6, this avoids the necessity of simultaneously solving L equations.

Cr : : Be,

these methods were to be applicable to the solution of stiff systems the conditions B_y = 0 for j > 0 and ER non-singular were imposed. These conditions, though not a neceas ty, were set in order to guarantez asymptotic stability at i h = », That this property is achieved can be seen as follows: J As the step-size increases, the zeros of the characteristic polynomial tend to those of p, (2) , as can be seen from (3.8). Now, Pp, (2) is given as

p, (2) = det{ ] a B (4.4)

i Requiring B_4 = 0 for j > 0 yields 0 : py (8) = det {Bgl (4.5) which, as By is non-singular, is equivalent to setting all zeros of py (2) to zero. This guarantees stability at h = =, independent of the coefficient settings. Under these conditions the first of tre iL equations (4.1) becomes

I 65,5 Ynees ~ he) 3 Ypea1 = 0 (4.6) j=-k+1

{ If (4.6) is to be of maximal order k the backward differen-

{ tiation formula (BDF) results [HEN62). By imposing the ] l I | condition that all ¢ equations have at least order k and i

CC. ®

« 7 2 ~—

further restricting each equation to span only (k+l) time points, for efficiency considerations, it follows that i-1, i=1, ..., 2, free parameters exist in the i-t! equa- tion of (4.1), each of whica can be adjusted to improve the stability characteristics of the method. Summarizing, the class of composite linsar multistep methods considered in this research was for inl, ..s0 Ls restricted to [R1]: 5,5 = Bi,3 =0, j> i; [R2]: Pilg ™ 9 3 50 md By,3 ¥ 0:* [(R3]: 5,5 = 0, j < -k+i; and I [R4]: Pj > k, where Py is the order of tne i-th equation. Conditions [Rl] and [R3] lead to the cyclic methods | proposed by ponelson and Hansen [DON71] and used by Chesler and Pierce [CHE71], and Dyer, pierce, Haney, and Chesler [DYE72] to obtain numerical integration processes suitable for use in orbit computation. In their methods free param-= eters were adjusted to improve stability properties and ob- { tain higher accuracy as h + 0. None of these methods are stiffly stable. . In this research effort the criterion used for choosing

i ———————————————————— *As By is restricted to a triangular matrix this is equivalent to requiring By to be non-singular.

5 \ 4 c= ;

i 2 — ~ . : (

72 the free parameters of a method was: Maximize the Widlund wedge angle a of pefinition 2.3. The backward differentia- tion formulas for orders 1 and 2 are, respectively,

Yps1 = Yn id ep (4.7) .

ps1 = Wp "Yn * hy - (4.8) They are both A-stable processes and, therefore, no further i improvement in a can result. The backward differentiation formulas for orders 3 through 6 are stiffly stable, though ] for a < ®/2. For orders 7 through 15 they are not stiffly stable [GEA71b]. [ If conditions [Rl], [R3], and [R4] are first imposed, then it can be shown that there exist k+l linearly independent multistep formulas for each i, from linear combinations of which can be composed % linearly independent formulas of a CLMM satisfying condition {R2]. Sets of such linearly in- dependent multistep formulas for k=3, ..., 9 are given in Table 4.1. Using an interactive version of subroutine ACCEPT* in conjunction with the plotting program described in Sec- A tion 3.5, various linear combinations of these formulas were BE Ad #a brief description and a listing can be found in Appendix B.

P=k=3 1 v2 =i 0 1 121 0 1 Joi, §2 KA wang by 1500000 «11 18 9 21 6 0 0 0 «2 *3 6 #11 0 6 0 0 ) §{ «6 3 21 0 0 6 © 2 9.18 111 0 0 0 6 eed Ne i ge B ’ iii iy-i P=k=4 1 «3 #2 1 +] 11 «23 2 -i 0 1 .25 48 «36 16 <3 | 12 0 0 o © «3 «10 18 #6 1 o 12 o 0 © 1 «8 0 8 1 eo 0 12 0 © «1 6 +18 10 31 o 0 ©o0 12 © I 3 «16 36 -48 25) o © Oo 0 12 I P=k=5 [ oa #3 +2 =i ) 11 =& =3 2 «1 0 1 ’ «137 300 300 200 =75 12 | 60 o Oo Oo Oo © 12 =65 120 <=60 20 <3 | o &¢ © ©0 ©0 © f EUR Sw gt 0 om 0b O° 9 e2 15 «60 20 30 «31 o 0 OO 6 0 © 3 =~20 60 #120 65 12 1 0 0 0 0 60 0 -12 75 =200 300 =300 137 | o 0 OO O 0 60 P=k:=6 5 «4 3 2 -1 0 1 “5 -4 =3 =2 = 92 1 .187 360 =4S0 400 -225 72 =10 1 60 0 0 oOo nn Nn 0 «10 =77 150 «100 50 =15 2 sa 66 © © 9% 9 © { 2 24 «35 80 =30 R “1 0 0 60 0 n 0 0 *q 9 «45 [+] 45 «9 11 [s] 0 5 60 2 n 0 1 8 30 =80 35 24 2 | 0 0 0 0 69 b] 0 A «2? 15 S50 100 -150 77 10 | > 0 o 0 nn é& 0 » 10 «72 225 «400 450 «360 147 | 0 0 LJ 0 0 n 60 Table 4.1: Linearly Independent Multistep Formulas for P=3, cess 9

Big, -

l 3 —~

idl) VY

aia ae ew «8 .« # oo BA @ 1089 9% +0010 S900 -37S 17s 90 40 1 2S , " i “80 +409 1260 +1080 700 ns » 10 » 3 ° M 40 4% ite To mo 0 wm + 8 “Rw | 19 18 Ths ios are ci2e a od 3 amo S a8 ihe coo jo8 mE ceo » : J BB ee “30 70 39 eo 10 HR eR a 3 9 1 os MS -700 ios veo sod 40 HEE HET a 30% NTS amp wd <7Se0 1049 FEZ 233 8 a i p— Cee EET TR ee i i ' sta stp o3iTe0 SAS wisf00 ics comms Mo mgt: Me 8S 5 5 3 a. B03 4700 Vijay [ie Peso -ier0 sas .1e3 13 halk 2 2 2 8% 3 ¢ 8 0 +1338 E00 “hese eivso Seo ceo 8 8. 8 5 #0 0 eo + o ©

3 oe

] 4] tried for each of the equations (4.1) for different L and i orders 3 through 7,* subject to restriction [R2]. Table 4.2 gives the resulting composite linear multistep methods d ~- rived in this manner.’ For p = 2 and 4 the new methods con- sist of three equations (=); for p = 5, 6, and 7 the new : methods consist of four equations (i=4]. Table 4.3 compares the Widlund wedge angle of these methods with that of the corresponding backward differentiation formula of the same order. Table 4.4 provides a corresponding comparison for the value y for the region S, of Definition 3.6. As can be seen from Table 4.3 the Widlund wedge angles for the new methods have been significantly improved over the corresponding ones for the backward differentiation formulas. Table 4.4 indicates that for all of the new methods presented the half plane, defined by the parameter Y of Definition 3.6, J f wherein a composite linear multistep method is asymptotically stable has been enlarged when compared to the corresponding backward differentiation formulas.

] ————————— { #*An unsuccessful attempt was made to obtain a stiffly stable , order 8 composite linear multistep method. Failure can be / reasonably attributed to the inadequacy of the "manual® search of the free parameter settings and not to the non- ( existence of such a method. tin Table 4.2—printed by subroutine ACCEPT—A(I,Jj) and i ] B(I,3), J = k+l, «coy Zs correspond to ¢; 4 and Bi,4° respectively, of (4.1). .

: ‘ Ie

£ #5 ~— “5 76 Pek=3 1 1 2 3 Atlee) 2 [] 0 Atls=1) bd v2 0 Ati» 0) wil 9 ° Atl, 1) 11 18 9 Alls 2) 0 1 12 Ally 0 1] 3 Blls,=2) J 0 2 Btls,=1) 0 0 0 > Btls 0) 0 ° 2 8tl, 1 6 0 os Bil, 2) 0 6 Btl» 3 ] 0 2 P=k:=4 1 1 2 3 Atls=3) 3 0 0 Alls=2) ~16 3 0 Atl,=1) 36 «16 11 Alls 0) wad 36 ol] Atle 1) 25 «8 216 Atl» 2) [2] 25 ~272 Atl, 3 0 0 93 Bitls3) [2] 0 [] B(ls2) [1] 0 0 Btls*1) [4] 0 0 81, 0 [4] [4] 0 8(l, 1) 12 0-8 . 81, 2) 0 12 cas Bil,» 3 J [2] “3 Pzk=5 ' i 1 1 ? 3 » Algs=4) ei2 [1 0 0 Alls") 75 e12 ] 0 Alls=2) =200 75 118 0 Allse1) 300 =200 735-133 Atls 0) #300 300 =1940 780 Atl» 1) 137 +300 2980 ~1680 Atl, 2) 0 137 ~3030 5470 Fy Atls 3 0 0 1373 +5595 Alls &) ° 0 oD 1158 Blls=4) [4 0 [1 [4 B(l,=3) 0 [] 0 0 Btl.*2) 0 0 [4 0 Bils=1) [4] [ [4 [4 Bil, » 0 [:] 0 0 B(l, 1) [1] ° -60 30 Bl, 2) 0 60 0 ~1860 ( Btls 3) 0 [ 600 =1530 ’ Btl, &) [] 0 0 600 Table 4.2: New Cyclic Composite Linear Multistep Methods

: \ y

B ” Pzk=6 1 1 2 3 Bl H Atls"5) 10 0 0 0 A(ls=d) 72 202 e 0 Alls=3) 225 +1455 195 0 Alls=2) «400 4550 ~i299 28% Atls=1) 450 ~8100 4340 «2039 Atls 0) #360 915) «7540 6225 Alls 1) 147 «7277 8505 <=10360 Alls 2 0 2930 «74a5 18455 4 Atl, 3) 0 ¢] 2944 «14865 Atl, &) 0 0 0 2299 Bl1,*5) 0 0 0 0 Btls») 0 0 [¢] o ’ | B(1,=3) 0 0 ] 0 Btl,=2) [°] 0 [+] 0 B(l,"1) 0 0 0 0 B(1s» 0) 0 0 0 0 Btls» 1) 60 w60 420 180 Bil, 2) 4] 1200 «60 -4080 Btls» 3) 0 0 1200 %680 H Btls &) 4] 4] 0 1200 P=k=7 1 1 2 3 i . AlL1276) =60 1] 0 0 E A125) 490 =60 0 0 A(lsaws) «1764 490 «210 0 Atls=3) 3675 1764 1722 «774 A(Is2) =4900 3675 =6235 6349 Atls=1) 4410 =4900 13100 ~22988 A(ls 0) #2942 “al0 «17650 48160 ( All» 1) 1089 «2940 17710 ~66290 { Alls 2) 0 1089 ~11297 68i5° Alls 3) ] 0 2860 =42364 Alls &) 0 0 0 9748 Blls=6) 0 0 0 0 B(1,°5) 0 0 0 0 | B(ls=%) le] 0 +) 0 ’ { B(1,"3) fe} 0 0 0 B(l,=2) 0 0 0 7) B(l,-1) 0 0 0 0 8tl,» 0) 0 0 0 0 . Bel, 1) 420 0 600 840 B8(l, 2) 0 #20 ~1860 =2100 ’ B(l,» 3) 9 0 1200 -8400 Bll, 4) 0 0 0 4200 | rable 4.2: New Cyclic Composite Linear Multistep Methods (continued) |

EY VERN

] CE \

‘ ® 2 —~ — \
NEW 1 ORDER METHOD BDF 3 89.42% 86.03° 4 80.68° | 73.35° 5 77.48° 51.84° . 6 63.25° | 17.84 |

7 33.53° | » | ! FREITEPIREIR ) I #The order 7 BDF is not stiffly stable. i Table 4.3: Comparison of Widlund Wedge Angles

[ ORDER METHOD BDF 3 -0.0048 | -0.083 4 -0.24 -0.67 5 -1.4 | -2.3 3 -2.9 | -6.1 ‘ 7 -10.2 | * I—— #The order 7 BDF is not stiffly stable. Table 4.4: Comparison of y-values.

~~ J ~

: ©. K #3 —

’ : (

The effect of these improvements and the introduction i of a new order seven method is to permit a numerical integra- i 4 tion process based on these new methods to increase th. inte- l gration step-size at a higher order much more rapidly than | processes based on the baciw.rd differentiation formulas, | for a larger class of problems. This effect will be espe- cially evident in stiff problens wherein the eigenvalues of | the Jacobian matrix are complex. The Lambda Loci for the backward differentiation | formulas of orders 1 and 2 are shown in Figure 4.1.* In | Figure 4.2 the Lambda Loci for the new methods, as well as : | those for the corresponding packward differentiation form- | ulas of the same order, are shown. The Lambda Locus for the new order 7 method is shown in Pigure 4.3. The effect | of the optimization criterion used in selecting the free ; parameters—that is, maximize the Widlund wedge angle—is [ especially apparent in the Lambda Loci for the orders 5, 6, ' | and 7 methods. The “humps” in the Lambda Loci for the new methods are almost all tangential to the rays forming the [ boundary of the set S of Definition 2.3. : These new methods, together with the backward JE ———— *An identifier consisting of 2 letters and 3 digits, is ! used to label each plot. The three digits are, respec- tively, the values for p, Lt, and k corresponding to the methods for which the Lambda Loci were plotted. .

CT. 3

Bg \

ruc, +

SEE EEE

/ : a | /

l R12,

« i LY

ae Ca

Mo rd

Figure 4.1: Lambda Loci for Order 1 and 2 Methods }

Cre : a= :

Sa

~.
i BOF , NEW METHOO \ \

: PURE ESP Er IPE

es oe | Tramp -6 R434C I 1 vg Ln : NL \ BOF N N
NEW METHOD \ \ ' -1 ) 3 / / rd / nrg RE a

Figure 4.2: Comparison of Lambda Loci Between New Methods and BDF for Order Three Through Six

, ; \

— \ x

RSAC

— a hp She OF % \

NEW MET? \ \ . FPIPIIPSPIPIPUPUPITIrR Seer ©

"4 /

cg” ”

I -15 RE46C 8 SR Tn Ph .

BOF 5
. f NEW METHOD \ !

= Lr

Nd py

( AJ Figure 4.2: Comparison of Lambda Loci Between New Methods ’ and BDF for Order Three Through Six (continued)

d > ~ 83 R747 2

1 sass | Ne Pd

[ Figure 4.3: Lambda Locus for New Order Seven Method

[€ “w (4 \

_ ® 3 ——

differentiation formulas of order 1 and 2 form the basis for a new variable order, variable step algorithm described in the next chapter.

d #5 —

CHAPTER V A NEW VARIABLE ORDER/VARIABLE STEP-SIZE INTEGRATION ALGORITHM 5.1 Introduction A computer program embcdying the cyclic composite ’ linear multistep methods described in Chapter IV has been written for solving a system of first order ordinary differ- ential equations. In the following a description of the com- puter program organization and logic as well as results ob- tained from some sample problems are described. 5.2 Computer Program Organization and Logic The computer program, called DIFJMT, was written as ’ a subroutine to be called by another program. DIFJMT, modeled i after Gear's algorithm DIFSUB {GEA71a), is a variable order, variable step-size algorithm. The step-size and order— -! manner as to achieve as large a step-size as possible while { not introducing a single step error greater than a user specified value. The computer program has four major sections which ’ will be discussed individually. They are: i) Starting procedure ii) Implicit equation solution ’

bh : |

y 2 —~

iii) Error control \ iv) Step-size and order selection logic An organizational flow chart of DIFJMT is shown "2 / Figure 5.1. A detailed flow chart and a computer program § listing can be found in Appe ix C. The k-th order composite linear multistep method is ’ i given as ;

1 I (8), ¥nges ~ ny s¥nees’ = 0, ju-k+i i with i=1, ..., %s (5.1) i where the a ,3'® and Bs,5'% are as in (4.7) and (4.8) for A k=1and k = 2, respectively, and for k = 3, ..., 7 are as I contained in Table 4.2. I The process proceeds as follows: The backpoints : Tabeiorl® °° Yni and their first derivatives are assumed known. On entry to DIFJMT the index i is set to 1 and (5.1) is used to solve for y ,., and Yoee1? i is incremented by 1 and (5.1) is used again. The cycle is repeated 1 times— { 2 being a function of the order of the method used. The error test is then performed on the last term computed— Ypeeg—and if the estimated error is less than a user speci- fied value the cycle is accepted; otherwise it is rejected. In either case control is transferred to the step-size and order selection section wherein the largest step-sizes, ’

| | | Co ;

I 87 Perform initialization and start-up procedure is jequired. p | Predict and correct | L new mesh point solutions. Perform error test on last point computed. i

Set step-size and order for next cycle, if error test passed, or for redoing present cycle, if error test failed. ' Pid | AP. No Set up to redo i —— present cycle.

: : . : ’ Figure 5.1: Organizational Flow Chart for Subroutine DIFJMT |

Cr -

l + -

consistent with a user specified maximum allowable local error, at the current order, one order lower, and one order higher* are determined. If the error test indicates { llure the entire cycle is repeated at the largest of these step- sizes and the corresponding order. Otherwise control is returned to the calling program and on the next entry the largest step-size and corresponding order previously com- puted is used for the new cycle. The following sub-sections describe in detail the various segments of DIFJMT. The solution of the implicit equations (5.1) is treated first. Subsequent sub-sections describe the other main segments of the computer program. 5.2.1 Solution of the Implicit Equations’ Each of the equations in (5.1) is implicit. The con- = ventional iteration for solving (5.1) uses Picard's succes- sive substitution method which proceeds as follows: If the

A . (0) : py predicted value for Yno+i is Yo94i the iteration JF -—— *If the error test indicated failure no increase in order is attempted. Logic is also built in to limit the order to a maximum value of 7 or a user specified lower value, namely, the value of the parameter MAXDER in the variable list of the calling sequence for DIFJNT as given in Appen- dix C. The treatment in the first part of this s ih-section is similar to that of Gear [GCAG7].

re, % IN 1")

{ eS ~

(m+1) o (m) y %5,i¥ne+i ws hy sEWngsq tarsi

for any norm, in a neighborhood of the solution Unless h is of the order of the reciprocal of the largest eigenvalue of 2£/3y, (5.3) will not be satisfied. This clearly restricts [ the maximum allowable step-size that can be used. A Newton-Rhapson iteration can be performed to solve ‘ (5.1) without restricting the maximum size of h. The pro- i cedure is described below. :

Equation (5.1) for each i may be writcea as

Pl 2% )) (05 5¥n2+3 hB; s¥nges)” (5.4) '* je-k+i

x \ H

l + —

where the right hand side is a sum of known values and hence is a constant. Equivalently, a value for Yneei is desired such that Glypg45) = © (5.5) . where dod { z - iil ’ - hB v 5 ce 51 1 (0; 5¥n2+3 8; ¥ness) (5.6a) rT jm-k+i

is a constant and Bia = - ht: f(y t GY 045) = Ynesi ~ PE, £(¥noeir Enos!” (5.6b)

: (m)._ . ; a Ei. oc 28 (m+1) Assuming Yne+il® an estimate of the solution of (5.5), Yno+i is to be sought such that ety iT) = C. Upon expanding [ (5.5) in a Taylor series about “her and keeping the first two terms, it is found that aly™) + wry ™D - y™) = ciy™) =c, (5.1 where the subscripts "nf+i" have been dropped for notational simplicity and where w = 2¢| (5.8a) m - 3! (m) : ™y

— .

! 3) ;

{ eS gy

or, from (5.6b), Bi,i LL" =1 = . "WY I? (5. b)

I being the identity matuix and .

terms, the desired result; namely, Bias y yo opm dd lee) - aj, (5.11) i,i ’

oe 1, . \ “

where a® - pe ™D 6) + nay™ - y™1, (5.12) In order to obtain a recurrence relation for atm , substitute y (5.11) into (5.12) with m + 1 replacing m in the latter; thus, aml) ope ™ ey Bi,i 1 (m) (m)

i Then, adding and subtracting wlmey™,6) - a™) and re- i arranging terms gives the desired result; namely, 1 am) | gm, gm lpe ™ g) - al™), (5.14) In order to start the recursion process—obtain y@ y and a?) consider the k-step explicit formula

of,i Yneei * I (af,5 Ynees ~ PBL5 Yates) ® 9 j=-k+i with i=l, ..., 2. (5.15) Now, 4] may be taken as

¢ #* —~

Ee a * - h8t . ¥ Yne+i C5 394 I (e],5 Yages = PBL, Ynges)- (5-26) 'd jmek+i da (0) is chosen so as to make (5.1il) correct for m = J, ni+i 8 . . _, 10 sing (1 hedth (0) correct” meaning Yne+i Y.ag4i® Adding [I . INYos4i

to both sides of (5.10) with m = 0 and setting the left hand y side to zero, as the condition y (1) =y (0) is being imposed, nL+l ni+i results in

{ Bs 4 a, (0 0 = —ti[hfly (0) 1) = anit y »0 95,4 ni+i ni+i Bs,i ni+i

Bi I (85,5 Ye me Ynges) i r* jm-k+i By substituting (5.16) and then examining (5.11) with m= 0 and y (m1) = y®, it is seen that choosing q

(0) & 1 by B - a¥* ) Gnesi = GF 1B; 5 I [lef 5 25,5 ~%1,5 %,1 Ynesi ‘ " re jmek+i Pe * - B* v 1 hiaf 4 84,9 81.3 5,4) Yneei’ (5.18) makes (5.11) "correct" for m = 0. Note: Both y (0) and ni+i 1 2 are linear combinations of earlier—known—function .

values and derivatives. Further, if +1 as given by (5.16) oan rE p (0) _ , (0) = satisfies (5.1), then, from (5.11), a si = hE(Y, psi” topei! - . (0) a a aa 2 i hy g+i® Therefore, doo+i can be identified as the predicted .

value of hy pei” Moreover, when convergence is achieved ’ $ (m) * m .. . . % Lye] = hy 43 that is, just as y,.; 1s considered the

for i=l, ..., &: { 1. A predictor formula for suitably chosen 5 3 and

. : ((}] c

  1. The Jacobian matrix J is evaluated and used to compute By 5 _.- wlez-ndd 57th

  2. Form=1, 2, ... (5.11) and (5.14) are recur- sively used to solve for Yne+i and 4 osi” the latter being taken as hy 94i° As the predictor formula used to obtain y 20 , and : { hence ot. does not affect the stability properties of the method-—since each of the formulas (5.1) is iteyated to con- ’ vergence—the sane form of (5.15) is used for each i. More precisely, in (5.15) for i =2, soep and j = -k#i, cece i,

95

l $1," 82,31" (5.19b) Then, for i = l—the initial predictor—(5.15) is chosen to be of order k and of the form

» | J wl is not re-evaluated.* Convergence of (5.11) and (5.14) is determined as follows: 1f, form m = Ay ¢™ ¢ [IW time ™, £) - a™j|| < up, (5.21) where BND is chosen as’ e/(2k + 4), € being the user speci- fied accuracy, then (5.1) is said to have converged. For m > 1 the inequality i s™ » mina, 26™/s™ 1) < BD (5.22)

] ~ 1 —————————————————— *There seems to be no explanation for this anomaly. However, this backs up the experience of Gear [GEA68]. ’ trhe Euclidean norm is used, as in (5.32). Sonis choice was made on the recommendation of Gear [private communication] and works satisfactorily.

function evaluations.* In order to save computer time a maximum of 3 itera- tions of (5.11) and (5.14) is performed. If convergence has not been achieved by the end of the third iteration it is assumed that (5.11) and (5.14) are not going to converge. ~ At this point wl is re-evaluated. If the iteration fails ’ a second time at the same mesh point’ the step-size is de- creased by a factor of 4 and the entire cycle is repeated. There is no theoretical justification for the choice of a maximum of 3 iterations. However, numerical results indi- 3 cate that, if 4 or more iterations are permitted, the number [ of times the corrector fails to converge 1s not appreciably less. if only 2 iterations are permitted, the number of < times convergence failure occurs is considerably greater, resulting in significantly more evaluations of wl and more function evaluations. 5.2.2 Error Control ry - With reference to (2.24), Lyly(t g.4)r h] is the . i | amount by which the solution of the differential equation . os £2ils to satisfy the i-th equation of a composite linear ERE — *This look-ahead feature is a variation of a suggestion made by Gear [private communication]. TNote: 1 > 1 mesh points are computed for each call to DIFJMT.

multistep method. As such it is a measure of the error $4 ) | introduced by the i-th formula. Since the solution of the j-th formula of (5.1) equals the truncated series of the . true solution through terms in n*, it is known as the local truncation error.* Now, for (5.1) é IN k+l _ (k+1) 1 Lyly(tpo, de B) = Cy pn BY (Epges)s (5-23 — k+2 . where = means to within terms of O(h™ °) and where

1 k+l Kk C xs = TFDT I G3" Tey, 5 (k+1)378; 4) (5.24) J==k+i Further, Dy; k+l’ tiie local discretization error constant for {

the i-th formula of (5.1), is given as D x41 = Cs ,x+1/%4,i" Ld =1, coop 2. (5.25) Now, for i = 1, «..s %, SH |

’ S k+l _ (k+l) Yt pei) “ Topes "Pra BY (thgeq)- (5-20)

JE 2 : *Strictly speaking this is true if, for a k-th order method, the true solution is (k + 1) times continuously differen- k+1 : tiable and Yno+j+i YE page) + o(h ) for j==K, «ss, =1l. Although in general this cannot be shown the technique used

99 Equation (5.26) forms the basis for the error criterion used . 4 | in DIFJMT. It proceeds as follows: At the completion of ¥ each cycle a is computed, where V is the backward difference operator defined as Vx, z vik, - i WOO for 3 = 1, voes (5.27a) 3 and vexn H] Xe (5.27b)

k+1 o pk¥l, (k+1) 5 VC Yne4n hy (1) (5.28) for some 1 lt t ], the quantit |= | is - i ne+2-k+1’ ‘ne+2’’ q yl Yne+t compared with =/Dysrr where ¢ is the user specified maximum J single step error permitted and F

Der = i=l,nanst Pinal? = (5.29) re k+1 ’ 19°" Yagael < &/Pkiy (5.30)

the cycle is accepted.* Otherwise the cycle is rejected, the “. ————— } *To obtain (5.30), it was tacitly assumed that the (k+1l)-st derivative remains constant over the interval [t ssp-k+1" thee) in general this is not true. However, ie is a reasonable estimate for pktly (k+l) (t) anywhere ! in the interval [t ,.o ,.1/ tigen)” As the . future points ! ’ |

step-size is decreased as outlined in the next section, and $ the cycle repeated. If the error test fails on three succes=- sive tries the order is set to 1 and the cycle repeated being treated as if it is the first call to DIFJMT.* This proce- dure is followed for two reasons as follows:

  1. It is possible that the error test failed due to

| the eigenvalues of the Jacobian matrix associated

| with the differential equation being in an un~ stable region of the A-plane. Since the order 1 i method is A-stable—as can be seen from Figure 4.,1—it is used. 2. The error test may have failed due to accumulated | round-off error in Yne+s” j= ~k+l, «oes 0. Hence the integration process is restarted. | If the error test failed and the step-size cannot be decreased’ then the order is immediately set to 2 if it is greater than of. ipo | solved for by (5.1)—for all order methods—span (tog+1” tos) E ' [to o4p-k+1’ LP using Dps1 as defined by (5.29), the / the error test (5.30) yields a reasonable estimate of the A, maximum error introduced at any of the ¢ mesh points ¢ Us togeit i=1, «ov 2. *Special procedures followed during start-up are described in Section 5.2.4. ’ » » Tne step-size is not decreased below a user specified mimi- | mum value, namely, the value of the parameter HMIN in the i variable list of the calling sequence for DIFJMT as given I in Appendix C.

H 101 2, or to 1 if it is already 2.* 1f the order is already set $4 . [ to 1 the cycle is accepted and the user notified. ¥ As DIPJMT can be used to solve a system of N first order differential equations the error test performed is not really that given by (5.30) but rather that given as k+1 3 H Dy lv Yngee!l < Es (5.31) ! where®

k+1 2 _ k+l p HY" Yngenll I HT pgp) 5/eg)" (5.32) i -~ (75+ : . > n— . I v Yne+’ § being the j-th component of the N-dimensional 4 k+l S Th - vector V Yna+e and wy being a weighting function for the I j-th component . ¥ JE —— Note: The order 1 and order 2 methods are A-stable. ’ tany other suitable norm could have been used but (5.32) is ' easily implemented in a FORTRAN program. Note: If (5.32) y is used it is only necessary to form N72 an] 1 and ' : 2 - p { compare it to (e/Dy,,)"- .

102 The discretization error constants for the methods H used in DIFJMT are given in Table 5.1. § 5.2.3 Step-Size and Order Selection Logic The step control mechanism is activated in each cycle following solution of the algebraic equations (5.1) and ap- i plication of the error test (5.31). The step-size to use in i the next cycle or when repeating a rejected cycle is cetinated to be nh, h being the current step-size. The value n to be used is determined as the maximum of ny_q+ Myo and Ny qr y where

( u 117 GG+1) EIN PIA" wl +l Ing+2!!

I with j=k-1, k, k+l, (5.33) Note: Y is a constant such that 0 <y <1, and k is the current order. Also note that if the step-size nsh were f used with an order J method and if the error were exactly : proportional to hIl that is, 93% y oll remained con- pr / stant—the error test (5.31) would be satisfied, and satis- B |» fied exactly if y were 1. As 1199y oul] is not usually constant the factor y is used to provide a "safety" margin.*

i FE ——— *a value of y = 1/1.2 is used in DIFJMT. t/

I 103 4 Ld < 3 : © oo o >

oN — | [+] i A NON a | y-1 3 a oo oo | + .. mo = | Q X NN ww mm | = |! | . | | |]

” [=] © “ ~ mm © +»

-— [)

od = nw ~~ nw — i NN mM Oo o vs N = = Nn © = 2 SS NN NN NN o 3 : m 8 © A © = h I = A NN ™ Pe) J a f} 1 — 1 < H ’ -

— 8

nw ~~ mw ®

—- Fel BN

S » od PY - wn ©

i 104 Now, given ny for j=k-1, k, and k+l, the next step-size and . 2 | order would be those associated with the largest of the three step-size ratios.* The computer program listing given in ! Appendix C does noc increase the step or change the order— assuming the error test passed—if n < 1.1 since it is felt the increase is not worth the extra computer time required X to interpolate the mesh point values at the new step-size. ; As the step control mechanism requires quay WOR in order to determine Nyy requiring information from k+3 dif- ferent time points in its evaluation—and as only Y,, > J=-k+l, ..., 0—information at only k different time points— i is available at the start of DIFJMT, the solution at a mini- - mum of three mesh pcints is computed at each carn.’ More- I over, as both the error control and step-size selection logic [ are dependent on the backward differences of Yne+o! the i-th : equation of (5.1) is formulated’ in terms of jg TPN instead So STR 4 *1f the error test failed ny, is not computed and no attempt ' is made to increase the order. J Tie an order 1 or 2 method is used the same equation—(4.7) ‘ ! or (4.8), respectively—is used for all three mesh points. . For higher order methods 2 > 3. $n terms of backward differences (5.1) can be written as k ! (@5 6" ness BE n8; o¥ngsi-o - 9

with i=1, ..., 2.

Il 105 : of Y,9+i-0 for 0=0, ...., K. . i when the step-size is changed the present value of time the remains fixed, but the past values of time totes’ Jo-k¥l, «ev -1, must be adjusted to reflect the new step~ size h = nh by which they differ. Note: kx denotes the possibly new order. Thus, it also follows that it is ne- 3 cessary to compute k - 1 new values of the solution, Yne+i’ i je=k+l, «oer -1. The value of hy needed for the predictor is trivially adjusted by introducing as a multiplicative factor; thus 1 hyn, = n(hy,)- (5.34) J—— i where, Zor 0=0, «++» k and i=l, ..-¢

= o v %i,0 = (-1) ) (“ay 5-0"

and ~ Bio Bi,i-0"
The 85,5 ° and B;,5'° (i=1, ..., * and ju=k¥i, coer i) are as in (5.1). Formulating (5.1) in terms of packward dif- ferences does not affect the solution of the implicit equa- ay ya tions as described in Section 5.2.1 except that VY ead” ¥ 2 i.
¥ | g=1l, «.-s k, must now also be corrected. Note: VY na4i - . ed Yoe4i® In order to save computer time, the computation of o { VOY gait ¢ > 0, is performed only once convergence has been i achieved. » >

i 106 Two cases need be considered: . % i Case 1. The error test passed on the previous cycle i and n > 1.1. Case 2. The error test failed and the current cycle is to be redone. case 1: A k-th order interpolating polynomi 11 is used to b! compute the k-1 new values of the solution Yoga? j=-k+1,

wees =1. In order to keep the error introduced by the in- terpolation procedure to a minimum the k+1 values required i to define the interpolating polynomial are chosen at k+l consecutive mesh points. These mesh points, selected from i ____among the k+f mesh points available from the previous - i cycle, are chosen in such a manner that, for a given index . j, they bracket t , + jh. More precisely, an index v is I chosen such that, for a given js J=-k+l, ..., 1, te jh € [t, - vh =kh, t , vh] = I, (5.35a) for some v = {0, ..., 2-1}.% Equivalently, an index Vv is p— chosen such that for a given j, j=-k#l, ..., -1, ’ v <=jn £v +k. (5.35b) J——— 3 A . At the start of DIFJMT only the values y (n-1) 243° j=-k+l, ..., L are available from the previous call.

i 107 Note: Using this procedure, for differing values of J, v H may assume a different index. Using Newton's backward dif- i ference formula* [HEN64] an interpolating polynomial which interpolates ¥ ,_ ,.7 y==k, ..., 0, is used to compute } Ynpeie Note: If (5.35) is satisfied then’ A vy k+l, H Yne+j = yl, + jh) + oth Do (5.36) 1 Requiring (5.35) to be satisfied restricts the maximum allow- i able value of n to 1 +k - D/G =D. Npax = (2 k 1)/(k 1) (5.37) Table 5.2 shows the maximum allowable values for ng_;, My» { I and n,,, subject to (5.37). size is as follows: J Step 1: The factor by which the step-size is changed, ' n, is set equal to max {ny} as given by (5.33).

1 { 4 (a { ] 2 #As DIFJMT operates on Voy o.:, 0=0, .... ki $s oese Bs

this form of the interpolating polynomial is used. TEquation (5.36) assumes that the values Y ,_.4, 3F€ exact ’ or at least accurate to order k+l in h, and that y(t) is 4 (k+1) times continuously differentiable for t € I, As h the error test was passed on the previous cycle the values | ¥og-vey BFS, at least locally, accurate to order k+l in h.

i 108 | ] = SE EEEGEENEE BESPNCIIT [x] : Me-1 | x | "ker | — — 1 3 | * | = | 3 2 3 | - | a | 2 3 3 5 s/;2 | 5/3 4 3 3 | 2 | 372 |

5 4 8/3 { 2 | 8/5 we | } 6 4 9/4 | 9/5 3/2 7 4 2 | 5/3 * | 2 *These values are not computed. Tprogram limited to 104. maple 5.2: Maximum Allowable Step-Size Changes aA 4 A

¥ 109 If n < 1.1 the step-size and order are not i changed and n is reset to the value of 1. 3 Otherwise, the new order k is set accord- ingly and n is limited as given by (5.37). } Step 2: Using (5.34) hy, is computed. Step 3: The index j is set to -1 and, as h > h, 8 vy is set to 1. [| Step 4: 1f v does not satisfy (5.35b) then Vv is set to 2-1.* H Step 5: The Newton backward interpolation formula, i interpolating Y ,_ 4.’ y==k, sess 0 is used to compute Y ,.4° § step 6: The index j is decremented by 1 and if j 2 -k+1 Steps 4 and 5 are repeated. ; I Otherwise the process terminates. case 2: The procedure followed when h is decreased is simi- lar to Case l—an interpolating polynomial is used to compute } Yne+i for j=-k+l, «eer -1. Note: Here, ~ denotes values to be established for repeating a cycle, and * denotes values N—— p established prior to performing the first of the rejected

) g cycles.’ The solution values used to form the interpolating :

1 Shae a i *An examiniation of Table 5.2 will show that if v = 1 does not ’ satisfy (5.35b) then v = 2-1 will, for all k. | f ¥ TNote: The ~ values reflect any interpolation due to a pos- h sible step-size increase.

§ 110 polynomial are those available prior to the first rejected NE 4 H cycle*—that is, the " values. The order of the interpolat- ing polynomial used is k. Note: k < k. If no step-size increase was performed then the Vy, o=0, ..., k, are available to use in Newton's backward difference formula, as in Case 1. If the step-size was initially increased then = only the V¥per o=0, ..., k-1, are available. In this case, PS - I together with hy ov a k-th order interpolating polynomial is constructed. It is readily verified that such an inter- i polating polynomial is op) = ag + au + ajulu + 1) 4 oer datulu + Deru ek = 1), (5.38) where us (t~- tag) Me (5.39%a) a = -L ¢% , for o=0 k-1 (5.39) o cl nt’ Fame. . ‘

5 and fod ap = [hy,, - I (v=-21r1a))/k=-1t. (5.39¢) 3

Note: &(u) satisfies EE — { *These values are retained in a special work area for this purpose.

11 oj) = Yne+j’ with § = k+l, ..., 0, (5.40a) [| and $(0) = hy (5.40b t Ye . ) i Hence, $(u), with y = -k+1, ..., -1, can be used to com- pute Yagtn with the resulting values accurate to order k+l B in h.* 5.2.4 Start-Up Procedure Starting is almost automatic. The value of hyo is computed from a user specified initial condition Yor a user specified initial step-size h, and a user supplied subroutine I to evaluate f(y, t)—the right hand side of the differential equation. This is sufficient to allow the first order method to be used. From there on the error control process and the step-size/order selection logic, functioning as previously described, maintain the normal cycle processes. BE *The assumption that the values y_ ._.. 0=0, ..., k-1 and _ ~p - - “... ( hyn, are accurate to at least order k+l in h, and that

B 12 5.3 Sample Problems and Comparisons H In this section five problems are presented which were \ H solved using DIFJMT. For comparative purposes these problems were also solved using a computer program which was a slight i modification* of Gear's algorithm DIFSUB [GEA7la). Gear's algorithm is based on the backward differentiation formulas P [HEN62] of order 1 through 6 and uses the Nordsieck vector [NOR62] as the means of storing information. The Nordsiech vector a, is defined as } . n%y, a ; I : a, = (yg, hy» —37 + poey Sa] o (5.41) I where h is the current step-size and k is the current order. This choice was motivated by the ease with which one can predict Yael and can change the step-size h. However, it has been observed that the choice of the Nordsieck vector, az implemented in Gear's algorithm, does affect stability properties when h varies rapidly [BRA72a). This undesirable ' — P ST *The modification involved removing all coding in DIFSUB that . 7 5 corresponded to the MF = 0 and MF = 2 options in the calling : sequence. The former permitted DIFSUB to use a variable order, variable step Adams-Moulton algorithm. Note: Adams type methods are not stiffly stable. The latter option per- mitted DIFSUB to evaluate the Jacobian matrix by numerical F Ba methods rather than calling a subroutine to supply it. As : neither of these options affect the basic properties of the LF method but only add to computer run time, they were removed.

¢ ] As

13 property was diminished in a third computer program which i was a slight modification of DIPIJMT—referred to as DIFBDF. | This latter program is identical to DIFJMT except that all t formulas in (5.1) are identical* and correspond to the ' backward differentiation formulas. Note: As the order seven BDF is not stiffly stable the maximum order was lim- ~ { ited to six in DIFBDF—as in DIFSUB. i The five problems and their solutions are- contained in Tables 5.3 through 5.7. In addition, each table includes the properties of the solution obtained by each of the three algorithms as a function of the local error requested. The 1 maximum global errors listed were established for each solu- 1 tion as the maximum of the computed local errors at all mesh points, where, for a system of N first order differential ¢ + equations, (y 0, = (g(t ))y i (local error) = ._ a —ni 1 = , (5.42a) n i=l,...,N (w ); ~Nn 1 and ’ = max 1 3 i (wy)y = max{l, 5.0 "a Pgydglde (5.42b)

P + ——————— 8 { , 7 *Observe from Table 4.2, for i = 1 the formulas of the new I methods, for all orders, are the BDF of the corresponding order. Hence the modification was readily implemented.

tRecall: All three algorithms estimate the local error ) using the norm as defined by (5.32). (F

114 i Ty A B...... io een _

vio = 3 |

CoE — BR ping al grror Bound aw 1077 wt | 107? 1w0™ wt wi? ST CC a pry—p— of Function Cvaluations 182 pLlY 01 | 09 ni 90 28 of Jacobian Lvaluations 1 10 19 | 3s 1n n ls inal Step-Size (sec) 0.3% 0.304 0.265 0.210 | s.160 0.126 |o.005¢ o tae fain = 35°% ’ . 1 [ 10 2 re | p—" ION SE —— ! far Na Re RE —————

imum Global Eres TY RY I) PERT dl ERE Ud CE RE ad Lic wi? ' of Function sealuations | 134 tN 22 » | 3 oe “ of Jacobian Lvaluations LJ Ld . ls | * . . { inal Stop-Sise (sec) 0.505 0.365 0.202 | 0.10 s.100 0.3 som Dn f Pu Time (min * 107%) . ? 0 | 2 | 16 3 n p Pr | SSE IE — E——— Py ——— : ay 4 ximus Global -rror 3.06 = 107 3.15 = 107? | 1.35 = wt 2.23 « 10”? 9 «0 Sl PRE Cy 6.60 w of Time Steps © » am 140 mn £31) ne J of Function Evaluations 138 193 206 264 354 om [$1] of dacssise Cvalsstioss | . ‘ ‘ 1 ‘ p N inal Step-Size (sec) 0.554 0.326 | 0.246 | 0.200 | 9.240 0.114 0.0845 7’ Table 5.3: Problem 1 oF

ee e—— [] SirTERENTIAL GCUATION worm — EE ————————————————— p j= wwoy + oy Sp — vien —.

3 rs ware 1-s0le 2 3 11) | we) = NEIEERY | =) vellr aa 1] | Jone TOE 111) | ’ £00) rs “00 © 0 0 1 | 1-1.000e" 0 0-800 0 0 | 's | eo ow oo | | l eo © 9-0.00] | | 2 [ | s | 5 | 2 | | | 2 *2| | ve and x = 9 i > J

5 =) { PROPERTIES OF TEE NUMERICAL SOLUTION Local Error Bound 10 10”? | 10 | 0”? | 207% - we H SOLUTION VIA DIFINT - Maxim Global Error URRY ORT PR PR PRE CRE fbi § of Time Steps 9 | m oe | e90 as | nn 1641 {of rumction Cvalustions | 58 ne 1s [ a72e 1903 | 26m wwe |

—— I SE—— SE —— SOLUTION VIA DITSUS wrimus Global Error 2.48 = Ud 3.98 » 10”? | 2.08 = wt | 6.35 » 207° “nm ad . en. 10?

—— — - —T Ton] —_ Maximum Global Error | 2.20 = Ud 3.19 = 0 7 2.48 = wt 3.97 » 10 » 08 = 107° as2 = 0207-0 - "of Gustion Stanton | 9 po | 10s 1394 10 2026 nm

Ey | DIPFIRENTIAL COUATION 1 souimion p gous eve git) = oie wien wre [* [ -21s,08, 1 -2 | : - = \ wea | 1% [2] | 201-3 08 ' where | > 8 a1) | ate) » EE | TL a IEEE AY { i ™ panne [sal i | rE 1 1 ° 0 1 | dad J \ ae oo | | sel 0 oa 0 | PPT Ep—— J Leo oo ~0.001 | 3 ana 1 r | uy = e719 (5 ain 10 + 10 cos 10¢ \ ra - 4 [=] |

— {=| { | ge a maze v9 | { 3 1%] | i : | { L = =) PROPERTIES OF THE WUMERICAL SOLUTION a1 Error Bound 107 rd | 10* | 207? 10730 wi 173 WRaERER bP 8 | mum Global Error 7.11 «207 | 20-2 = 10”? | 20.2 = 307° | 25.0 « 207 | 92.3 + 20720 | 50.4 107 | ge. «20702 of Time Steps | us | 4 | soe | ns 1003 0m 1760 | of Function Evaluations | 775 | 2036 | 1383 | 1632 2m na wn of Jacobian Evaluations | 70 |» |» |» 124 14 163 inal Step-Size (sec) 2.26 = 10° 1.04 « 10° | mm 51 | 5» “is un SOLUTION VIA DIPSUB EC EE — se ——— oo ____———— imum Global Fizor 7.47 = Td 5.07 » 1077 | 1. 10° 13.0 = 2 19.2 = 10 AAR 1074 | e0.0 » Md Tr Zl F-Pt Fl Jl I of Punction Evaluations 2687 1226 1592 14 ”iy “20 320) of Jacobian Evaluations » » “ |e © © “ inal Step-Size (sec) 1.63 = 10? | 1.97 = 10? | 1.53 = 20° “a wm me an vu rime (min « 107%) or | 207 73 | 2 3 0 "0 B Nn | | 1 ME EERE: i LTS ¥ ¢ Cs : ximum Clooal Error a. ad | 5.03 = 1077 5.82 » Tad 1.0 = ad wa 9122.8 wise. wi? of Time Steps 2s se | so8 ml un 1554 2038 of Function Evaluations | 743 10 | 1m 1750 | 2300 3266 an of Jscavien Evaluations | 48 5 @ Po) I ss 54 . inal Step-Size (sec) { ne 1.98 = 20° | 1.35 «207 | 60s a2 nr 258 ’ IE SE E— { ¥ Table 5.5: Problem 3 i

ee —————————————————————————— H Sree Ti se

i goa {

} PROPERTIES OF TUE NUMERICAL SOLUTION i at eres somnd [ae 1? [ 300 10? ot wht wi? llaximm Global Lrror 9.03 = 107 | 1.0 = 10”? na» wr] wo. Td 20.8 = 207% fae . fe “ese 2

Tree SATION rr PO coven I—.

[==] | seus 0 | oo -oy |

oy = 10 wy = 0, tan(ss?) = 14.3 | : and

oe 0.1 SROPERTILS OF THE NUMERICAL SOLUTION R——— —— TT IRIS—— 1 Lrror Bound 107 107? | 10° | 20° | 10730 wi wi? l _— . SOLUTION VIA DIFINT IRENE ee ——————————————————————————

of Time Steps | ae | ane 575 523 386 193 1084 | of Pusction valuations | 433 [ on | 1220 | 222s | 1295 108s 57 of Jscouian Evaluations | 36 “ w“ i | 107 160 inal Step-Size (sec) {aes ans | 1.07 520° | am 1.2210" [ase oe u Time (min + 107%) 3) 100 176 | 168 Ith | 260 4

s { SOLUTION VIA DIFSUR : of Time Stepr sz | 208 | 052 162 2610 2700 2020 i { of Function valuations 2540 | es | 2269 200% 264 7303 ”»n podpramerspmrsmerrnll [i | a | 20 ” » x » inal Step-size (sec) mn 0.3 | e.5 0. 22 ms %.3 Pi cise (min = 107%) wm ” mn i 7s nz os | / ¥ SE BUS— J EE ts SOLUTION VIA DIFBOF . ESSE ee — EE memes cers . uximm Global Error URRY DRY PRY DR ERE] DERE a EL of Time Steps | 1m | 205m | 300 % i 10.906 eo bof Function Evaluations | 400 | 52,70 ne 1973 1524 31,969 10,597 of Jaccbisn Evaluations | 27 | 1 | »s » 1» n inal Step-5ize (sec) 1.2% 30% | 0s w2073 | 12020” | 00 1.03 + 3 v.58 = 207 | 0.2% , >

FR EE E———— “Maximum stop-size permitted. p

Table 5.7: Problem 5 4

i 119 No CPU time is listed for the numerical solutions obtained using DIFBDF as this computer program was used only in a I= diagnostic mode making any time comparisons meaningless. The conditions for which the numerical solutions were ob- tained are contained in Table 5.8. Data tabulated in Tables 5.3 through 5.7 are for the mesh point solutions obtained up to and including the first mesh point beyond the final time specified in Table 5.8.% Problem 1, whose solution is a decaying exponential, [| was chosen in order to obtain comparative data for the three algorithms without any interfering factors such as stiffness \ | and non-linearities. As can be seen from Table 5.3, the numerical solutions as obtained from DIFSUB and DIFBDF are 1 comparable in accuracy, number of function evaluations and / | number of Jacobian evaluations. The numerical solution as obtained using DIFJMT, however, 1s not as accurate, for a H given requested local error. Moreoever, for local error bounds less than 10710, DIFJMT took more time steps and correspondingly more function evaluations. For higher error a | bounds, the opposite is true. In all cases DIFJMT required ’ more Jacobian evaluations than both DIFSUB and DIFBDF, a point to be discussed later.

Po 2 *Por DIFJMT and DIFBDF the final mesh point included in the summaries may be 2 or 3 points beyond the final time, since, i for each call, the solution at 3 or 4 mesh points is com- puted.

i 120

© t @| x ) [ [ Ll ol 'ol © wv — - 4 i ”™

I=} al wo | ~ * | 1 ) ool m | a ol © . of 5 p= =| BE BE °

. i =

= | pf i |e E of ~ Z = EERE z x ol of © Po = = I °

/ ; g

ol» -

x| © A % Nl - 4

/ : |

— — < v| © y- ol © &

| © S ol © © sl » rr] Pe. E . = P) E 7 " 1 u 3 w of © ” : al A A | 3 | ol o © y : ! nl on 8] 4 wil | A © < wl wl wl © = : 1 } 1 0 al aol [0] Q o » +» FY Q h wl wn un =] ’ gl E [2 Oo © g FE # E Bl ~ ” i Bt | I ol A 5 % = ul sl df 8

] 121 Problem 2, <uggested by Krogh [KRO70], is both non- $ | linear and stiff. The Jacobian matrix J is given as -1000 + 2x, 0 0 0 0 -800 + 2x. 0 0 J =U" x - x U, (5.43) 0 0 10 + 2x, 0 | 0 bh} 0 -0.001 + 2x4) where U and x are as in Table 5.4. The Jacobian matrix at t=20is -1002 0 0 0 | \ 0 -802 0 o | ) J, =U x x u, (5.44) / 0 0 g 0 | 0 0 0 po

and as t + @ is ol 0 0 0 | [) -800 0 HE J, =U" «x |x uv. (5.45) ' 0 [¢] -12 0 ] [1] 0 -0.001 4 | ‘ . E y { Again, the numerical solutions obtained using DIFSUB and ’ { DIFBDF are comparable* with DIFSUB requiring more function evaluations but fewer Jacobian evaluations. Using DIFJMT ’

*Por an unexplainable reason DIFSUB aborted when a local error bound of 10711 was requested.

the number of time steps required was comparable but the num- ber of function and Jacobian evaluations were significantly ( higher. The maximum global errors for DIFJMT are worse by

p approximately a factor of three for all test cases run. As the eigenvalues of the Jacobian matrix are all real it was expected that the numerical results would be comparable for all three algorithms. Yote: The extra function and Jacobian { evaluations used by DIFJMT will be discussed later. t Problem 3, also suggested by Krogh [KRO70], is similar to Problem 2 except that the Jacobian matrix has a pair of complex eigenvalues. The Jacobian matrix, J, is given as 10 + xy 10 - x, 0 0 ~ -10 + X, 10 + x 0 0 J =U" x x U, (5.46) 0 0 -1000 + 2x, 0 0 0 0 -0.001 + 2x,

where U and x are as in Table 5.5. As t increases from 0 the Jacobian matrix goes from J to J, where

| 8 10 0 0 ' -10 8 0 0 3, = u' x |x u, (5.47) 0 0 -1002 0 | ’ 0 0 0 =-2.001) LF

i 123

p H -10 -10 0 0 10 10 0 0 | 3 =u : |x v. (5.48) [1] 0 -1000 0 | 0 0 0 ~0.001) As DIFJMT is more ideally suited to the numerical integration of differential equations with complex eigenvalues than DIFSUB and DIFBDF* "better" results obtaine d with DIFJMT were expected. Again, the maximum global errors contained i I in Table 5.5 are higher for DIFJMT than for the other two bh i algorithms. Moreoever, the average number of Jacobian and function evaluations per time step are higher for DIFJMT. ) Problem 4 was specifically chosen to exhibit the im- 4

i 124 stable for all step-sizes. However, neither the order 5 nor } BN the order 6 backward differentiation formulas are asympto- t tically stable for all h for this problem. For a given step h, the point in the -plane corresponding to the eigenvalue =o + Juy is, with ¢ = Re{)} and = Im{A}, 0 + Ju = -o,h + ju;h. The locus of points for bh ¢ (0, =) is the oper ray in the 1 second quadrant of the ~omplex plane defined by w= - | co. (5.49)

This ray together with the Lambda Loci for the order 5 and i order 6 BDF and the Lambda Loci for the new methods of orders I 5, 6, and 7, is shown in Figure 5.2. As the step-size in- § creases the operating point in the A-plane "moves out" along N the ray defined by (5.49). If, for a given method and a given step-size h, o,h + ju,h is in an unstable region of the )-plane—that is, the local error would increase due to numerical instability—the algorithm should detect the con- dition and either decrease the step, descrease the order, i i | or both to satisfy the error criterion. From Figure $.2, ’ for the order 6 BDF the critical step-size for Problem 4 is approximately 0.04 whereas for the order 5 BDF it is aporox- , > imately 0.1. Table 5.6 indicates that for local error bounds ‘ greater than 10-7 the final step-size used by DIFSUB is ap- 3 proximately 0.04. Though not shown in Table 5.6, for these

2E |

. : .

RAY DEFINED / BY (5.49) pa

4 . = /

i Ny [ ie. ~~ | —— Sy ' — ee ea yi " | m———— — gan Nig 1 ” [ i / } [{ r ORDER © { ORDER § t er | \L r

i = f

Ni. § : I oN / a) Backward Differentiation Formulas

” / / [

RAY DEFINED . 8Y (549) » —— eal “N i . jig "

—_—t ee eee PRE . I 7 0 (TV | pr C

_.. <a b { Ni \ f

ORDER 7 _ / i Fa ORDER 6 —__[ ORDER 5 2.f : ¢ \ + ' h Le ( b) New Methods { Figure 5.2: Lambda Loci and Ray pefined by (5.49) for Problem 4

126 local error bounds an order 6 method was in use by DIFSUB at the final time—that is, DIFSUB was *locked in", properly detecting that an increase in h would cause the local error to increase beyond acceptable limits It did not, however, recognize that the step could be increased with a decrease in order. Similarly, DIFBDF got "jocked in" with an order 5 method for local error bounds specified at 1071 ana 10712, + ¢ However, in all of the test cases run, DIFJMT was able to successfully increase h after the transient effects due to the eigenvalues at -o, + ju, decayed. Note: Though the new order 7 method suffers the same disadvantage as the order 5 and 6 backward differentiation formulas, DIFJMT was not “caught”. This is partially due to the order selection logic integral to DIFJMI. It is this same logic that most likely [ accounts for DIFBDF not getting "caught" with an order 6 method, though it did fail to decrease the order below 5. [ As expected, the numerical solution for Problem 4 was obtained in significantly fewer time steps using DIFJMT than with DIFSUB. An improvement over DIFBDF is also apparent, 2 though it is not as dramatic. As in other problems, the i maximum global errors were not as low for solutions obtained by DIFJMT, when compared to the other two algorithms. Again, ,o the average number of function and Jacobian evaluations per J— *DIFBDF also got temporarily "locked in" when a local error bound of 10~2 was specified but was able to eventually in- crease h by decreasing the order.

= ~ |

1 wm time point were higher with Dreamer. A ee RE So mien tn, the Samian asin 3 io 1 pred t oy wy © ’ i bv [] “hese U, oy, uj, and 0, are as in Table 5.7. Note: o, i BRE iat oo 2 EER Lay voters ERT See AEE We ne tn pe art PRE ti I with DIFSUB or DIFEDF. The improvement, though, is not as taal ih I The following conclusions can be dravn from the re- : bia Ak Aaa Ei pA. sesteds I 3. Bac mea-ekis pushioms ant stir pecions saaratn the eigenvalues of the Jacobian matrix are real, numerical Sie ail dy wien ifhctsion T3gntoa 8 com pin pede RI LT Ee 4 ‘tive real parts DIFJNT requires fewer solution points. BEE i ctor ont mectan sc uations per solution point than do either DIPSUB or DIFEDF. | |

  1. Numerical solutions obtained using DIFJMT have

a somewhat higher global error than do numerical solutions 3 obtained with either DIFSUB or DIFBDF. 4. The means used for storing information has an effect on the properties of the numerical solution. It is felt that the higher number of function evalua- tions and Jacobian evaluations per mesh point required by DIFJMT to numerically integrate a differential equation are

p due to the predictor formulas used. Though the predictor does not affect stability—as the corrector is iterated to convergence—it does affect the convergence properties of ) | the corrector. The modified Newton-Rhapson iteration des- I eribed in Section 5.2.1 for solving the implicit equations is dependent on a predicted value that is "close" to the true solution. Using the same predictor for all 1 equations in (5.1) may not be the optimal strategy. A study of diag- nostic output obtained from DIFIMT revealed that the correc- tor failed to converge significantly more often than for DIFSUB and DIFBDF* resulting in more Jacobian evaluations 4 and correspondingly more function evaluations. In turn, .

this detracted from the results tabulated for DIFJMT.

FT — Pp *Both DIFSUB and DIFBDF use the same scheme as DIFJMT for i solving the implicit equations.

Tote: The step-size is decreased by a factor of 4 if the corrector fails to converge at the same mesh point on two successive attempts. Hence the number of solution points required is also increased if this condition occurs.

H 129 The author has not been able to explain the higher maxirum global errors observed when DIFJMT was used as com- pared to the numerical results obtained using either DIFSUB or DIFBDF. On the average, the maximum global errors ob- served were worse by a factor of approximately 3% for DIFJMT. That numerical results obtained from DIFSUB and | DIFBDF—two algorithms based on the same methods, namely the backward differentiation formulas—are not comparable for all test cases implies that the means for storing infor- mation affects the properties of the numerical solution. Note: DIFSUB uses the Nordsieck vector whereas DIFBDF stores information in the form of backward differences.

CHAPTER VI ” CONCLUSIONS AND FUTURE WORK New techniques were obtained by which the stability properties of composite linear multistep methods may be accurately characterized. The relationship of the Lambda and Zeta Loci to those stability properties was established. A computer program for plotting both the Lambda and Zeta Loci was presented and described. Using an interactive version of this program new and efficient cyclic composite linear multistep methods of orders three through seven were obtained These new methods form the basis for a new vari- able step, variable order integration algorithm described. This new algorithm, designed for the numerical inte- gration of stiff systems of first order ordinary differen- tial equations, exhibits better stability properties than available algorithms based on the backward differentiation formulas. It is ideally suited for the numerical solution of stiff systems wherein the Jacobian matrix has complex 4 eigenvalues. Numerical results to support this claim were Fe j presented, including a comparison of the new algorithm with Gear's popular computer program DIFSUB and with a third algorithm based on the backward differentiation formulas. 4 The new algorithm, however, exhibits a slightly higher maximum global error than the two computer programs

used as a baseline. Moreover, the predictor formulas used in the new algorithm affect the convergence properties of ’ the component formulas and detract from the results obtained. Purther research is required to obtain "compatible® predic- tor formulas. As the new algorithm presented herein must interpolate new mesh points each time the integration step is changed-- a cost that goes up linearly with the number of differential equations being solved and the order method being used--ef-

fort should be expended to reformulate these new methods in terms of modified backward differences (equivalently,

divided differences) wherein interpolation is not required

each time the integration step-size is changed--an approach used by Brayton, Gustavson, and Hachtel [BRA72a)] and also

by Krogh [KRO72].

Stability properties of a composite linear multistep

method were considered only in terms of a fixed step-size.

As this is rarely the situation research is required to investigate the stability characteristics when the step-

4 size varies--see, for example, the recent paper by Brayton a and Conley [BRA72b].

The new methods presented were obtained solely on the basis of optimizing the Widlund wedge angle associated ’ with the method. It is felt that attention should also be focused on simultaneously minimizing the local discretization

132

error. In addition, higher order stiffly stable methods j should be derived. The test for stiff stability of a composite linear p multistep method presented herein is entirely graphical. For linear multistep method--that is, multistep methods with 3 L = l--Ngrsett [NOR69] has presented an analytic condition ( » that can be tested to determine if a method is A(a)-stable for some ae(0,n/2). Rubin and Bickart [RUB72] reported an analytic condition for testing whether or not a composite linear multistep method is A-stable. An expanded treatment of that test and of related topics can be found in [RUB73]. Research is required to obtain analytic conditions by which to test for stiff stability.

APPENDIX A

Lib SS Sinan . *On the XDS SIGMA-5 computer double precision corresponds to 56 binary bits of mantissa, or approximately 16.8 decimal digits of accuracy. z =
rR 1

i 134 If it is desired to implement corollary 3.12 to test H if a method is stiffly stable for a particular value of Y of Definition 3.6, card number 77 of the MAIN program seg- t ment (deck identification "MAIN") should be changed to rcad i »yAR = (y, 0.)". Mo other changes are re yuired.

2353353535358 S RAR SE RESETS LITRES ARARECECALIRALI INTIS srprrezssveryrressrry zy EsEErEEIRILEILLL LLL IE EERE 2RRRR00 HE EEE EEE EEE EEE CEE EE EERE EE EEE EEE EEE EE LEE LEE EE . a BS : Fee

g os SE° : g . IE Ht { £ 2 A 8 E ¢ ll Pl z 3 . 3 § #231 TiS 3 £ H st E:T, It v Fore : ¢ > . ow 2%P.C%. a ¢ § =f *5 @ £ Sh.dglizt : H : 85 sit 3 § 8783225528 { H H Ee a3 BB ) ; Jrisosts. E i SEER 3B: Hck HEA cf BiB Bored EULER Log §, pil oe BE ELLE RS ET EI ER Eat NE Bt 3 5 3% gicitipds seine Sedpritizaia pos Sf BRifioRd, $ Re DU ES ER ie i LE Sets I 4%: Cliffe Soe Tata Tita Co fs fRiEITie SENET oe gi5e BEse.%3 ve - < -_ - 8 3® EER 2X22 ov ov mus caveversecoararsnue gs Ta IRERRAARARARTASIILININIILINESLID z2zzzzrzErrzIrIrazrEaErzszzE 22 ERE AE EIR IRR ER LLI000 FEE EH EEE EE HE EE EEE CHEE EOE CE EL EAE EAE Sh 2 Ff: : si 5 alas 2 H FI J : EE LL FE oT £ ETRE 3 2 yy Z FMisid 1: fRicils ie HS S 5273.75 : = ES SSz232y p32 PoE if fiitig 3p b Bgidnd ii 12iig *§ 2 s Egotatii FHET TREE" : 2 f 0; Blatt iN Canin £5... 8 Beda iE il.pn-idd SEE EZ Paps FEEEEE ool 32% 2 Eh: as $ 8 RiavacPats: tL $ oxB895E5282" § °O RR TR Rl KA PCR I cr

srsescassrcisesasasiiRtisifefEbiRisiIIRIIIzREbREbRRERRENNE pryzarssryrzare srry ERZEZr EEL ii ln srrrarrIEIEIZIIZ ERLEREL CE EEE EEE Cr EEE EEE EE EE EE Eh

137 F333 zzrzzzrrzzr [FEFTTELT P— i § $ 3 ry Et § | } ital |

= ER ae : Bz 4 : 2.514 : §2 4 HRT BE ol | Post H ty Pk

PoAEst 3 235%LS

§ 5{% £28, { .3 288 "si s 28 | 82 ER: —- : saaspimiiiiEcEfdnfiplamEEnMMEOIILIILENITIINIMN SerzsrrrrrzzzzzzrzzrazirErrzzIiiIIriziiiiiiiiiiiiiiiiil EE EE EEE EEE EE EEL ChE EL COE CE DLL ELLE ELE ” 3 £: H tit : 2 $ 3 iii i ij i = iE : gi ig Ha ; Te H SE EN: : £ fed $ £3 EY F Sa ded : ii HE 7 1 d§ 5 iii iF. te} i 3 ES $8: S11 % §2 ii ig = c 3 5 s+ 3 1 HE 3 SE: 3 Ser p oa | && Ss ££ $.31:¢ is gf: 2 i% i t= ER 3 32 $EIFZ $f. 3a: &::23% . E +8; B ifi:° HEIR EEE fer f 285 EB iii: $I: Fesf iad. iT i 1 §5, 8 i1:iES 133 § op ignnioHe ©, £ 285 TF 38; it $8 7 3e3EF.E 1 ER% $F FE OP irigt o $e fF 353 ies ih $27 sar pie iP 2: 3 ie. 35 iiafvii iss ’ SotBey Egg SSF 3 3 OYpe 353 #88 Ff ili 3 AES H [T3Y zat TE t 3718 <3 ed & & 7.215% fogplaerfaise i PEE EI § iorsb3E ¢ 1 - 28552 E50: £280 ¢ L Sees = 3% d 2 %LT0%e : 5 =e 3. 8 H =. . s 1 JE J oS gi%e 53% f...7 533 > Sesssccsecse sh 2 Zeeen st 2: 2 2 H g § gk | ww © we” wusuuuee JUTE T d |

3 i y

4 b>

233353352 s33%28s y oN SLZRRZERRRES ——————— 222855021 20%% PAPE E —— H eee SE5ssEEsEss ye : gi PEATE sa : _ gil PATTIE rr o£ § Ii seer 7 Tr massszEszezssiy , i a teeseieeeeese " 1 feet i gf I2iet : Lg

—— : NO— :

i 139 erspssacsrazIsssssssIsResaEtaRRREL ETE IRE ERRRRLERARERINE EBS ata EEE EE ERC a SR RE ees ) BT TL LETT I A A LL i A EE rr

B ; : £ a 3 € = is TN : : : : 2 Eo fs $5 P

H se } £ z er B I (o%e, 2s : Sor 22 £ Zee FH s : x 22% gi.is rx} F z = Fy ga- B2°28 <5 14 i 7% §6i Lui. oz B #58, Hg i: sss Bpin oe MeL 325288 vd FE o gi: REC se EH 2 SEE ys £ = 258s 5.0.8: $.% 2 +L fi:2c. o. as gE er dE os. SEFC C0 gisds 885 Torres . EUS SI,2I% Solois §r22S52 2. me _ Eden EoCiPidiriEa; nc Bg ziiiiEyiiiiiiass caiisigiiEiiif Figg RTC BH Pa DE Err a NTI FEEEs gfhnnadt Ht aul FL JP : Ee rT " 8 L28 = | ET Lis wszsr oo uf 35.0 2 ¢

’ \

140 rib aEREsRARIASEsR RRs ERRARRRE RRR ARAL PT 0. .

3 A A Hr Er EE Ee 4 : i $ 1 i § % 3 H : H 2 . : gx : § I i [4 = > . 13 : gz is i 8” . = 8 f+] 3 gs - : 1 i iu Ee. 2 3 g i HE EB. = "ii : Sts flan Bs Ed é =n ge o 222° Se oF B g : Ss g37. g Se ptEh 33-oof ¢ Bx 1 teilaesd TUE BAMHI. GG 8% $233 T 3= [Rgls feAlf fase 3 29300 2 Bo 5 58° 28 siiEfe § Bigasis AEE 59 533 § Either oi ET le TL LR LE oF Booo. BRBEiofesiS: BED ZS.2T ZEEE ZT FE S...B RRs 0, B23 tiie i TERRE Lo oeeiisusiiis Cac 2000 eeeeg® EE 3 3333 fiistirassete pe seer tues BEER = eo: we ik a gz #222 2 §3 = : cas zsrasiasis EEiEaE ant EcEsEaRaRatREREs RT EREMSRIREMAS * sasiziziiiasiiifEeesfifiiizaaaninaiimmEInsnaaIEIesRsaatn ay SetiisssessirrrsrTITITriiistsssriistiiiiifiiiiiiiiiiiiy

s i

: e ) 3 4 3 s 1 B 3 H 3 i E35 i { § ° 5é o H : . . N : ‘ E 8: § - ks & te £2: £5 i 3.8 f= £ . a2 8 3 35. 8 £ aeE x oe =a ®o Ty s20%s 3 - aot 8 H 2 apfad 21 gE EE; GC ga Eg dees 3 : { "22 wks fof IF EE 238. 3 $2 Epica f } nighedd JE Regt 80 BILE: WE E37 338 sof: 1 sbEEc.- 33 s SeasBEEE Ct IS i Bs TEs 23t B50 LO SSHEs Viz so FITSREIICE: Yrs 2. $.7 33°cacy fio BeR3TEoian 20.00, PeiZRERRECIRaLETNTS, 2 &F chs insl s.oge = BT SL et 1- H BE Te i 151 3 E3gLressstooaznisiotEE T 7£- 2.iged2 8.08. Sef, dTu¥egust, opt Tio tanalouinanaly o B31 ise Foeiietsts Sufseiizs Fs EERE iz s f & sg g 2 g 2 =: :8 2 2 2 g E38 E & LS g 2 : $2 8:3: § LE

£35505333598KERRAR SRERUSACITILILGEILINLSILESILAALRE ZI Td A A ii H i r 2 - H i i: : Ea z : i i H i i : : > i bs = - » 4 13 1 H = t §.. ¢ . £ E 3 H 33 5 3 i y P55 . wf FP OI Lj Pr 1B ff. 1 3 i if Efe Bg 2 = £5 fa juz ic Msi 2 3 3 § Libsged igi e aidl: 8 FB 33 OB RRA fail : st =°2 ° yd a 8 «"8e7 2% = tT PSL Lt £EY £38 2 8g ; Figs of 8 : gy | “ES s BBs s IF 3 PIgu3T CT HR : #5% “33 ¢ ? s8; ¢ IHR. GL HT Bet a Ex, 23s _& H gi S «EYE: © : Te o 25° 3 33g 8f.sf. ou Z-ilg i S22E7 5 7 1 Nf 3 —od 39 PA FL Ta Te fava & =" 3® 522 3 s %E 2 »° Be 8E22%27..%. 3.2. 82:570.5 E85 $3780 15y EY JH8. 822," o80p 35. SES 3. F To plee zg iE zy [tT Teel SEfoTiEIE 2257 o5f% * 7] CoS t2% 2° *°

Er ES SehpeLsiEt eu FE “Fie IE S985 ~el 22 27 2 149 $35 SS — I. | ——

I mamemenssosazsrsn seers Lit ItRARRATRRRRITIS ITINITEAIRIIIT 3 NOs EE TB HTH

HE & 2 HE 5. F iE SB cs Rog : E11 25, 5 3 %8 tt ¥ g { ii Boos: p 9 B23 £ = $8; ¢ Sir Cop 3 i § : © $8 & 2acls- FSF $C = { $.F 8B E"%ut = $ 2s NT 4 Ss . $y: . 33585. 8 5, 3 £ Ei ps H i: § Frres: TL. CE POT F 2 $80 £ Chpoti Bs C. ¥Ei- i od £ 38 Tob 3 piiBL: tE 5: Laon i gg .¥ iE op ERELCI REED sips 2 oi = y $3: F oimROf gp 502 § 3 < i. jii , } p58 © SEE, 50 28 £85. 3 so = x 33 iE : . thf tire Eiy i 38 GE Ct. g $2:0% LEEecpyitiabEi os: 33 ge* 53 a2 § ied GEE aasiiiogl oi Bin eg -t : rir Eisfsays? Ly Le H ho] « os 23 3538s LZ REE HT RE fis 3 cif £3 PREECC os ¢ § Pmsmediii (Blige EGIL I oo Hi ¢ pp digg cE DRIES EGE NES BoaS S§ gEiigeem Lod ET FEY P08 232 sf 33 ERC

» 3-5 yee SS. 735 : 9T¢ TFS 33 S22 ivy Tal 4 PREzf fvs 22 33 Looks { GUL SUL UP EPL BL LLL UL BUDO BUDD BLE v “ wwe ov

carerersenosrsrane tsa LITIRALAREEIRARARITIIININIIEIIC H 3 H

fF 3 4 1if

7s T is Ha . ss £5 & 38; 238 5:7 ° pe g 225 = FIC 4 ig = 238 ;¢ Ci¢ ¢ HH eo ru ru £3. . - $= § 238 33 3:1" ¢ i fe Gp Achim fpi OV : : iBg% § § peitizi. Beli i ii ¢ $888 TT FUSE eED 2 : H - $PTT Sf GREiicERi: Lif fz if 3 YL 3 @ Elio EhEEiel $02: 00 ZF : 3 eo. sp BepEniie Ri: #.f iI i 8 iE SScaaciyaess Buife 32 TY = POC OS fifi igameiteptd RIG S Ef 2 cf 2s § FL BRL C.aRus ans $52 se = a 2 3 § SE 3 eo SEEUTCRK et 25 i - #5 ¥F L y $e 277555 LEE REY 427 Pefse B2fis SF BSERIL 5.7 %,0 Bit Gif Flin Flloze 1 iE §* °* poi fs 23° oon Too ovErhl 2 . s eo pis 4 sesecenaretE Rana RAns REI sIEItIRaRR2E200023322322222S

EE ELE EERE RRR LEER LEER S

H £3 i. -3 | i i oe Es -% | : HE a mee EEE EEE i 3 SR sts tLe i § © gE tg CENET

) Y

. 143

23555238580 8E2 0 RRRRRCaT I 2 ILSREIRSIIRIRNALIRRISEREZEE Bees nr a a TT CEERREETLa Reet aieaiciatiieiibenteaetiinitatiotiiiniste 8 . 'S. 9 E- 4 H ge ®t : i 2 — ¥ : 5 , E 0 3 & ts : H es s 3 - H : <x * § 2; - 3 g ¥ . E 2 :. Fg £ i Es ¥ z 1: <5 3 . tial H HE = . 5 Se = § ig } HIE 2: cE g $s ites 8s Et 3-2 LEE RR 3 3 » = $f Si.87 ; gE EEE SEL Tt 5: 1 &f gd dhl s -3% FEY if Zi FE = - ZeuS 5 5 222% 2.2 3 $3 0%: ESF cf zo.8k 31g: | 3 PR 1 3 15 i ¥:29%,2 227.2752 2 To #0555, 3 2.780 £ ES.aEPTESL ER Ye giiid oops ® SpEESHE ° Tonal SoiHIMNEC i 3. ES EREIR 27: EEEELaER a sf fii fener de Ue RRR Jl wd 3. ss i car erenserosrane ane LAT LERiLRRRRARERARLIITITININANIINIT | JE restart TT 111 1 1 I aH EH RH HE TI HEE { Beppe Eee iiidsE Eta zis tert abRRisetitinstistitnontiioniis H 3 E : z : . 3 H 3 gi & A H $ 81 i | I H : Fo & $s : gs B 2: : tf ¢ £ EL EE gs 7} bk Rok ARE Sk ER : . E tS:54¢8 = ¥ 3 gst EF ® iff Ne tL. I; cs ed BE res ts LET 9 4 H $82.3 23 3p, 2 pt. FEE § :7fci fF eit : E5°af [43 PORES osm: As oH 4 TET ONE TEE Het. BT : ‘ £ ia og oBRPl pEEEEE I g575a EZ 4 ® . "3 ef geniz, LEI ZT 3 RAE RE 3 s s PUI El ERLE: f.T38 Ck 2 : HH HE PIT Eo : vad. 134 s $ 1 © iS SERSEOCMLEIST § NEN is : sé = s 3 SFT E_CSERESTIEMLE : BREuC. 5 ZED 3 E & E32 Foe TeTEEIID MEI MR TEL OY : sf 5 N = Ft 1 : - eds * 23FTeS © La - g <8 Fa E :C¥ so..i. 2d : SE: ! : £2:7 y . § Est ® Ei 8g iffpey 2 yz i s&s £ £5: 3:8 SAE Bie 1 $55.2 Fr :oiY By 3 BIFIfi cf Givu oo ® OEZELE ti: fess; 3 oEitli Ra LUPE URS PY UU PELE BL

sgzszaiizis EEEEErETIIIG i ETHEY {1171 4 ; H ¥ i

had baa 3 s P=S% H i Pak | ! 5 obiiaiat F.uliisiict i $Me Tiled? |} 2 H : 1 srsespsarssatssszamEaIiiiiiIINIILEEZIIIILAIEIIIIEIILOLLD ; I Err TR HR . fiiiiiiirizisciicieieisiredsebedst it ibERIARIARERILLGNL I EpEEREEEIEEEERERRRRERERERRRRRRRRRRRAREERIARIARILILRIIEL

2 4 . bd - : H g 3 § g A & : i £ g : z ‘ : 3 = : : £ i: «¥F © EB 3 ce “i s » » 5 ek & 3 3c : 2: 5 - 2 i p32 EE BE » 2 gz 8 3 - : ¢ Fest A Sar 23% o-8¢ #. rd RE $1 go Ef gr E3i3 CE TE ¥ efs = Ek] 52s o% Ee3t-2n . . ss == iz i EERE 2 rr *3 33 ’ EEE 2 go oueyEzi YS TIUUEEE : gC EE if po mRba bh OECARi Taal hag | EEETpned B SI grlcBEtt UU peRstIRiifEil “Ee wn 1 sesleslt Feeflarfite: HS Ted TRC Se EE EE & gs ser 2 3 3 s 2 8 vou vou pe wu www |

2B |

145 2595853055 SERRE AL ARNLESALLLT LEETIESLILSLELEELIEEELEE IEE . SEESSLEELELLABELLLEROCLLLLY

SEERRREIEENRERNR SREY IIY SREY

H § Se t tg i s : . .: i pt s :

z2 . H

2% = H

1H $s ;

ez Ee oe

3? s2c® 3%}

$8 see£ 8; = 1

£8 2% FC 77 i

£2 Sis Et i

: 38 E8228 22 i

: _§% 5.27; £2 i ests EE lB joker Bi: ni gions mnilabiy v Tged RR HE RH Ta

rast PlesdiiTeeme fetatiteenis zs 2 3s 23% §

§ careremsvcoansninse oa 0RRIL RFR RAREETIIITINIISANEINIS Jpeg TTT re a RE A HT LCE tt LE CLE

i: yy = H 2 EC Py af fF. 018 13 1 Wh; i igs 5, Bb ELE S851 3 $f: BE i Lad

E 3 Hu Foi Sf : fs 8% a6 pli pEit I igpicioszB3 ski EA Th FRERSES CEE tpi ISTHE 9

3 Egibpsmgfil HG dEInE oe : es | PYOP gREEeSE eg} pf 1 ZI. 2B os 5 .

: § ig EiEiiec CEG, dei. fill Bit : v g 34 BEETEL ZEETSsr io43:5p oo; 3 3 £3 25 SpEtXil §3.03328 § SFE. Bes $ Kh : Z¥T RINt.TE §7¢507 § suzaxeg SaE 7 < < 134 Poin Nba Lonrifaldn cy qo : i : I¥S . 20.2 ’ $ saste $s : 2 =2

$5 o Bao : £22 $s : =: T0 A LL ses 13 Lye. l $: & #337 HET B 3. (¥.fee Fetes

APPENDIX B ’ .

4 COMPUTER PROGRAM LISTING FOR INTERACTIVE VERSION OF SUBROUTINE ACCEPT i This appendix contains a brief description and a computer program listing of the interactive version of sub- 2 t routine ACCEPT used in deriving the new methods contained $ in Chapter IV. Though referred to as a single program the inter- [4 active version of ACCEPT consists of several different units. A block diagram showing the six main segments is i given in Figure B.1. Each segment is a separate FORTRAN subroutine. NOTE: Though several of the subroutines shown in Figure B.l, and indeed the entire package, could have I been combined into a single subroutine they were purposely written as independent units. This was done in order to provide the overall program with greater flexibility. Other versions of the plotting program use some of these units. ' As an example, another form of the plotting program uses a different ACCEPT subroutine, one which reads a method from cards and then calls the same FILL and DETERMIN sub- \ routines shown in Figure B.l to compute the characteristic polynomial. The computer programs listed on the following pages | were written in XDS EXTENDED FORTRAN-IV for execution on an i] |

i $8

® 0 & ~ . Lodi i E23< us 8 = we va ale # @ 1 | EI

TE 5 La [3]

El

i : i 8 : > s ad fs s 34 £3 8 1 Bal | | = " > | i £1 5 u 2 | a =< E goga 3735 E ¢ w 1 i = ols ga CHE st nwo E IA a0 { [ 84m gE 4 FE | 5 = Ct 2 s gh" | E 2 < 588383 $ ¢ Ee 6 5% £ EB - | = 5 Ew sd xn 28g 4 3 Bl 3 0 LRT Z Hy = ° & glssg"”® i © © a 09 - 3 ol a aA Qu" 3 AR FET HE © $e w Pm § bi [4 3 - °2 a 3 z MD HTB $e I: S998 88c3 | TYE [FE $$ = 0 | ° © 5 n - |

oY i { 80» 2 - ” py | S Blo 5 § 8 25 3 $3 a ! oem 8lg 3 RET $3 a : g » = i ® 38%. Fhe 2 . £ ¥ 8 0 $ EP £. 5 : g 3856 a x ™ - alo ee 9 » a8 9% Bl3gss | ERE ’

cg © 0®

a Ed

go 8

i 148 XDS SIGMA-5 digital computer. They were compiled using the i ADP—automatic double precision—option discussed in Appen- B dix A. Interactive communication was via a teletypewriter, the only device available for this purpose at the computer 1] facility used.

y { /

’ . 149 $34382358558E FRESRARERETASASILILNRISASIIILISARERIEEERZAES i Basra pt ane EOE REE EL EEE RARER EEE EE BEELEE BEE EE BREE BEE RUTR RIERA IRA RRRIRRRARR anna ’ A FRE , a. i g & 3 r 3 ay HS : s i i 6 z £ . A . ’ : Se 2 : - r . & = £ i Fos

¥ Eg A £5 $ £2. ° gd A 1 : $ z g £ EF ¥ - 1 5: : df dnb . H 31 3 " ‘ : tpn a A : df by Boot il Biss Nas § Poni Hie PMEBEINISSE: Go £ #3 +3 ¥Bp diBmzad waa fi.x 27 TREY EH LB _dre BEES Rime gst HER 2 Ls | 0 BE fPii ffm nas By HCG AN ud , ‘3 s= & EB LEh pommmR. BURY Too: § Bess & CIS TS ep Fat Ey EE pr graeai Bn:E bE f LaFISI Ei T25% 5; ob BE cer frivews £305 1308 5 JEL C0 7 7daiZ 2 Be gougit £ 300% = Lv 28F§s2 “2 =33°° i: BE scr sur ememarw c30e PIU BETLESE of Eze. TL VUE VA SY SUSU UU UU o

150 ERE asITIs ILE sRssI EARL EIIILsIRI RE 2 1 LE EH

i 3 f — 5 : E32 a = i ¢ E . H 2 A EY | 3 PT EE. SR 3" H sg 3%; 3 OT OF § PF Bi % 4B BO: ifs of af BB HEE IES i i=l Bea bs Sp. fs EE EE. Boi Vi 1 t5; 23 202% 2 23% FI O30 53 Zp sf Safi plfint DE Sn ld BFL BE SEE 2p cfsid Tip ET ITs Dif I Ti: Bp Bo Bed BOE Blam de BE Mz gd 8 FF Fes FINES Sl EIS LC” Lie 3 Bi £2 pe §7ifras £27 £10 fh {i RS Be a le CR £20 S23% Ti 3% Fp FE 35 PZ iF it I dp BEB ED ORDER BOE Ho Er TE Nw

§ 1

: ” e z § 2 Ss 3 § cE 3

151 S853552TS5 SLE AR LARTER LRILE HHT TTT fr Sos o008s0%59%% RESERREEEALE ey 333 3 18 % | DO FLL EH HHH SERERIRFTIEIINY » . 2 a Ex =f £ : E iE H €:- 3 z A ¢ > »¥ 225% . i s i gd * 3 c Tz M2 s : § g 28 i‘ B S Ee : a 5 22 * 53 $20 2 i . Ho br § Gd t sf & : i, £ - 1} # LE cs $1: a3; Tq. §2 . 0S $2 1% = 2 Bf 228: 2 s§: 5 : Fr 3 ip :oZR ofmmfis FE-PE TR 83 IEE fF Oo: 2 s 3&5 i Srila hr if an Eg [oid Hn FR 286 433520% ese 337 £5 "355, sy =f T gees 23 2% a:00 0s $7658 I323333F 35Teeiid BY diwiavd-i.o] "7: Banpeiis Toph FLEET A siTgieate, ind IE i Fh Sok NER RTH a Ad LAE J 2s 5s 8 LIFES P22 gs 22 ho SRE ig I.-"8 Fas wwe * Tw ao 2 2 32 3% 328 i JAR RA eS roe ERRLEBEALRIALEL + TATILAXRRRAREARARARLIITININILEINNIIAL . EAEIAEALKIARERARERILAMAMALBAERRIELS 1TITILEFRIIEAT SEeeriier FERALIASMAMMALAIERALEIRIZANERIAIEALALE sexs i t si mppEE EL EI RH = == $id ole w 9 Fs $F pct Te SEE o #5, « Et; EE Tat Fd E:f 34 1 gz: zk *poesle 12° = | FEE pREIINEcE 2 of Ez B Hires i: x wbh tf 2S : s SE - —~ 22'S Tol 1 ] cs vibes p55 22 $ Sf BE ak i: = i of mR adil 2s 2 2 Sef o:zel 1 A ghEssdEs 23. 3 +s ™M:1 St] g H

152

=ofsEarIiTEsRRafanafusITiTITIISSNAsiiiiiiiiiiiiic ERERRRERRiARKARARLARARBRLEREEAREERERDERADIEEETEIEES SP ELUPELEEELILLTEEEES ESSERE IEEIERIR IRIIINIIIIINNST

ii | s £2 2%3 2% g i = £3: €:

. cemerenrerouoor en t2gITIRALALARARATARRALT

22222222202 222222 HHH HUH HE I : : : { a PE. Po} : iE PE . z H -s 4 is Pu i g 2 : | ! sf I g . i F< 3 52 H TIE PN =n2zg od 8 $2 fs 2 Hyd : : $e Fs 3 dts: ff § 2 Lit Fas. 8 fo =x Tf, BE = iY Of if $7T..40 5% 3. FEO: Ug Paggasze’ 4-3 og af g2c ¥e ES Pooahfifed ders BEES REIT: S EREEAIIUEC ifhielrytiRzivy ti HH ¥ a» FIECITUERT SentesocoRil ts 1d 1 ALLTEL MLN a FP { ¥ Shiearadic of.o5hs Eres SLm E0ES, & Stipsipfsl s2.v fF x wiv CJ gs 23 ovo - . wmv evercessansE at jit : i HHH

A 3 y v H Hl . i | i y Pa e = : & x $s ek SE gas ’ saF es. Be » 253 § ep fis 3. AE PR Thier SEARLS pd I - 222 |

154

3353333835358 RRAR AEE 333333332323333233332 { gegEsesesEaiiiiitt

3 : s 5 . Ef

: - ” 8 * 3 fe i pe 28 2° £3 siiezy i FERS TE ERE 8 £38 Ie 38 FH 8°33 22° i 3 25% $3 23, Zoi. Eo hr Ls . Ss TS rT 23% 3-8l2%E2, i ex” “2

: se 22 I camemensorar mene TaD LR RRAREARRRRSLISSTTIIIEINNANID 43243344233333 33 22433333333 233332232222222333 332423332322 1 2d A BE RE ERIE IIR Tee

  • & it: £: [TE 22 AR Fe 3 : 2 | [41 1 $s 3p = : 33 . ’ i. 5; = > 3 H3.3 23° § 2 iz i: £ 33 ro J - BH » ec 8 $8 E38; deg diz 8 z: sts: ESC. % $2 23 } gs TREE EE s s22 3 z ¥ E83 2 882 LW 2 14 ed 2 PF. = sr ‘ AER EE THEE tT: 37a Pi : zr § rr. Si.3F 335 3 3 z3 5338.E° - ¥ or * oo3 SHEER 33. opt & 38 353i § 2 i 5.0 haley 3K BPE $8 BEE 3 fii i [ BelER Headly fo TEER HG
    Lo 2 BRYRLELEF SEFRVSS 2. gr; Enid HE 2 2 Hd 33°28 a ZENS a Tr is © eS 232 SF "5 2» £ 2 EEEEYP §asEea~ ‘t= 22. “~.%=:s IB 2 il 1

155 PTT IT I CL EE A LER SA A AA

FRR a EE § £255 LT rs 2 2857 B asrat 4 F§:3i13 a ' § fessk. E -4 Xr £ : TE i H Z Fad ¥ oz ® . ] Pag: = oC : Ef gE: om 13. WB SRE oo. 2 s : 385% zk 2 i Eee B25 LT 3.7 BE oer U5 Sia iE GL Fo 1 gE DosR °F a 22 gErzz £537 o3 3 BT 7 BR 7, oBzg, tt BE REFER Rep RT {013 s2 2 =; tlgge~t* Toe ae 22 kengz BF CptRia § pr Sepp 73 T5022% Te S834 § £8 58s Tl. 33873s Pigg $31 wEIE oop dogpile i Z 3% 38k SEE gos LE BE EE IE 2%:

156 PT re A rr TH HR HH HR a HE EH BEER RSEACAKEEBEAbEErEFEED BEEERESTRLT KER ELEBEERET K 4 ’ : ¥ §F $3. S 5; 50, : B§Eye z § = *® 3uxf H ° z = : 3 § 5s. : : = : = FE . e 5 § fii ts i Of EB. :oBlpmt sz . ai. Beli ID gd SiEEEE - & TBE-=e 92f: FIT «38 #° BEEifii : 2 : $°9ise ois Eebodiiwes aaiitii £ i, F $i335E3e L73e 2o3;-202529585% eens o£ i FoB.35c7 832-27-33a TURIN" BIR § > £ e33.2335 B2.P. CS SReeg3ciT 0. %e § 2 : so - Uf goegtapslg” “ped Ei ERCF is ii _f S:2"s2 8 s It 2333 B3iy =F _ Sgr, tis

i z sx 2 32% BB “es” 3 epznrmaniz tx H i 3 333 8% ow 3 smamzmasis | 3 ms EazigEnItsAs tan RET aE ssn EER0IT atin ItINIIIIANIY HT a Ea i BE RRR ER Raat es Bae BER ES £2: 2 % $88 H § if TAR . gE = E §5°.¢% 3 H Et UW - : : bog g* &2% fe § i E, ' i - Ba 357 "gst £ : = H 4 £2 £723. Sse uieted 2 .

157 SE 3p3LEsTEE

EEEEEEEEEEEE 4 PRaEEEEERREE

Eo: : } £

“Sa

pe H

§ §

i. ¢ 8 Syi-a, ¢

Ligetisiss 2 ERarirs, i £2%e

252 : es g83 £ 8 rassrssiiipiiiinIEeERRLEbRAEERREEEERSEERMLERNILEIENNE pH HL DE EE LE aaa BE EEE Ee Ee sses sonsarnar ap crasb EE BED £8 £ 4 : a h . 4 H 1 =7 =z a H Hl El Ze $8. H F : 3 1 8 3: 2 2} : : 3 de: ° vg; 3 - z sks 2 §- ELE £ 3 - £5 ERE %: 3% . z H * PE “Z: § 2s @ 5 ¢& : § t#% . gi Els 5 § xt tz 2 °F ? 3 : 3; 3% 3: FE:r Gf Eat g 2 H © ui- 7 §i-¢ : 3 g-t . i £ SHIRE Ts Em J) NERS. ® ) : i €§ 3 FE =p Ft snae.. ( 3 : F z gf 557 EBay BE BLE. Co 3 $5 88. - § B28 330 oof of; Edel 3 8 I 4a S wf 2:2 PIS §.8 $5838 pi oC

E 158 careronsnpzss seine eLEtR2TRtE PE PTR i nh RS REESE YY fELEE feed ee LeeEeaReeiaeass

ot ,

a : & E * iE 1 . £E £ EE 2 eo se § &% EF 2:3

es

i ie if £: e 3 ; g : § 33 :

aS. ° s a - - od ‘ 3 | oii e = i 2 Ree s s 8 - $c : g 5s elEt c Pia Fy < - 2

g § £ : Etre ¢ gE 8 gE SEs. = 3 se st: © gd 2 g€ Fiz z H 52 < 8375 = I 5 Ei Ee Z 82: 2. : 338% 2%. § & » £322 228.%2 Tesls 222.%5-C § 5 Zmgs & ESogE- FRITZ STOoILE $8 30a 2 ATURE NE Sela SEER FT oo-eEREegefo-to.miiilos 2b 7 veraTitoioTEni.viaziiEt gs al Ma Ez Ee s858er 2s eel vuuu Luu LY wi oad v

APPENDIX C

Y ’ FLOW CHART AND COMPUTER

31 @ I TS —— H MH . vo|EE = —

= g i | {22 | El id | 3 z 133 | ir —

— ist] 1 [3 fiz 0 = 1B'2) fs? | : ) “4 rs I . ey = H.2 ER HEE) Fit A TE gta | 4d TT esis |S A | Es os iE dag gd, TEI L§3d:% F le gp XS 35108 | f | d [77 8 =| BEd 9 KEE RTH PRY ANAY HE STH ANE SE of 35 57 500 dp DANG TG MET JE] pa emiin VV odie yl Ena Sais Praia of ER i=l H FFE PEisaz {<3 | : ro SE FT SRIANE. ._. [#7 | 1 ait! Becond

5 0 : . ) . ) EES. 2 kB o% ¥ i kid A en [A NAOH nai : INE /AANEL dL -] Sali :

— a od 5 3 $7 o

mr 8a i — IF xa £:4: ALGER) - i Ri TEE [3 fit] rErEEiiE = A de . HH Boga! | gig : i Eel : HH Hell | i gE © 1 A3- 2/2 32 55.82, {| 4° vel [J] : £%L v CE: RTH ANE JR £1 2g £3 " Bid 13,18 Ls PD fix 2) 3 : i 33 a | 18505%2 a \7/ PF 2 “Ny ’ | 15, ;22| EGY J “ Bd fozezd- DY © Ts a SEEM EE See NAT EET lm 8 mn 2H HEM CRY weve | HEE) I DAV ENTITLES $2 5% Ei 3 - [3:|D>—— (DHiFDlii 12 a | - | WE a I i co i | i -r.. : Nel] MARA A: LET Bh P on 7 ean HF A He tan A 7 i] i [URE 1B) EVE 7 aliaiififade % ed ci IEICE, 2 EY

h ’

lt 3% Ro ry 3 = pif) HH I=. 3 TRH INR 4 red sil _Lsfva 1 | 4 '1 Progal . [i] i ir iH | §1 i F 1d EA FEIT LH $i Lid | 3 1387 [sg 45 23.15 sit Fal EI TH pe [hd] Bl {| Rp Bay Ba = Yr |B . BR IRE FH | | 4 §| — RE g alo [tl ANA ES As f[d 2 : SHB [NH ih GY c Hit] = IEE IE IRENE ART isd a TF —— YS : | 4 ; 4 q Wii we OF i Alf H TE Ae L/NE 8 Ani Ten MG Es a LaF vam VASE IEA LIL it Hl 2 : il” p— E i ig - [7] ivi PER B $A [+] HEATHEN : [25 - FET Eh c o | & Filgrlisiiiiiag TN. El EET it i § orig aie ’ ol = : PEs 8 tasty 8 t Ti [3 0 [1335-% om 3 [25:5, [1] . mn GY PY Af oe iF: thet figs) oF i fre m= nail £5 i pre) ry) Mm LE 3 ik ~~ t B D\o/8 | 1 .- I A | ! 5 Vz &l A zs? m2 0 T4 d he > i BVANE! iV s Le L\2 Li elie ots fE 0 oH HH HOE0 AVAEE {LEB HAV 5 : BEALE BE VERY Vx ° : : : 1 3 .- [4 ELA AR EE wl oe i ft | a [ 22h i tH) EH hE! Elis] [283 == >> a La , ENA aE a RE 1 EE } Bey.2* 1d

163 ssersgiten at RRsR RR EReTataTIRILEIEI FARSI ISLICIRAREETEE2 REEL EEEEREEL ERLE RREEREEEERAEALEE ELE] LLEEEEREENEEIRAREREE $335335033504 5009040400 F SAF R ERE Tana FREE 7313225053220 9%;

  1. 2302253239505700005 0000904504343 000 0004; Seas rEra rans Re ae pt oe 2 H si oaME Ss ged: EO. I, 328 yz 23%: iF mE.n Fr3 ev Sol 3i¥s.E "E $230:

BAER ITEoarInIlsssEaNSRRsEnan IR nIRE san InSb IRE EARNER EREREE HI AT AAA AAA CAAA AEE R RE AAR TAR foto def ff fe fe ! EE Err rr Ee TL LE $ | H $ $ $i 2 $ % Bal 2 HO : 2 leg : 2 £382 | Le ite $a. IFEii.: 7 3 § F533 7% § ut 35 ete 3 ¢ g 5.2382 3% 3 143 $3% $. 8550, £ § £ $:7323 3 $2 i= 123 $28-3.2 1 eg E353 i3 3: ied 12% $SE.8%= 2 § 8 GE.3%: 83 3s ind 228 3 23 $ 2 8:32:23 38 HH) iz7 $2528 = $ 8 Popses 323% ..:3 ne $88 ....> $5 Safe = 3 £a8saz3.8 £73.358388 if: sade Iysits $ PP2E0R223 A.3p30R3%% INL: $goRefss -Iz%.: 2 $3 PRI RS83 .Co2RITRIIZ IRE HIER ERE g235..283 pEfTRIziic iiss $3EES=55 og0s £3 g53s22552 Ln02TIRNIIC 83S $eSSsags TELI8 $ © £23869: 7 go aise, LED $8538258C fgTc.tT 3: § . Ht ER at ELH $iEEEFC SSIs ey fi % S52227a%58s TroLpIsATECT IRL g dsaBiise EUL,S 03st S232 .8368 .PPIL3TRTIC. ilk. § iz8zz2: of? 38353 32728350 AFIT R.en: FNL § 33825722 3135s Stic FREELY SIT355255558 pth Ss 2372997 INLLSEILia FRELTEIST SIIILRIINITT iNT rE HEHE EARTH FETT 327 VT NESS oy S 2% .n.eas 3°5.58.-585238 Raves: SooEEL VF .*y a3B33233235 le SEs nih $22225230 Saa3%222.38, 15.7 g 233388758 5. 0v%. V0.0 383023322 "3333500800% eo: TaRisisEIT WrEN,. 307i: 383233333 piidicierocs ize” $igi33gsor: ibotzatElIc” IRIR3I333 S3ITISIINACI ieY. Ec8isakezas 3°L.° $58.0 Hi tt HR EE 303830 eEr $0500 -2p2) 248 iR32332332733732222008 Ent JGrseszechd ¥ESOoyS Ill: HR Rt RE EL $.23°3855%2 £3558 303.8 P3023 2322223830200008 oY HEHE: 11 BE Ol TTR pEEBITEE2ISIEAIoEAREE RIF ZRARIRIRERLE lerl.IiiiFiTT FEE RY EA LT RT TRSS 130558: srog ores RLTISG IE SIE IRIN 5F. QETLYITIINNIG ELiZ.c-Told SERPITEERIAfIITITYS ir VIC IEESTarozil e $3. Ladi SSERSSEFSISE °°

PEEIIUVEY LR 18. LE Tes SPRATT .A2SRS2 S.58520230 e Zuo oiEEL 2L32 IT 1&8 STC iil caREigts Tons -fi Aes: ? Spo EEX Sa8f ITE SM _S0T iz. CoatEIRIRiTRY T1838 z ESETFLIS: FUT ITPIAIY hye IRRITTLTTALSY (TIININ £02 S.igpis 3583 59: .t..88.;. 883% Ps xpfPIf 3M TUtSirtaiRNs.fRiE Hl I tT 8% YTIe FEY ~.s8TRTRSCRTSINGC H Eg igor iepEmiz 3087 S00. FNPCR 23 7 iS H £0: ff SEEPS 3.5. S3T80T80CA% 2% T 3% 5 Fao ir, PREETI bpla.NTRRERRITTT OTT TiS, | sossssnssstesisie LASSIE LOSSILLLLILSLITILIIZLLLL |

cent tet ALPES EE SSenESRCRE I PIR Las Rn ILE es ISR Iy FER 333 aa ed TE ELRRRELRRRAR RR EER iREini sian ibibbibiiiiisiinibabibithsn in pr EE EE rr ee EL

g rg 2 F bE | H 2 st

x a. ' ’ iy . z 2-% c | i . 223 . s H . oo . H H : £i. - s ) y M4 £ HH . § ? . £5 - = . E s $3: F s p H | | $ E> s £ got eo EE & H iia Se hc N . HY 5s s- =s53 : H v 2 = « a= er’ eo ir 223 H H g 3 gs fe 3 gf. 8 a 2s § S$ - FH £3 3 FN Eis 2 88 iF § s $= s 2 £2 pot SET = y= §s3c c= E32 . gz $s £4 ws 285 © £2 2. = £3 . Fy . Ey - &ee es? eS" < - 5 Ld Sz <3 - 3 oo © STRELEPI2a” z Bes 23 ssfonze $2 o BF. piE.c ERENT 3. 23s i 2S ediiila- e 28s’ & ViF%e Tohicepsebgd "0k" 3 %-_£ 35° $5 gipisEiz of TRe-2 yPTS LLS¥BliogPeniassviiili SIR ® Tare rE are ri) Se ret Ci fame TTS seme Teen € E4%%..o%0000% 1 B’r felp saflechs, FES TI SEeP Tije c spbocidionE; FoRFs eBEs FIInELSTT 88: °8s 27 i z LI s s g gE 3 H | i pe “ |

ares snser oan ren ERR nana RARLT ITI] rE tit sti t HH EE Ea E ERLREREERLRERERERLARERERARIRAERAAIAAS AARNE, LERLLELEELERLEE SE 210205a705550005550000 a00abebaREEE0 DEEL EE 00

  1. 5 3 38 H ¢ ¥ i § z ; ¥ Ee go H : 2°8 H 5 gE § H : g if Ea i s :; 2s 2 | v : it = 5 § £ & é* . ir Ef: » “ ¥ £2 Eos . ££ BH 2 a i 1) gE SE A g 2 3 - =% E35 : £ . 55 - : . og SN. 3% °F d 2 : $3 2 iv an. 2 5 - [9 ¥y it : TE A TE Es H = "eo Te £2 ~» 3 = s $ _— = 2 82 so Fz Biz 2 3 § .s - : § Fd £33 3: 2 TI = = ig 2 £« = 2. ic B EZEc&t ¢ 5 : 3 . £2 3 3%: fo: %= ¢ g E of 2 52 = 3 $2 £ » Sez 3 2 se 14 $ . eT = B It = Zz. 5 3g «sf § uf == LB SE 2% cf EB : 3% % 3% ged 5 cE Saf . 3 § & 23 1 BE FB, Cp foe TE Ex 92 £3: Ey Si scsesb. FP $:2% § EF OE YB Fi. NIN 3.5 PefsEys FB TEU $3 OF ST. 5 Soff fISS £5. ¥r227:8 7 0% 28.7% % od $2 ,oPzr S i837 teen gies, E3fiiaze’ § Zee I IC $8 Iziil oy 83%7 2iIiz " $58le EFL 3005, T 5.585% Bois $% 35 E°%c0% ~Bog 3SURZ.2 "esa olsosel Be Suk . sop et STE TET PO e2Z7e 25-55% 8o0%.33% err en” 2" .5 380” 0" %%0 22d Sgteils :Poi a: 2327 SITE pi | FI S + ils *RRLTT FEedTip 20TTUE 25527133 szed” Boats 33 2 =8 oF Cr PALE LT SRL SS Shlens Is os £2 2111 Sa IY tgtsisee fplls E870 TES:iliE=~r IE i 5 ssc: § zs es £2 ] 2 3 : iizz H Zi: = == £ gz wou wv we [I “= TT TAs Tr rrr rr ah FREIRREizIARARRIERRARARRRARRRARERIELERnEnRaraiziiiiiici HR A ELE ESRLRLELL ERR EE RE EE EE EEE EE LEE EES Cr i o i ¥ H =] t

2 Eg

22 ig = Es 3 5% £5 Ze nd i ne 5 . ™ 2 fs, —-— = EH 3 - i-8 FT H gE £5 4 35 fz & H 7 ot 11 4 ct? ' & Beg KH = ES 2 : gis § Ei: 3 £2: 2 3 £22 € £55 J 5 id : ges 3 £8 . 5 se ge ¢ 8252 i g { : 2e2 5 & £3e0,2 £3 - 3 { i 535 EB §7:'38c ei. IE 13 == © a = 55-209 Es cs % 4 St © 20% SeSgied o- ¢ % sg 2: cg s S¥sbc Rs £3 £3 ft Se . at 3 LITTER Sef 25. 1 “wr ET =" % Sesolte kt: 2® 2. » hE = p= TET BepEEts VHS JZe3: & E + = Ete M3 cis e333 £ } &3 33: ot 8 d8 1 ,3zecsz; sg: 233253 fas ° S32 gz, 3.8F = Af3elS5s 23ns gi=23. ! £23 SST ®eictlLT o. SMS-lTER ESE: rs $53. I 2, f8%3E, iT ieleifsesi LSEf oEzi7o; i RAL a Fits LX Th Z8%3 SRFR, E-EIIESEREETOOIIIELS 4 FITIV — IZ e032: i EI Dr Ms Pe SEE ol LIB J ee Chl a5: REvas GERRICoCIoTIoSB.Siialal SRT. 3 TAFE NEC Et top r 11 1 ht £0 TI + LI EL 2 gEbss-3.939323 m3 a8 sra8) 7 Ziedlelrlalr g 1 108 EN LL Ne hips J 17 pt Ie $ s23.83e3lsTe In Es Ra ’ 58 3% =eJ2E8PETRES 1 eg Sh o£ 2 8 8 e 22 £2 3 ° i s2 ss E33 = H ® Hl Ss 2: =3 - RB wo wo |

i ar srenspEear ESE at SESS S0SRSsETarararsgenanaz ante RBI RSS 202 TT ra ER EE a HR a a [rte TUE ERASE PERE REE B | EE

£ ES # : 3 pt ¥ : J B TE = EH H p: g ws ® ° 3 : } H : GE H EE z § . - $ £2 4 : I FE e : 87 Hi g : 4 3 gz . - HS 3 'S / . 2 . 3 . oo: . 3 3 $ 3 : fe B e® 2 po 3 = . a . ¥ i Ele 3 A ° 3 z y c 1 i Sok = T° 2 : : § : : 87 1 : Zz es 3 = < ‘os . I 383 8 8 3% £5 ¢ 5 £ = = 1 31 $88, = 3 8% of T = £ g 8s Nz § 1 1] Sa® oo S&F Tos Se = = es Sch of =a Tif x FE < 3: T= zE 38 ie gesat 2° $c! ! ddd; ob 2c iEtd L333 o£ Eel Ee pIRCH x 111 JS BE Bgeds Cg ob. 25a iREel FRG ’ feE= F 3 EC £ vr FrZlRoERaifeT iE Fe B0S0 LT So 33 EET: i SSeSTh. LeFit EF ge TF weedTos2.iEagic nf Bie. 0iad BRFECE.2EE TEPAIieIiniced Big d EE c3REILIS TESTE SEt.32% Folic Bocn CEEfIasIaizesi-gis BELT 3 TF hr Ppt tet 3 PRE Pe BR TINE eEaTIAYIRERIL.~R32 LgroTEsT 83 PE at Tt BEA bk 1 1 EERIE TLIC IIE 1 CPt Sec =F St .} a Re a EE PT LF TF Ii Io. ti oBITSLERaBTNtElRNE Tol08B JB01eR0 T5700 abe cig Fess so, er n03E8{0F J.22 BERT B r 278 gil $235 § 2% €% § HE £2s gs? © £21 | 232% & 23 2: 8 $8 233 822 2st |

REFERENCES » [AHL66]) L. V. Ahlfors, Complex Analysis, Second Edition, McGraw-Hill Book Company, New York, 1966. [BIC71] T. A. Bickart, D. A. Burgess, and H. M. Sloate, "High Order A-Stable Composite Multistep Methods for Numerical Integration of stiff Differential Equations", Proc. Ninth Annual Allertown Confer- ence on Circuit and System Theory, university of Illinois, 1971. [BIC72] 7, A. Bickart, and Z. picel, High Order Stiffly Stable Composite Multistep Methods for Numerical Integration Of Stiff Differential Equations, presented at SIAM 20th Anniversary Meeting, Philadelphia, June 1972; also to appear in B.I.T. [BJU70] G. Bjurel, G. pahlguist, B. Lindberg, S. Linde, and L. 0dén, Survey of Stiff Ordinary pifferen- tial Equations, Department of Information Pro- cessing, Royal Institute of Technology, stock~- holm, Report NA70.11, 1970. " [BRA72a] R. K. Brayton, F. G. Gustavson, and G. D. Hach- tel, "A New Efficient Algorithm for Solving i pifferential-Algebraic Systems Using Implicit pif ferentiation Formulas", Proceedings of the IEEE, 60, 1972, pp 98-108. i [BRA72b) R. K. Brayton, C. C. Conley, Some Results on the Stability and Instability of the Backward Differentiation Methods with Non-Uniform Time Steps, IBM Research Report, RC-3064, IBM Watson » Research Center, Yorktown Heights, New York, . July 1972. i [BUT64] J. C. Butcher, "Implicit Runge-Kutta Processes", Mathematics of Computation, 18, 1964, pp 50-64. i [CHE71] L. Chesler and S. Pierce, The A plication of R Generalized, Cyclic and hodified Numerical Integration Algorithms to Problems of Satellite Orbit Computation, System Development Corpora- tion Document TM-4717, March 1971.

169 [DAH63] G. G. Dahlguist, "A Special Stability Problem for Linear Multistep Methods", B.I.T., 3, 1963, $4 H pp 27-43. [DEJ67) B. Dejon, "Numerical Stability of Difference Equations with Matrix Coefficients", SIAM J. Numer. Anal., 4, 1967, pp 119-128. [DON71]) J. Donelson III, and E. Hansen, "Cyclic Com- posite Multistep Predic tor-Corrector Methods”, - SIAM J. Numer. Anal., 8, 1971, pp 137-157. [DYE72] J. Dyer, S. Pierce, R. F. Haney, and L. Chesler, Generalized Multistep Methods in Orbit Computa- tion: Studies in Existence Theory, Efficiency, t Optimization, Systems Development Corporation Document TM-4888, February 1972. i [EHL68] B. L. Ehle, "High Order A-Stable Methods for the Numerical Solution of Systems of Differen- i tial Equations", B.I.T., 8, 1968, pp 276-278. I [ENR72] W. Enright, Studies in the Numerical Solution of Stiff Ordinary Differential Bguations, Dis-

[HEN62]) P. Henrici, Discrete variable Methods in Ordi- nary pifferential Equations, John Wiley and 4 , Sons, Inc., New York, 1962. [HEN64) , Elements of Numerical Analysis, John Wiley and Sons, Inc., New York, 1964. . [JAI70] M. K. Jain, and V. K. Srivastava, High Order stiffly Stable Methods for Ordinary Differential 4 Equations, Department of Computer Science, Uni- B versity of Illinois at Urbana-Champaign, Report No. 394, April 1970. [KRO70] F. T. Krogh, On Testing a Subroutine for the Numerical Integration [24 Ordinary Differential Equations, Jet Propulsion Laboratory, Pasadena, California, Technical Memorandum 217, October 1970. [KRO72] , Changing Step-Size in the Integration

11 _— A Tr "Blu det T of Differential Equations Using Modified Divided Differences, Jet Propulsion Laboratory, Pasa- dena, California, Technical Memorandum 312, October 1972. [NOR62] A. Nordsieck, "On Numerical Integration of Or- dinary Differential Equations”, Math. of Comp.., 16, 1962, pp 22-49. [NOR69] Ss. P. Ngrsett, "A Condition for Ala)-Stability of Linear Multistep Methods", B.I.T., 9, 1969, I p 259. [ODE71] F. Odeh and W. Liniger, "A Note on Unconditional Fixed-h Stability of Linear Multistep Formulae", i Computing, 7, 1971, pp 240-253. / [RUB72] W. B. Rubin and T. A. Bickart, A-Stability of ( i Composite Multistep Methods, presented at the SIAM-SIGNUM 1972 Fall Meeting, Austin, Texas, October 1972. [RUB73] W. B. Rubin, A-Stability of Composite Multistep Methods for the Numerical Solution of Stiff ’ Ordinary Differential Equations, Dissertation, Department of Electrical and Computer Engi- neering, Syracuse University, February 1973.

; |

17 . [sLo71a) H. M. Sloate, Simultaneous Implicit Formulas for the Solution of Stiff systems of Differen- p | tial Equations, Department of Electrical and Computer Engineering, Syracuse University, i Technical Report TR-71-4, May 1971. {SLO71b] H. M. Sloate, and T. A. Bickart, "A-Stable Composite Multistep Methods", presented at The i 1971 SIAM National Meeting; also to appear in Ss BEX. [WAT71] H. A. Watts, A-Stable Block Implicit One-Step Methods, Dissertation, University of New Mexico, 1971. i [WID67] 0. B. Widlund, "A Note on Unconditionally Stable Linear Multistep Methods”, B.I.T., 5, 1967, I pp 65-70.

172 BIOGRAPHICAL DATA Name: Joel Marvin Tendler pate and Place of Birth: May 18, 1943; New York, New York Elementary School: Rabbi Jacob Joseph School , New York, New York, Graduated 1956 High School: Rabbi Jacob Joseph School, New York, New York, Graduated 1960 College: The Cooper Union for the Advance- ment of Science and Art, New York, New York, B.E., 1964 Graduate Work: Syracuse University, Syracuse, New York, Graduate Research Assistant, 1964-1966, National Science Foundation Trainee, 1966- 1969

4 ‘