## TendlerDiss1973.pdf
i A STIFFLY STABLE INTEGRATION PROCESS
USING CYCLIC COMPOSITE METHODS
by
| JOEL MARVIN TENDLER
B.E., The Cooper Union, 1964
| DISSERTATION N
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
‘ .
~
[
5
v


|
i
| A STIFFLY STABLE INTEGRATION PROCESS | .
[
| USING CYCLIC COMPOSITE METHODS i
| 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
| -
5
-


; (
i
i
ABSTRACT
i
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
ii
«
|
!
.


: (
i -
4 )
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 i
| 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
111
5
- .
- 3


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
¥
Es


’ ¢ ~
| A STIFFLY STABLE INTEGRATION PROCESS
USING CYCLIC COMPOSITE METHODS
by
| 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
i
{ A
| AJ
| {
|
bh}
v -


£ \
/{
[1] )
| To the two Philips
in my life—
| one of whom has affected me
i as I pray I can the other }
AJ
| ii |
|
%
H v


(
i
i
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
Ii
. L :
g -


’ (
|
|
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
I i.
|
i
ff iv
(
|
A!
. pe


: (
i
i
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
i
i |
A!
wv Fe


< “go -—
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
{
: .
5
& [*] .


. » a
: (
|
1
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 .
| ;
|
{
AJ
vii
X o
¥ o a


: (
i
¥
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 ¢
I
[
.
|:
| viii 3
Il
h
b .
1


« 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
=
]
¥
v


¢ TN -
: (
I
¥ 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 (
.
I
[ |
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 *


< i a
: (
| 'e
H |
| 5
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.
|
I f
l
:
.
I
[
Ii
5
2 .


<< Ny =
(
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.
»
5
o Fay .
.


< i) ~~
i ’
i 7
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
I i
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. }
A!
LJ


¢ i -
» !
8
2.2 Linear Multistep Methods (LMM)*
1
i { Linear multistep methods may be written in the
form
1
| 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
1
| 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=o

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.

5
v


« i“ ~~
5 4
9
where
|
l 1
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;
- p+l wr E
| TURE Yltn,) + omP*t), j=-k+1, ..., 0; and
iii) y(t) is (p+l) times continuously differentiable
[ p+l
then y .q = y(t.) + O(n 3»
[ In the light of the relationship between order and ’
local discretization error one would like to maximize the
order of a multistep formula. '
| 2.3 A-Stability
In 1963 Dahlquist [DAH63] introduced the concept
of A-stability; namely.
pefinition 2.2: A linear multistep method is said ’
to be A-stable if all solutions of (2.5) tend to
| zero as n+ when the method is used to integrate
the differential equation ’
\!
| . |


‘ Ww ~~
A 5 {
10
| )
|
| y(t) = gq y(t) (2.9)
with any fixed positive h where Re{g} < 0.
Substituting (2.9) into (2.5) yields
i
. = qhB, . = 0. (2.10
; I (oy -ahBy) Yn 10)
| j==k+1
The solution of (2.10) is given as
| 3
n
| yo= I 2% (2.11)
i=1
| where, with gh = ), the t's are the roots* of
| 1
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]
sm
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).
AJ
|
|! -
E v


< Ww ——
4 [
i 1
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)
lal
———
1 #*The characteristic equation is discussed in detail in ’
| | Section 2.7. |
5
Ed
¥ sd


‘ » ~~
A {
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
a 1
| 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.
|
|
5
| pe
i v r


¢ Ww ~~
Li
1 [
- pen -
| 13
circumvent Dahlquist's result.
2.4 Ala)-Stability
| Widlund [WID67] introduced the concept of Ala)~-
| stability; namely,
Definition 2.3: A linear multistep method is said
| to be A(a)-stable, ac(0,7/2), if all solutions of
(2.5) tend to zero as n+= when the method is used
| to integrate the differential equation
y(t) - q y(t) (2.17)
with any fixed positive h where gq is a complex
| constant which lies in the set
5, = {z:]arg(-2)| <a, 2 # 0). (2.18)
A method is A(n/2)-stable if it is A(a)-stable
for all ae(0,n/2), and A(0)-stable if it is
i A(a)-stable for some (sufficiently small)
ae (0,7/2).
As the set Sq defines a wedge in the complex domain, the
angle a in Definition 2.3 is referred to as the Widlund
wedge angle.
h
EY ie
© ;
- Ee


< iL ~~
/ (
BN ; : (
3
14
|
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


< Pel —
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
1 a
ky = fly, +b I Ajgkyr tq * czh). (2.19b)
i=1 i
i
| with i = 1, «..r 9»
and )
q
| c= 1 Ay for i = 1, «ees 9- (2.19¢)
j=1
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 |
de


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
|
|
A!
Ld ao
1 ~


< oF i
17
l
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
aa?
¢ q
| *Both implicit LMM and implicit RKM may have to be solved by
an iterative method resulting in a higher number of function
- evaluations per time point—1 function evaluation for an im-
| plicit LMM and g function evaluations for an implicit RKM,
per iteration. However, on a per iteration basis, an impli-
‘ cit LMM offers a distinct advantage over an implicit RKM.
I tA variation of (2.22) has been used in which (2.22) is used
to solve for & future time points, keeping only the first
| m (<t) of these values for Lhe next iteration. As an extra
(t-m) function evaluations are performed at each stage, this
form of the CLMM was not considered in this research effort.
| SNote: The index for n = 0 and j = -k + 1 is not zero. To
make it zero results in an unnecessarily cumbersome equation.
It is left to the reader to recognize that to make the equa- ’
| tions precice k - 1 should be added to all the indices on t
and y.
5
a,


« 4 TN
§ 18
L
¥ 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
as
I p
> i, (3)
Ly ly (tp) hl Icy, 407y 7 (Ey) .
I 3=0
I with i =1, ..., %, (2.24) ‘
: where
| L
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 .
5
3 : 0 3
pe Send |


[ oF
19
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;
nt
. 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 ’
py
|
|
|
| | |
hi
i C
9 a


¢ 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.
|
a
C
gis,


« of ~—
. p :
l
2
| 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
1
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 ’
5
: 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
0
I Ry (NY, 5 = 0, (2.29)
=x
| where Ry (1) H Ay - AB. Define
0 2
| 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 '
KL
n N
{ , = 1 Cie (A), (2.31)
i=1
where the vectors Civ i=1, ..., KL are functions of the
’
A!
CTE
k my


a. re ~~
2
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]
0
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)
=J
| where D4 (2,2) is the greatest common divisor of all
g minors of order j of the determinant of R(gz,)) and
- D,(z,2) = 1. The general solution is then given as
I [DEJ67]
8s L m
i,]
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).
Ad
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.
|
|
| ’
|
|
x .
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
5
. “3 a
iE a


‘ rg —
/ : (
4
2
' 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
i .
pg.) = 2(2) P(E.) ADV, (3.4)
a1
| 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). ’
i
(
: .
| 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
3 >
defined, corresponding to a zero of A(X).
Let B()A) be defined as
|
( B(x) = 50 Ns) - (a: AN = 0}, (3.62) po
]
| 3
|
|
. 4
‘ a
1 a :


‘ 5s ~~
29
: 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
i3
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. |
|
|
5
.
= Bn .


p oy ~~
i 30
or alternately as
i :
ple.) = I pytont, (3.8)
| so
| where
k
| pyle) = I agyed, wien i=o0, oot (3.9)
i=0
| Similarly, one can express pl(g,2) as ;
| k
plz, 2) = [ aye, (3.10)
| -
I where
L
I ayn) = a ts with 3 = 0, ..., k. (3.11)
i=0
| 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.
|
5
@&
1 une


p 4 ~~
/ : (
31
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.
| ~
5
Te
he k


a. 5g —~
/ g (
33
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). '
A!
< ©. a i
- ts.


_ 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
L
Cn
1


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
I
| 650) = 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. ’
| |
| {
i
5
Car
~ §


“__ 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).
and
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). ’
v
Moreover, if N,N) and N, (0) are two such regions
:
; 3 | €4 = |
\ |


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


<< bg ~~
| 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)
or
p(g,A) = ¢l(z=1) = (Z+1)A1(1%A). (3.23b)
In the notation of (3.4), '
5
@l
Ek hk E


_ E38 ~~
/ : : (
| i
»
| z(g) = ¢» |
| pg.) = (g=1) = (g+1)X, (3.24) 5
and
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 I
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 : . (
|
40
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
| |
7 .
Cw ;
8 —


p F's ——
— \ <\
41
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,-
’
5
CIN


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
| a
#If the CLMM is to be at least of order zero, then each row
’ i of the matrix
0
i a= J] A
i=-K
| must sum to zero; hence,
0
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.
| |
Ad
5
@F
hs


l x ~~
. 4
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
T d
| 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.
h
CT
\ |


a F's yy

i 44
[] —

. i

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
{
{
A!
k Cr
Wy


_ 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,.
|
SE
cw
: 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 |
5
= 2 CT 2
2 NX


- 53 ~~
4
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
m
| | 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.
(
A!
Cn
iE


< i =
/ : (
| {
48
|
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
| [
I
5
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. ’
[
| I
5
CT. :
7


‘_ ox ——
— \.. 4
50
| )
-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. |
bh)
; CT


‘ « ~
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.
h
! SE
Ek : ¥ 4


{LS
: 52
- . []
T4888:
' 1ji:qii1,
MEE ER ER
THLE EEM
Fitz i573:
Sextfilol,
§. 83,439
Sfus EEE
Hit
Es5pg%° .
i |
i a
[+]
s
2
® <
i 3 =
280 28; I]
$2. ; $2 4
= = $3:
| gli: ji: -
& Ba 8.3
i LEP dgiz :
- a
Exif 1H g
148 zt a
| iad o
=
Eel
» : |
Be
| £1 8
2x .
: ,
I 3 ob &
EEE J
“ea on
3
| tH :
bi :
: -~
i)
i. of
$2885
1883
| 2883
S &
cand 3
§iic3
HA FE
| EME
58%;
st gx® '
FERRE
i
I
|
1
“ g
5
5
$ hi.


< i ——
— \ . \ 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.
{
{
5
@h
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
\J
.
5
Cw
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+) !
BE
: (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
{
A!
Ei
3 es.


#* (
» =

| 56
| 3] :
53) 3
—%-i ®
| $1: | [4
I | rE | g
4 s
, | e
; A les) | -
t\d IT = pel
‘ A EA 2
- 57] [9] J wi) | 2
! zed fig : | EC
HIS t EE — 3
id: HIT 2
I HALEY fizz] Jia i/ 5

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\ at [9 »
RENEE : EA FRA Lie ER
ETT 79 gl Ee) 13 dell z
pS E80 ¥ 8 py we boom nd

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
- 5
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-
3
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.
-t
F puring the start-up procedure, shown in flow chart
N
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. ’ |
3 [
o ———————————
a f
a i *Indexing starts at zero.
Dp
g The locus for © varying from 0-7 is the complex conjugate
| of the locus corresponding to ® varying from 0+-7 and hence
- need not be computed. This can be seen as follows: The
* characteristic polynomial can be written as
p | kK 2
BE ple, = I I agnied
a ( §=0 i=0
where the a;4's are real numbers. Now, with ~ denoting
| complex conjugate, consider
gE kt ’
pg.) = 1 1 ag sae. ’
| §=0 i=0
5
@i
Bs


‘ 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 “ *) )
<T mews | a0 w
a
WE x I—
| 5 ele SL 1
Bs | - | 3 .| 1
. J
i a: | ro = 3 | Pp
| - J | Gremd=e
hr | <P
oN | 1 | = .
f =] | wiry
| TerAR Te | —
} EE { Wer) nd
|
| Jo——— i )
[ S = ee
| Zo Y
E ey. ©)
Voie Genes
Lee—y 4 =
Gar> : | Seane |
A C5 © ~~
har a Nv
Tor ©
u
E
0" a ds 3
ar wary Suny, SaOnn y
| Lei
| Figure 3.4: Main Program——Start-up Procedure
{
l
5
Cw
EL v


|
59
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.
i :
AJ
|
5
o & V
3 bs


# —
\
: (
60
“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. . |
|
( |
4 5
£ CC ® :
4 \


1 \ ( {
61
|
- is further refined by a second rootfinder, REFCROOT. (This
| rootfinder also uses a Newton-Rhapson iteration and improves
on the guess of the root locations* as given by COMR™ Ts.)
ttn
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.
J
{ *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.
[
13
a :
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. ’
5
i. v
3 i


KL. HE _
64
- E
t : .
ig \ Lal Gli o
| off
| |
oF
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 -
it
= \/ $35] | B
K |
ih ) 8 I fl Bt
© ii hi 31) | 135]
i" r=
: :
Te
Bo,


\
— \ \
|
| |
| 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
i
- 1
f
. .
. 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—
i
| 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 (
{
= (
|
(
|
5
GF
\


\
4 CHAPTER IV
A NEW INTEGRATION PROCESS
[§ ’ 4.1 Free Parameters
A composite linear multistep method may be written as
L
| ) (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
K
I plg, A) =det{ J (Ay ~ a_i")
i=0
| As 3 exists, the latter may be written as
K
plz, A) = det{agl § (Aj - B_0e5 3)
| © | |
3=0
K
| - det{A,) det{ J (A_; - AB 25)
0 w -3)% },
3=0
{ 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.
|
A!
CT
3 is.


\
(
t
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.
1
5
Cr :
: Be,
J |


+ ~
: N
| 7
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
K
| p, (2) = det{ ] a B (4.4)
i=0
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
1 N
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-
Bd
{ tiation formula (BDF) results [HEN62). By imposing the ]
l I |
condition that all ¢ equations have at least order k and i
c
|
5
CC. ®
~


« 7 2 ~—
(
!
| n
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,
I
| Yps1 = Yn id ep (4.7) .
and
| 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.
{
|
A!
i :
\


73
| 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
|
|
5
el
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 ©
- LE «378 1080  -a20 140 -3 3 [] 9 0 Wee °c ° [] °
3 «32 isk 472 2 #72 e168 32 3 H ®s 0 ee * H 2 ®
.3 30 1s orp -1080 37 a0 che . : 9 ° 5 a peo ° HY El
2 BW oie 0% ome cee vm me a3: 2 5 $$ 3 3% HE
Li BID 08 N80 "Bao oases 1a od 8 8 8 3 5 5"5 ee: 0
15 80 39m0: heck 1a700 -asemc 11763 -4720 2283 : 2 3 5 8 0° ¢ om
| pees
PT WE NE TW Pa RE de ow ws 8S
27529 22683 +as60 70840 7938c 4330e 33280 123070 .o83% 280 E520 2 a - 0 ° . s
“tar Ta329 10080 11767 1176- -882p Tos “1880 ez 28 ° 2: : H H H H H H
Ed 4% 1905 ‘Seay -es1- 29g e1sd0 See wig% 19 ° 3 ese s ° ¢ 0 3 ° 0
ls 139 oto -155e 78 19> Ged c270 se 5 ° ° 3 2 EH] ° ° ° °
3 133 T1080 cles 5% wars geo Ben oe . ° 0 ° 2520 ¢ ° H 3 °
. WS egep Be 282  S0s 1680 +380 ww 0 © ° vee 0 2 0 °
s “5s 270 +80 109-378) 19% 1282 13 19 0 H El ¢ KH c ase 3 El 1
esp 308 eB0s 1870 290 seis aso 270 S90 35 E 2 : H H ° o #23 o ¢
. 12 300 1680 oars Sa2c 1176y 11760 -10%0 AZ) 245 ° ° ° > ° ° ° ° seo H ?
«280 203% «12960 35280 ~35ce 79385-70560 +5340 «22680 7129 ° - 3 B ° ° ° H t ese
Table 4.1: Linearly Independent Multistep Formulas for p=3, «..s 9 BA
| (continued) J
A!
TN
3 oe


(
i i
] 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). .
i
{
: ‘
| 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
Bitls*3) [2] 0 []
B(ls*2) [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
1
=.
ho.


< + ~~
>
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
~~
5
] CE \
\


‘ ® 2 —~
— \ \
; 78
| 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 |
| 1
| 7 33.53° | » |
! FREITEPIREIR
) I #The order 7 BDF is not stiffly stable.
i Table 4.3: Comparison of Widlund Wedge Angles
NEW
[ 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 ~
4 -
5 |
: ©. K
#3 —


’ : (
3
i
”
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. .
5
CT. 3
\


Bg \
|
| ruc, +
| \
1 \
SEE EEE
I
| / :
a | /
Ly
l R12,
|
a
« i
| LY
\
ae Ca
Pd
Mo rd
-3
Figure 4.1: Lambda Loci for Order 1 and 2 Methods }
: :
Cre :
a= :


81
I i
3
| 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
pe
-8
AJ
Figure 4.2: Comparison of Lambda Loci Between New Methods
and BDF for Order Three Through Six
; 5
- :
- -
, ; \


< ® . ~~
— \ x
82
RSAC
15
~~
— a
hp She
OF  % \
N \
\ \
\ \
NEW MET? \ \ .
FPIPIIPSPIPIPUPUPITIrR Seer ©
- " | r
/
/
"4 /
d /
cg” ”
/
I -15
RE46C
8 SR
Tn Ph .
+
BOF 5 \
\. f
NEW METHOD \ !
=
| = Lr
/ /
| /
/
/ /
Nd py
£
| -
-20
( AJ
| Figure 4.2: Comparison of Lambda Loci Between New Methods
’ and BDF for Order Three Through Six (continued)
|
5
Cre
N


d > ~
| 83
R747 2
—
| -
| 1 sass |
| Ne Pd
ARE
| =
[ Figure 4.3: Lambda Locus for New Order Seven Method
5
[€ “w
(4 \


_ ® 3 ——
> [
84
differentiation formulas of order 1 and 2 form the basis
| for a new variable order, variable step algorithm described
| in the next chapter.
| :
| -
I
[
A *


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— -!
; from 1 through 7-—used at each step are chosen in such a
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 ’
h |
| bh : |
: \


y 2 —~
. (
'
86
|
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 ;
i
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, ’
i
I
] "
| | | |
Co ;
1


< © 2 —
~ k
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.
Yes
: : . : ’
Figure 5.1: Organizational Flow Chart for Subroutine DIFJMT |
|
5
Cr -
~N
\


l + -
88
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].
v
x
re, %
IN 1")
; )


{ eS ~
89
(m+1) o (m) y
%5,i¥ne+i ws hy sEWngsq tarsi
i-1
- I (0; ,5¥ne+j - h8; 5¥np+s’
Juki.
| m=l, 2, «cos (5.2)
is performed, where y = £(y,t) is the differential equation
, i being solved. This process will converge if
$. &
i,i of "
| ||n = =i <1 (5.3)
Sod
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
i
B:
- h—22 £(
Yng+i Re 3 Engi thei! :
i,i
i-1
1 .
Pl 2% )) (05 5¥n2+3 hB; s¥nges)” (5.4)
'* je-k+i
’
5
Ce
N
x \ H
\ |


l + —
90
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)
ii
: (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
AJ
| n
%
— .
\ -
! 3) ;
|


{ eS gy
9
or, from (5.6b),
Bi,i
LL" =1 = . "WY I? (5. b)
’
I being the identity matuix and .
- of ( = &
| g_ = 20E) 8 (5.9)
{ m 2 bo)
y _— )
Substituting (5.6) and (5.8b) into (5.7) and re-arranging
| terms results in
Bi,4 m1) _ Bid (mw) (m)
| i - Re. InlY = vamp 4d Wt) = Jy ]
ii ii .
+ C. (5.10)
| Note: If (5.10) converges it will converge to tne correct
solution. Hence, the Jacobian matrix J need not be pre-
| cisely known; in fact, a very poor approximation may be
adequate. It can therefore be correctly inferred that In
need not of necessity be recomputed at each iteration.
Now, subtracting (5.10) with m - 1 replacing m, from
p (5.10) and letting LR ha 5” J yields, after re-arranging R
{
| terms, the desired result; namely,
Bias y
yo opm dd lee) - aj, (5.11)
i,i ’
(
|
5
oe 1, .
\ “
A 3


92
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)
* hts Jw imey'™ 0) - ail. (5.13)
Ed
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
i-1
| 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
AJ
h
ws
: 5


¢ #* —~
93
i-1
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
i
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
1
1 .
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
1
(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 .
|
|
il
\!
Cr
SN


\
\/
94
value of hy pei” Moreover, when convergence is achieved ’
$ (m) * m ._. . .
% Lye] = hy 43 that is, just as y_,.; 1s considered the
- (m) : E
m-th iterate for Yaeei 5° too can 404i be considered the
m-th iterate for BY oss”
Summarizing, the procedure used in solving (5.1) is,
1
for i=l, ..., &:
{ 1. A predictor formula for suitably chosen 5 3 and
’
. : ((}] c
* a s in (5
| 81.5 is used to obtain Yne4i’ as in (5.16), and
(0) e
Agi” as in (5.18).
| 2. The Jacobian matrix J is evaluated and used to
compute
By 5 _.-
| wlez-ndd 57th
i,i
3. 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,
* Fs a
°i.3 = of 3,5-17 (5.19a)
and
3
| |


IN
)
| 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
| 1
* - B v = 0
I of,5 Ynys hg 0 Yn = O° (5.20)
| j=-k+1
Note: To within a multiplicative constant the order condi-
tion serves to fully specify the a} 5 and gf ,. Also an
’ 1,0
| order k formula was chosen sO that the predictor might en-
hance the convergence properties of (5.1).*
] The calculation of Ww l—required for (5.11) and
(5.14) —is time consuming for a large system of differential
[ equations. Therefore, it is only re-evaluated when the New- - .
| ton iteration fails to converge. It has also been found
worthwhile to re-evaluate Re whenever an order change is
[ made. However, if only the step-size is changed the conver-
| gence properties do not seem to be significantly affected and
ci ———
J #The order 1 and 2 predictor formulas are, respectively, .
Yne1 = Yn +h Yn
and 2
Yns1 © Yn-1 +h, i
For k = 3, ..., 7 the predictor formulas used correspond to i |
those found in the second row from the bottom of Table 4.1 |
for each of the corresponding orders. |
| |
|
l 1


\
J
\
» |
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)
I
- is checked and if satisfied it is assumed that either con-
I vergence has been achieved on this iteration or would be \
| achieved on the next iteration. In either case the process
is terminated. Note: The look-ahead feature incorporated
| into (5.22) assumes linear convergence. The factor of 2 is
applied to insure that on the next iteration, if it were
(°F performed, (5.21) would be more than satisfied. If such is
A | the case the next iteration is not performed so as to save Ce
i
] ~ 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.
|
\
|


\
A
NJ
|
97
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.
: |
\


\
NS.
98
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
i
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
- for error control to be outlined, based on the local trunca-
h tion error, appears to work successfully.


A
J
| 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)
As
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
max
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 !
’ |
.
\


N
A
NJ
100
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.
|


N
A
J
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®
i N
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,,)"- .
- SNote: In the computer program listing given in Appendix C,
{ wy j=1, coer N, is get to the maximum of wg and | Wngeisle
i=1, ..., L, at the end of each cycle. Initializing uy to eo
max{l,| (y,) 41} for j=1, ..., N, provides an error test r
relative to the largest value in the history of 2 un- -
i less the user overrides it before each call. The YMAX param-
eter in the calling sequence corresponds to the weighting
array w described here.
—
|
\


\
J
|
| 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
i
( u 117 GG+1)
| EIN PIA" wl
+l Ing+2!!
A
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 hI*l 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 1199*y oul] is not usually
| constant the factor y is used to provide a "safety" margin.*
5
i FE ———
| *a value of y = 1/1.2 is used in DIFJMT. t/
|
I
\


+
A
J
I 103
4 Ld < 3
: © oo o >
- Oo = bed
o' mm Mm © Oo |
7 suse. 3 I ER
% ERX SSS]
a AN A cc oe ao oo |
S » ao ~ | []
~ © ™ | bo]
ll ® i
- | =»
mw | 2
© 2 o |
co = 2}
= m © oO | ©
> oN — | [+]
i A NON a | y-1
3 a oo oo | +
.. mo = | Q
X NN ww mm | =
|! |
. | | |]
| 1
<«
” [=] © “
~ mm © +»
- m m © © Pe
i ; vw oo oO oo 2 | 8
+ X BES S | »
: 4 co 4 =o | a
i fi -“ oo oo © =}
a 1 1 n un | [+]
[ J o |
. ]
-— [)
9
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
=
0
nw ~~ mw ®
- ~ ~ ™ < ~ -
! ? pra = i IE I a
ps ~ Ta a a Sw
F ~ od td mM NN © © 0 ’
~ oe Bis es BE BR ER 1 5
{ a 1 1 | | -
n
]
—-
Fel BN
[]
S »
od ~ PY - wn © ~
| |
Bo


N
Jd
i
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
o=0
with i=1, ..., 2.
|
7
\


N
A
J
1
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, ..-¢
I3
= o v
%i,0 = (-1) ) (“ay 5-0"
v=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. » >
g
|


N
A
J
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,
L+)
| 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.
| |
|
|
7


N
\
J
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).
; Summarizing, the procedure for increasing the step-
| 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.
\
2 |


\
A
J
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
’
&


N
A
J
¥ 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
vA
) 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.


N
A
J
§ 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
v=1
Note: &(u) satisfies
EE —
{ *These values are retained in a special work area for this
purpose.


\
A
J
}
| 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
- 5 y(t) is (k+l) times continuously differentiable for
2 ’ te tn, - (k=1)h, thy) has been made. Note: If the
| order was not initially increased then the step-size/order
selection logic and the interpolation procedure insure that ip >
the values A— and hy, are accurate, at least locally, >
to order k+l in h. Though no such statement can be made §
if the order was initially increased, the procedure, how-
| ever, does work satisfactorily.


\
A
J
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 '
<Q J
— 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
A
J
| 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)
| y
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.
Pa
tRecall: All three algorithms estimate the local error )
| using the norm as defined by (5.32). (F
|
| |
| |
-


N
J
| 114
i Ty A
B...... io een _
- f= | MOE
itn
vio = 3 |
i |
i |
1 |
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
|
v


N
|
\J
ee e——
[] SirTERENTIAL GCUATION worm
— EE —————————————————
p j= wwoy + oy Sp —
vien —.
- [ase 1
Parverm
B 2 PET
EA BE |
o J
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
2 |
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 |
+ of Jacobian Evaluations J 7 " | 1a 103 |e | 208
rinal Step-size (sec) 1.28 = 20° | 9m 51 | an a 3s 248
cpu Time (min * 107) 139 nm 52 | 200 | a2 | ss | ose
|
—_— 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?
+ of Time Steps 01 4s $30 | “Ss (13) . 1492
"of paction sisstions | 73 038 2a |e nas | - | en
peed ES FU ER SE EE
\ Final svep-size (sec) 2.96 » 10° | was os | 307 3s . | 220
\. ru rime (min * 107) wm wn PE | 202 m . i 4
Ph 4 —  oE JE EE S—— 4
“3 SOLUTION VIA DIFBOF .
> —— — - —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
+ of Jacobian Cvaluations | 42 5 | = 56 53 56 ss
Final Step-Size (sec) 11s» 20° | 2s = 20? | C1 | on ae 30 m 5
| SDITSUS aborted after © time steps indicating no solution vas Possible .
{ CF
Table 5.4: Problem 2
f
|
{
|
|
-
- \
Po


N
3
J
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
»
LY
o


N
\ J
ee ——————————————————————————
H Sree Ti se
- [ og ]
- fi = mp o Von ug
vies -t
13 gis =e © oeie we
go =| 25
i a J
ware
[es - *] |
sel mn 0 |
LY YN) |
So» |
|
ane
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
- of Time Steps 2s | aes a | sae | 109 " " |
I» of Punction Evaluations le | ss | 2¢ | un 1758 2046 2348
of Jacobian Tvaluations | 45 | a | es | & | a0 1s 134
inal Step-Size (sec) 12020” [1.200 [2s 519 | wos 244 1.2» 0"
- | { |
U Time (min + 107%) ” 5 12 | 1 199 | 250 266
| 0 A LE A
Viximum Global Lrror 2.32 = 10° 2.58 = 107? | 2.58 = 107° 3.99 = 2077 3 «10720 loses 00” T17 » 10°42
jo of rime steps | = ms | 23,55 | 26.000 23,602 1,00 23,997
| of Function Cvalustions | 568 “ | ea.700 | e000 6,203 1.802 0,907
) of Jacobian Cvaluations n n n |» 24 1u »
\ inal Step-Size (sec) mn | 193 | 6.30 = 2077 | €.19 = 207% | 0.29 = 207 Ja. 10? [3.00 5 207?
\ TE Ue 5 5 | sso | sae see seas 5728 v
- fA SOLUTION VIA DIF OF >
ei ! - - 3 - —" ™| SE
xisum Gloral Error 2.93 = 10 | 2.97 = Wo 2.9% = 10 21 = 6.02 » e620 2 12.1 = 0
of Time Steps m pre 266 ase sar 11,058 12,19
. of Function valuations | 402 sas m 0,947 2025 nm nan
of Jacobian Evalustions | 26 30 " 0 3 1 ’
inal Step-Size (sec) 388 1.2 = 20” 27s 0.245 261 LRU 10? 0.64» 2077 ’
I —
“maximum Step-Sise Permitted. /#
Table 5.6: Problem 4 =
\
>
v


N
J
J
eee
| Tree SATION rr
PO coven I—.
- fi = waugin) pin = vain
with ane -,t 1
! [.
P| wy |e cos wt
! | =n | rt
wn) PTR Pa pre
where | «2 |
i [=< 22 |
veyl aaa |
L2 2
|
[==] |
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
N \
- »


N
J
i |
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
§
. \


A)
A
J
i 120
.
o
© ~t
@|~ x
) [ [ Ll
ol 'ol ©
wv — - 4 i
| ”™
©
I=} al
wo | ~ * |
1 ) ool m | a
ol © . of 5
p= =| BE BE °
| - >
pt)
. i =
- | 1
i le | S
ol | ®
“lv ye |= |
x| of oo © =
1 m| oN A A °
[3]
-
—
= | pf
i |e E
of ~ Z
= EERE z
x ol of ©
Po = = I °
=
1 u
/ ; g
/ :
ol» -
Al
x| © A %
| Nl ~ ~- 4
| o
/ : |
[]
I]
— — <
v| © y-
ol © &
- e 2 El
3 — L
£
| | © 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
{
i |
|
v
1 |
|
|


A
J
] 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 ’
/
I C
*Por an unexplainable reason DIFSUB aborted when a local
| error bound of 10711 was requested.
\
-
|


\
A
J
\
122
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


\
<,
J
i
i 123
and
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
- proved stability properties of the new methods over the
backward differentiation formulas. The Jacobian matrix
| has a pair of complex eigenvalues at =o, + Ju, and a nega-
3 tive real eigenvalue at =O, where Oyr Wyo and 9, are as in J
Table 5.6. As the ratio 9,/9, = 100 the problem is stiff.
y Moreoever, since tan" (6, /0,) = 55° only those methods that i.
] ) are A(a)-stable for a > 55° are asymptotically stable for
all step-sizes. For this problem, and for the new methods
presented only the order 7 method is not asymptotically d
J——
( #The Lambda Loci for both of these algorithms are identical.
f
|


N
J
| |
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)
1
|
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 |


bY
J
. : .
125
i
ud
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 N
i = f
~ +
Ni. §
: I oN
/ a) Backward Differentiation Formulas
I P
” / / [
/ L
(
RAY DEFINED \.
8Y (549) »
| —— eal “N
i . jig "
| ==
—_—t ee eee PRE
| . I 7 0
; | —- { .
(TV | pr C
(
_.. <a b
{ Ni \ f
a?
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
\
a


bY
U
|
| 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.
\
w
4
= ~ |


1
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
; CT
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.
| | |
f{
| |


N
A
J
128
3. 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.

?
i \

»


A)
N
|
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.
|
LO
{
{
|
|


N
J
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
H
\
»


N
3
J
131
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

| |
y \


N
J
| 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.
\
|
J
|
|
(
|
I :
|


N
A
J
|
l
|
l
APPENDIX A
- COMPUTER PROGRAM LISTING FOR THE
BASIC PLOTTING PROGRAM
| This appendix contains a computer program listing for
the basic plotting program described in Section 3.5. The
| particular version of subcostine ACCEPT listed reads the
characteristic equation for a given method directly from
cards. The input card format is described with the listing
| of ACCEPT (deck identification "ACPCD") .
3 The program was written for execution on an XDS
SIGMA-5 computer. The program is in non-standard FORTRAN
making full use of the power of XDS EXTENDED FORTRAN-IV. 2
In addition, it was compiled using the "ADP" option wherein
all variables that would normally be typed real are auto-
matically considered double precision.* Similarly, vari-
| ables that are declared complex are automatically considered
( double precision. Under the "ADP" option, the compiler
[Y 4 automatically sets up the calls to the double precision
4 | FORTRAN intrinsic functions rather than the single precision A
intrinsics. As an example, calls to the double precision
| sine function can be coded as either SIN or DSIN.
5
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


\
J
|
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.
|
|
|
|
[
I
|
i
"Ye
(
(
|
|


N
)
WV
|
t .
135
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
ti
{ I
cps
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
< Yy .&5 2 H > 5 =
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
- $S5p “832 (7 IRR SSRI Ei FRE. ppiiss § dss
oe SpE. EEE 53 ITs PElppisiiy 2. DIS 5 08 CL 4 =
s 282353038 32 Spo Lob & wee oss 28 2:
shat icEpals $39 IEEE edleg
$s prisiziiii s2z $3353 F Sg: fi
SIMI IEL FEF BBE 2 OF RE f.....ceeesenned FE
| vou vv we wai —
\
.


N
A
136
srsescassrcisesasasiiRtisifefEbiRisiIIRIIIzREbREbRRERRENNE
pryzarssryrzare srry ERZEZr EEL ii ln srrrarrIEIEIZIIZ ERLEREL
CE EEE EEE Cr EEE EEE EE EE EE Eh
- " : i § i
§ ye i A {
i Fi E £4 ; :
Poel £E Ef F
g Pat es #3 & i
- Pog Egg 3 3
i ii t 38 2 : i
. x PE tf: BR OE OG
g 5 Ts § 2 $ 2: of EB £
5 4 PEE * bo 81 st?
s H Hd - i * 3 pr oe 5 2 =
§ £; :i = ¢ 0 ogi BORE Be
2 sib f fs EY GE os: § $ Zef fo +
» Sof PREY EF Te |. capes. § 23-20
H & a-g $5783 TE. Ea FL gp E%Te- : 5 32 33:
§ ob: fetd BE OES. pT Waitin: fpf die
{Bef FREE NEEL fe RRHESiD CCE fig:
on TTF £% D0 Sr Pe 5 Te *e¥Toetfg FILLET JSR I 31 Th)
7% Bt Bi BS I ETI i | EE i Sat
; al A De pe es PU EHEC
: ooEss tf IERETCIRT BITSS.ilESoftE. FIpRisg: sali x
4 Fig, 8 =%%=3 wesEoEe Tors es SETISRIT S355 &
“ £2 Hl E44 £2 2 822 ©
po wo we” - - _— re > — ey"
| messraIRinIsERanaRRsARSTEs ST STEER REE 0233332RE02E
BE eitissosoorpzsssseasprss snsssspasss F220LLEsezszest
EE EEE EE CE EE
i TE
: £ : 3
Fi § 3 H : 5 &
HN LF TR 2 B - p
| : 3: Ef 1} z i
| § s=& ¢e ES £ pd
FH ef BE 3
‘Hh mz €F 0 OH: o- 3 HE
w:. £1 1§ ss #1 > § $4
ern bt 5% ESF 4 ° i § 2s
Eseple i §g 13% 8 : OF $82
if gasffes pv INI: OL jit .
Sid: fy ME 3: Oot +38
HES HEE TEE A EE 58:
sp iR Ee ppkygd Bd; FEL (PDF
EB -is 88 3 oo 85.8 5: SB- O31 wee 52% AEE ER]
2 §:-82 §2 sr B95 ¥:§ EP FR Esl BEL sic "EE } N
SETI 3 RS ERP EE Ls BSR EE mcf :
ER Tl RE fie
Rb TN ER BOR PN Lo Cl Ne ES PE
8 EERlRTELT sts Eke Sx t © Eee sige Fod BiEgl
: =f NH fis 22 3 fF fa 2 HiE
| EER EEO Lt 8 8s WB #
[
i
\
.


N
EY
%.
J
| 137
F333
zzrzzzrrzzr
| [FEFTTELT  P—
i § $
3 ry
Et § |
} ital |
TE
= ER
| ae
: Bz 4
: 2.514
: §2 4
| HRT
BE ol |
Post
H ty
| Pk
Pb
| PoAEst
3 235%LS
HEE
§ 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
- sess: : i
: = : § . 2 FA
[431 | 2s 31 z : H : LH 3
gif oi : $$: 3 file: 1
i a T i! PP 5 Bcd 8° ¢
Syy 37.01 $e 5 P28 4 Hb i
2 Bit i: idi: RE i
sis £=ih sf fefdpcl it tis = :
| £2 zgeei zg 28 Re 5 :
Gigi ii, iti i: :
HE RIE itp BE Git ti ;
| oft Siri i3 2058857 2 seid ht 2 |
sit Spyi 2 Eo pe=oatt 2 A
i fp Syixi Ew AN EE PA ::
sf sel 3: Haas : Basin g: ? :
2 EE ERC Ff He ;
ip ric Rae A FEE EEE 52833 :
Hig PIERRE iz TE Ra Rl] Zot :
if lize i iE ede a3hf - ro 3
. H TET H yEipppots E3¥y 2 3: ;
et prfsgsialt " 1H Lo
| ra" 8s 2 £23 Ei Bodies. ;
ones grasps phos pe
Nose ; be
sEEnIERERENE ————
: sEIizissiiiaEiiil ———
; FRR _—-
H ETA LL ob
: H 5 PATIL _—_—_
| we PASE
PO 3 EO pre
{ HA Ea Te 61% 8 -
ix i EC F £3. o ies :
ig: £7 ° §£§ En eye S32 83 5, 5 i
’ i854 S¢ : BB Es fie $f 23p § :
| dps # £3 iz fj=r5% seacyv if ri ;
zg if BE £ s.. Cif 3tssifs fro¥s cfs Bi EE :
i Bak ue Br Lik Benet us HE
if Sohal Gd feiss £3 AG ; i
; iE Jeph hue sir seefial ions go
§ ic epet ist se §isbgmi- tiie HOHE i 3
i iF HEN ni Zejcanl ER TT ie
ig petiz fe §2= Zaiss, pe BRalt,e Lt ] |
; ; Br tts iE EaiitEl I CE TAT =. 5
| i sides ot: ini grit ILE Lech EES 3:
* 35 Sia fs gi*- fe Andean fe3- 4
Le i Hy HE 58. oF
i 3 Ef LA rey 2
yr id Hs He plicit wif |
ZI: i=% Bg 7 e753 Cine gig:
{ gif i5: siz: REEL Ca] 1 bt
$3 get ~~ Seg= IE HTH EE ;
==: ;
S— a
O— | :
wu
—— :
NO— :
.


.
A
J
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
p
| i
B ; :
| £ a 3 € =
| is TN : : :
: 2 Eo fs $5 P
- 1 | $4 =" a re Fe
J a6, 2 3ba_gz 3% _ 34. f Ei ES
: awzstets Sfgan T 7.380 8 EET = 7.0TES
| Y iiisssz eo i 5 S0583cf $8isaitas stoedad $583.00
FIR EE BIE BI Fr Er ive Cy EEE ES
§ Filles d sos Palipiigtiiies TE EEL et i]
o3 aaEIZIZS Le we GRSLCEEESSSEE eo TL IST eTIiIDEZI:E:
: TR prea Js SL LMF Pd wut) 0 4 CTLPPe LL
PI LE EE PEE lh ot Wowk Sait
H 3 oa 333Eeeziinaziil 3335s SSS EES CF
a2 22 E ss g 2 B :
LEE $= 8 83 gE 5 .
| sssssssapstzzsaasiBingaiIIziIIiIIEISAIALIIEEIIIRISNCEC
SessriirrrrrssEiiiiiiiisiEiSiErEiiETIEEEEEEITEEERsEEEsnnS
[ EE EC EEE EE
tH
3
r
| :
H
:
:
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 ¢
|
| >
4 4
’ \
.


N
I
J
| 140
rib aEREsRARIASEsR RRs ERRARRRE RRR ARAL
PT 0. .
a
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
EL
| 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
7
|
.


N
A
|
141
|
£35505333598KERRAR SRERUSACITILILGEILINLSILESILAALRE ZI
; (33334333333 5555555 EEEEs5ELEnE EEE EERE E5E5 EER ER050
| 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 ,
; fp. c MEEtcpsBpesp i 2 BF FETE
} 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
| |
i
|
\
8


.
A
3
J
} 3
142
carerersenosrsrane tsa LITIRALAREEIRARARITIIININIIEIIC
H 3 H
- 2 «sr
it z £ EP gs:
Bh 5 ¢ ¢ 3g
# 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
TT
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
- Romo 2g oR * fi: iH
| § t § . § I ge + £ § C82 :
H :' § 8 8 § gE . & 3 ¢* :
$s = - . 2 3 2 3% 4 » tT 22 H
fife alg 238% = ¢ FF aE 1}
p $:f % 8 TEg get 37 ER si ‘
: | Pi FE 8 8. 3: pce piiopH if :
IR IEEE EINE TT
$8 CS Re 03; (EY oppcn SRC aE
fo fish 3 ¢ TER CBs gl PF ogo 34F C5
Ef tags 80 cif Piel Ba yond {EH .
§ 0 TeSSEiE az 32 23% I Daely aati fiz Roh ,
iz EAE see 333:% 7 I. FCF ep 288i f.%ES gs ili
FEC US Sh Fr EP ri Ho PL Le M4
Ea TEEN Te eR TL LR da I
os orniness gizizie=Sizzil ld goosiisgd at s823lig
$s 3gzz 23 ma Fz Rozz: £23 |
“ - ws |
X \
.


N
A
 ) Y
{
i .
. 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
\
v


3
A
3
3
| - 144
é&
sgzszaiizis
EEEEErETIIIG
i ETHEY
{1171
4 ; H
¥ i
| :
H H
H 3
| 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
| y
4 '
| 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 |


N
A
J

| 145
2595853055 SERRE AL ARNLESALLLT
LEETIESLILSLELEELIEEELEE IEE .
SEESSLEELELLABELLLEROCLLLLY

| SEERRREIEENRERNR SREY IIY SREY

$ $
: i
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
. 2
bo.
.


N
J
)
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] |
bi
2 \
v


\
A
<A
J
147
y
0 ~
«=
i $8
82
® 0 & ~
. Lodi
i E23<
us 8 =
| we va
ale # @
1 | EI
. =
| TE 5
| La [3]
| Q
| :
| 3
| El
Eo
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 ’
How
cg © 0®
- a
9 @
a Ed
$85
go 8
©
v


J
|
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.
i .
H
|
|
} .
I
|
y { /
’


N
A
N
J
’ .
| 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
+ 2% $3: : FH Ez Ff
i & B CR | I
£3 gis: : £3 : Lot BHO
: #3, 28:3 i ges = & o° 8 §. ;
es 28552 PRE EF omz® 3 EF oE.T oeloc
i SI ER tions fof HERES Pf Spm B32
g c8icfr Siz 8.8 Bo BE cE id Spc £3 2
Epest-azifss : fd S230 CG pb oposef Zc
cp Ren phy of Noe eftgzal COEEC E70
- PE2TE FS: TL Vo 2.82 E330cs 5 % »3Ei% 837 3
£ = BIE 832 Tas Yo, SEL, Bosuiel.c BR BEE: BSS
T PBlizessh oh, 533 Silss E13785.5.% ER PRIFI OBEY
$18 EE hh AE Biz fas, giiic
es pREZLTIIT FC SiV oe 55." SeU2P.FLtE Se Sh-ce ETE
It FE het P3.553%" EiT= §23-~F~F-2 I &:CC0 ee.
i Se SErgesefe ,EPIPERRSE LUIND Efisspepal IT LTIERE LLL0S
zesrsg
LEI, me oe el LEE Lk mmam fT m3 52
i camemeneroon zens ne sgIaC RR ELARRRIRARARS SY TITIILIIIRILIT :
prrzzzIzzsrzazziiiizairIriiIIiIziiiIiIiifRifffffifil
BEIEEIEIEILI Loon Eton oObt obi bEbRatEiE Eee CURE RARIRINIE
get EEEEEEE tt EERE BEERER TURNER RRTIR RRR RRR RRRR RR Ras
:
¥ 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
‘
\
v


N
\
J
\
| 150
ERE asITIs ILE sRssI EARL EIIILsIRI RE 2
1 LE EH
i ;
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
- we wu ww w=" wt uut uu” wu”
| Bee fEACITITITSsaEaEInALEIITITITIINRAIEIIIINIiIIZIIIILIiELL
HHH HH EE HE ER RE
| H
bo
§ 1
x ;
| -
: ” e
z § 2
Ss 3 § cE 3
- ~ °c ev 3d
¥ S 8 5% . 2
cfg B os 7:
ef =z &. SB 3 et £i=-2 e ¢
Bz == 2 2 ESE £Té 2s * 3 .
2 FE 3X uf OE g2 £3 I xo § ‘
85 of E 5: ec FE af llr Coefe.tof :
- = Te eo = - a =
i § 0 om 3 osriele arses Nts
“ SE_ = = HES RTE EE TEE
§30 E-aic E st fr eS iF aii led Ep.
£00 2.0 a2 2% 23 gZ.d Biov,3T Bede Co Bfese o 3ed
I , 59% ifs 37 isis Ft TE NE NR |
_ EEC EIS BIC ESEdsi PIARE QRAGRES PRA 3 eIThS WITH
EAH 2 522 Tamas forte Marlins Veolia Turion iielr
© oR ip 0p HED BED BELMD BE RT IRD
: BEal- a . A § .
HE
|


N
A
\/
\
t :
| 151
S853552TS5 SLE AR LARTER LRILE
HHT TTT fr
Sos o008s0%59%% RESERREEEALE ey 333 3 18
; FEETEETIEEEEEE BEA ARRRERRRALERIRIRLERALIALE ELTRARARALL
% | 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
- RE 13 dade 10H §-1% 21 2 3 IEEE : < :
H faessiivis §zii 31:8 °F C7 gs; £ A
[| aE AE EC
RE ET fata gio AE : :
FREeRVLILLT I SETIETIE 3 83: ia® £3E 7 2 3 .
3 Bgaekifioiiz 20 333% JEPS HL FE - JE
| g Fe Che $5. 3% rl 's Ta
Mg i Re Ye $232 °F dEaizi. o if fs
§ SXireesiey Zoiiiiin Baily RET IE
re CT Hs iE a EEN ~o§ 23 ’
org asters rill git iE $aoofiel § ES -Siise
ei teesier p °° SeNESIE 22% FEET EH LN s “¥s% Sve
E> er & FooTis IT Soo.inieR E +1 YF Sep +
i v (253 ffi: Uf TiS T Fei o.gR
[ 383 ERE I TE NE
| Ls dana . .
:
i
|
|
.


N
A
O
\
\|
| 152
|
|
| =ofsEarIiTEsRRafanafusITiTITIISSNAsiiiiiiiiiiiiic
ERERRRERRiARKARARLARARBRLEREEAREERERDERADIEEETEIEES
SP ELUPELEEELILLTEEEES ESSERE IEEIERIR IRIIINIIIIINNST
| :
ii |
s £2 2%3 2%
g i = £3: €:
- $= $ ¥ eo fs:
| Ta - KN - = Teuess
35 Fa ef f:0558
z 33: 3 gE . 8 sz =3%033: ‘
y pl | 2s 8 EF 2 3 88 258-5: }
E | 2 223 £? * XX 3: 3 $E zEITss:
tr cE & FB 2 Te E&. 8222s =
FR A FER 5
EE: Eis.) fd 3c ilcuEr ©
+ ss: § 2S8EVE. es 28-72 2 2. SEescT. 1
frie: SgtEat FOF (SESS £3 EEL. 3 » >
bod a & Tires * < Fe.52 ¢ $s = REE 1 Fe
aT 5 ZT LPLLaB- E Yo. zPESY JY Eo. | 1 °o®
$F 3 sFsesfc BE Zoo EERE OT 2%.5- § ge
re 2 E2ac at 20% 223% goseatyi3lii Twi 3 § 5S E14
| 22 2 iigssgeg BRE RESATERRIL ETS
{ Se oIVErEre ny Eogel TRS ipse Siyiis 3 3 sac” 33
EF fees EErfzic oagE> eRlgiICiR Vier FF FtEpaii.L
IT Seedpiuinie Er YR $loptae®® SSF tr famziil
ss §tt sg 8B 3 £83 8 3: :E
| £2 gz¢_ & °° = s2 £ 2% 2 |
v


N
A
%,
J
\ a
153
. 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
£
{ i
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 |
|
R |
.


N
A
J
\
| 154
)
3353333835358 RRAR AEE
| 333333332323333233332
{ gegEsesesEaiiiiitt
8 3
| 3 :
s 5 .
| Ef
H :
: - ”
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
2 a
: 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
* & s == wr i3dy Inesssl 3 22 23I¥8C A3TLac = = 238 Se
rw $5 [798 TEIRTT pIiziet if Te BSE "ER 3.03 a2"2¢ 05% So
£og-Bf  gaets. erit.d 2s E3.oTpEIiEOCR] BROR. 3235 2
SEP: Feta pptuiBE Te to.cnifdZ.i.iY Ogle: ST
{5 FE33% Bia8isz oF TAs EEiizzii Dil EISSB
2ifFLif 2k x Figiz: ex 5
ddd - Cr sez J. : J.
7 |
i |
v


\
\
3
J
\ h
|
| 155
PTT IT I CL EE A LER SA A AA
ea
| 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%:
- &§ 3 2 2 8 2 28 ¢ 8:
we wv we © - ". x al w” nd - - Sou®®y
} ceremenass onsen et Ta KARE RRAREARARERLINIICINITRINI INIT
grpsrssszszsssasssssTrTsEIiiBSSELIITTInIiitifffifiii
BE LE LER ERR Rass ra pas pba BR LRS
i ; :
9 ¥ sss 5 2
VE Aa 222 EE i:
i jer oH Oi EF rd §j zi
geez ly 3 2 $e 5 £ 33 58. oF
E308 Bec iit i, iB gEEEfZ
3 REE EE EL ik ces gz FUE
Es yy E55 37 2% 22 Bal 222 &§ % S572.:0
re ? Sel fc ¢ 1 55 . ro. ses * 2 ses-=%
TEI. ¥Sei¥iimes 3 nz Da
sB.5T.33¢ 3832 i. 338 ¥ 28% §. ses § E ZEEit
c OFf REGERS SRC EES 5 RT. SY EE 7 OC 1334 p
£ §3% == 1°, 2202 = 1%7z: ¥ 25g: 85% © =r TEp.ef .
ERLE I ERR ET SE Yb J 8 FA:
. EERE EMRE elonag vp Rooofp Init gap PARR
§ 3% epEETUREE mize B00 Broaf Banal BD Gana
| :  BRithaeatgise UPS Patient fi imi
I FE 992.85°580E ¢ $a ERU3TTIREEE SED ZRIPECS 0
: IAAL TI HF En e2 pang iziZZZ Eos pLTILLL
ocEEsegipeitHipl rei. Fag pienersTTT BED OTC
y § SLPS, Eien CRE 30 BTC AWERNNSY sel HHL O2
Qa aie al § opee®5.2f 33E iE 8
H Z PFEz28r3 r: 32 2
"oe ¢ seLiswist gs 8°
| |
J
|
| I
v


A
\
J
\
| 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
PR
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 .
- vi et 23852. 5C-8% "e720 ¢ 3 » - on
| : teuEl BE. § ETBSI35ITE 3B. S9%: .: $3 i
¢ geet Poi & RapPsifsalofiesili: X - 5 RG
x IRPE33 gE | 5 23SEee3S FS. TESTO 38. pf 223
v si psaeize §of Ze7E2i800oRl5pT CRUEL aR Ps E.I.F ‘
| EHS fe LAr Tt BIsMay: Bo. 3z332F .
4 ELE oe SR Seri LpT 582 81 SEL
Ep raiip".® 2Eai. 050s 22,8852 30 Bagel
ceils ge yTorilz Sls" “iz cE 5%
| 3 ¥ OS=2 8 & ES =28%
[ rd
‘ - -
| | £ 22 2:2 zg 2 5 2% : 3
& $8 2:3 Eg 8 23 : 3 |
ou -
v


\
A
- \
J
\
| 157
SE 3p3LEsTEE
dH
EEEEEEEEEEEE
4 PRaEEEEERREE
Y :
:
| :
H
x 3
: i
Eo: :
} £
4
£2
“Sa
Bn
pe H
i I
| 3
§ §
I
i. ¢ 8
| Syi-a, ¢
A
Ligetisiss
2 ERarirs,
i £2%e
es?
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
; HH “a - Sof S_. wis °F e 8 ToT.ER - .
: 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
- Eh tI FN S58 Yer <n yo Nd Po £ ec zz
in BEE EE: IB PRE sgls Tn CARR
$28 3S Tal Bi Sap STIS TIT paeft. Eigld
! $3 TEED LER EER TRE Tey ofa PRE §iBNeC
sir 33.33 gid sil ial Boi: LaEED Ana
- Tei. ErivIesy rae TT N0 220 Ernie asT. FET pe Sy
ié-0 HM SET IR Ee Testo £3ERCE SEPr¥e « 2S3ge%
fefop $4osERRERi CRRRTigRaREc oT BollSt fisli: cin
Tele JEezebes rR, LT.ETE20 E228 TTT ae ERLVLT JTS
DI aa PEON HP Ei
i a: % ce 2ge 538 ge H
| | 3B $3 33: EH gz : |
| v —— —
v


\
J
\
E 158
careronsnpzss seine eLEtR2TRtE
PE PTR i
nh RS REESE YY
| fELEE feed ee LeeEeaReeiaeass
| -
F
3
3
=
ot ,
&
a : &
E * iE 1 .
£E £ EE 2 eo
se § &% EF 2:3
* Pep PF oc:z:
¢ = BE BE xs
s B® _o2%F tz tT: =:
Ed sb%® E% SEE
: $i: £- Z2%3
Popiied 3: oii
 «:t&P87 3 23
? ei%8Fz Titi:
S sZy®: SET-E
: y=¥i3c.tiii:
- $ “eo. 58 7 32 52.
§ooeprkfs 72eciey’
2 oa53t § 2¥ 3-i%rEiEi.
§ ge p ERR
EE ge § Fi-BEesRated.
2 E & & Poesia cicnfot
{ is 5 A 3
DT dd -
[ i camemenseoonninenntIRRLARIRARRALR
go85e908 $55555500000E05000080008008
§5500808 228002008 sunray ssaRntit
| T
¥
H
£
£
| es
<
Pe
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 <
< £257 : . LT - H >
8375 = I 5 Ei Ee
; Siz: & gz _§
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
|
|
)
v


\
A
a.
\
APPENDIX C
|
Y ’ FLOW CHART AND COMPUTER
- ' PROGRAM LISTING FOR DIFJMT
§ A detailed flow chart for DIFJMT is shown in Figure |
C.l. Unlike the other computer programs listed in Appen-
i dices A and B, DIFJMT is written in standard FORTRAN-IV.
i A computer program listing for DIFJMT can also be found in
this appendix.
;
|
|


\
\
J
\
160
; 2 Ee
- Ha — NY rT)
1 |: 2 FAM 37
i HE ENTE | 2%
1 — ——
38: nS ol vy BE
th ANCE CIE: . ¢
is VA AENAYA €:e
if VARS 8 J) C1753 J La $37, 5
ii p mcm— [2 £3.35 [1
i 21: 1 K y 2:75 2 we [=]
CO FOBN PRE
VN |; E-REVER get =Vel 1% Did a
| Fi * 5.082 ]
H £5352 °
LI 18372", 4
3
31 @
I TS ——
H MH .
vo|EE =
| —
|
| = g
i | {22 | El
id | 3 z
133 | ir —
- FH. rl zo
: B= Sw
HE =f iy | [23] Ne 0
= SE led | [ETNA po
i i fe] BH HAE °
H Ld
4 : 3
§\2 P1124 x 2 _® * se
} EET A EE dE | 82
als3331% = " |
oJ cls33E:d | OY
: P
— 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
ki!


\
\
U
\
161
¥
is
| 5 0 : .
) . ) EES. 2 kB o%
¥ i kid
| A en [A
NAOH
| nai :
INE /AANEL dL -]
Sali :
a :
-
— a
od 5 3
$7 o
| - = 5
| ta :33| [33 u
: G23] |igcE
11 (13 -
=7 MoE [1
i hf Ered be 2
He Hod | J V2) Beal) EET 3)
po Woe er 5
i AWE 1 a (7% 3 NAL \» LJ | ELF Sw
H RISER BER TALE ad 1-351 N/M 4 = O
§ E} 20
rE
=
mr 8a
i — IF xa £:4: ALGER) -
i Ri TEE [3 fit] rErEEiiE =
; +5 | [428% i) ii dn [af EEE -
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
|


\
J
h ’
162
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]
; -— : KH PB
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
- [s22:58E
gdiiss ce
|


N
A
UJ
\
| 163
ssersgiten at RRsR RR EReTataTIRILEIEI FARSI ISLICIRAREETEE2
REEL EEEEREEL ERLE RREEREEEERAEALEE ELE] LLEEEEREENEEIRAREREE
$335335033504 5009040400 F SAF R ERE Tana FREE 7313225053220 9%;
- M rgdehdar : £ :
4 gi. sg § fGpiipgoc. fet §.: | 4
y 55, o, PE § e271 bey CEE iE. XE H
LE BCE A FL Los y 5957: § 3.5002 2 »
Six wht Zee § 5° % zeit Jz. .% iv et Sr x
| 2o383,%05 £3 8 Fe oR aE: si=3E 8 £8 £2 7 52% OO
Bipmevs? 3a] PZ sf 2 BEiESrny §.708 1 yop FF CEES 33
‘1 frgiist roils 52 2g@ romidyeesii (i%eg 1 E53 ND 138.5 8
WEepZia SE 5 iB 2 SoRePTS FEa2E © LAL en 88T TO
BE aE Ife _ 3 caz © FE esEz? £53EE | £25 LEED 5%
. PETES T LT 5 T INE PF OEE. .5Fs E2055 5 2% 23L3T 2% 03 .
: 2. SEikle ® 23 wma & CEE .8F oss F © vo” B2e06t ww CF
SEE Ess in $73 58 uw pfiioclil 35:7 © dU oSeSiiS :
Eb Eiveien; 30 277 8 R23.EEC,3 008 oie SE fT ES
| PIPFLLG YL arf Rib liEoreb iE ETE 3 151 A
Baif.psiigBfC 23 gpdag oi i. aelnrc 255 ontbo PLY of
FRElefroizip, = FS EERaY Lael, SeiETOOTT LoEPp ETpollRE BS
RI EE NE Et tT CH ta Loy :
Sircizyiale. gucte Blu LoB:divitETET $ 5-7 Titles. of
YEEIET yiifis 0B. .sECES geEESeE LT tr TEE aDSTRTAE 38
| Fp aod TE TEMA i LT THI LE iiTrasii fi ly
Brians o B2T8yi%% meoPi EERTCEES ToElSLlIsFriiliiE
Sreerigogrmird © Nad To To i Fd ET ETM 7
gE, ye IisiE, "0 = nS ET ea SLOT pT, (RE SEER E00 Eas
LL SER Ll Mr Sr EE FE Loa
EEE iT Cr EAT TL PIA ME TTL EM Tl
Ett ee i AIAN LEE Ha 2 CEE
| | > Lo &a i it £5
03502 LE8508 S00 LE SESE E0880 LEELEESIILLLITELIIITLLLLE
TELAT LUCA LALLA ERELEE
225. 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:
- t wede® 32 & Ze g¥ecs of ¥22f.i0
§ §  Gesps ei. 52 5% RT wy ¥R.LEU PEF geSEliy
£3 oaaleiTss $558 Fr Zs Daf dg 2OFCYE LL 2Es.elSB
T 3 EEgEINE TERE OF RF «fp 0B 3nT.8 iy dTiSier
0% pePitiis $Fos 3x sx 20 i0f spit 22% TY
Poi egTioees Eye.t zg 2% Fzbeel-kriifyi;.- EUEL%eS
| 2 1S gro-5 09 IY (SSLEFS MITEL yaiizs y
£5 i gmail MIE OG TTL MEM Le
ITIL ssh. Rroli. J{ rte NR RE
sf fgisBesil ff gdel er od custi pEpaitbots, Tig.ili3
of: Isvieiel 3 83037 £2 ao BEreEiilReciiLiT gli;
§7 § S3afFafz 7 alii BE 2 S5rroliepRiie dit PREilt,
" HE tT SIRE J BLE LH IE PT RE ‘
i i re §gabootf F ofEY ag. BF ly PITRRETECS.C RELCCMEL
5% 3 wEZ258 § Beit” 287 fo wifi an3>obiOiisl BE. -2L .
SE i ioszesf § Boll RT Dp mil tiliERILCOEL onrEiii
PouptoiEt 5 ZBRLLGEL. SF TolifnoecESislelr 20. ER03
BER SHARP COR sion Be RIEt an LE Tonite
38 1 Tomah # hEcETeogE HE Ft MEP 11
iE § ESSly oocgttinitiEs 2°. PRainspagiie; B05 Dla)
$B § EfSs..or¥ oppfelilatE SFY BIO 33EINEIEIIvE® Yodedt
3 SE3T3% BC LotodsaEs BE: 802230 ,838530] FePRYLE
TF fsa, ROloEiayis JB ERSEIC.Gi3EEL Vr SREICE
§ 0% 5Eviaiic Tet eT Dis $,oeiyiiieiuSesly IIL. EY ’
£3 RRe3drd good BItieds DRL EEE - | pipiic
£3 JRESIEED RISES ETogell Lil iE - = E'Eyr.t
v § E2387.cF LepaEdcipt IEE gElY : I Eyiiys
y §RERieas BERIII..CEp get #235 13 E FF
T 3 gERliee STERUIERLC0E Ro. Lilie 22 F og
£3 oibizfr, TORTREEeutigl e233 3l.v 8 : g <
3 3 oEepiiler EUB38s%%y GuT ric: 2 8 H y 1
i a AR : + 7 .
:E 1! JP RT TT RRR CREE FERRE ELLER EEE EE EEE eee
: a
.


\
J
\
'
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 °°
- £27 28335: iat, 038 "3. $2
| ceseessssBosessseceses fRoyseessessa, ,, sFEINEE 23°. ©
| ESO mARSY SESE Sa RanSs aS OAR SE SASS SAR rane ar REE aT ey
tal LA EAE RT FLA SAN SS REE 3 pt
a do oo aa ff of eo oS 4
A ERR Lk R
2555523000027 5345005005028 §400050 2050028 FasL shuns Eager oon
FF Et § : 3 gf :
sz fy. fi $ ? 3 ow bz.
& gs 5 8 2 2: § 3° LB % SE £
oe PF LEgt i: fF I» 32 22:8 §
g: srg: FEI BE 1 Ha ERE
$e 23555: . $53. % if , #2, 828 ¢ 3
8 gp 23% zz iwi. i ce S327 3
$8: 85,85: $ 3p: in CS o58 35 2 £2 4 #
TIF (28 p2d: § i883 31% Seeds ®
S33. fe F:3 5 385.2 iD ess 2% F 2 22.8 '
$ Pi® Fz 2% oo: .%E5r i: $3585 #7 $ : 9st
+s Seo-a% = 28°88: HH 3 4] a 2 E23
rank the BE: £: Fait: tot: Th.5° .. § § 383
§F30F (EERE TOLER ites R387: £8 5 : gris
3 -¥: Fi 57 iris §.65 segcs RF T § Tei
c sey s08c o£ 3 BF 3f.28% $e Seti: at : %B3f
22:8 ECE or oo: LE ERLE HE =. 75te 8F 2 § <P5%
8 EER SFL F 3 ck RFECCE iE S§.7.8 ‘t % ¢ .Bis
3% Lier 5 EF 1 3% iEEaai.y EDS g22Fa8 .% = § SES y
voy. Ee 2 3 3. 0LxZiss iis 2¢.28 2. = § PPS.
© o¥e© Zp¥. BaF: fa VRIsifc iF, Se:€% $2 5 3 £23:
§ Ja TEx 35%: $s igift 3d Soi 28 7 3 %id ,
% J aF3 SEIT UST 2 IT L87T3sE fel ST.IVE SF 4% oNtE
¥ Ef E°3 IS.y 335 § SU ..eciRe RES 2.57.8 28% F 8 ‘3
FF TEC (ed) 15. 3 oa iE%eriac x.y P8878 Y. ¢ 3 gg:
BF o5- 85s “Ex 3 SE SEINLiE” 3d 'ES7SS .2 ¢ § 0833
FhoaSead Sk (5% 3 38 320% to costs. 2f 22% § SN3E
EVIECETRSL 200 § FTL fivpi.. iPox 22352232 $8 5 3 SS3%
siEFEoisE BIZ DI OSI5 ironniic igs SFFETIST RF 2.3 “02%
3 ¥3252,25% £25 § O05 paEpict ak COLETTI 33 TR% ..3%
epdotasuri, 8. 3 EB. 3025%8 $= 1G a.35. B5 °32 223°
TD en.nUPEay SUn 3 8. BPETI.. IIDC wp3T03375 38 :TixBE3R
TESliIETET BT, 82 32% De.nitT ERT JTe'ETT.E Ct Beisnhiic
LO-EITTEY, PTE 238 fEEEiL. it 25.65.22 .. 283%:338 ’
gE ETEZEF, 0 TRIE S32 HITITiac IM. WS ELSES £5 S;EIVIiE
SUISESRCILLEISLEY SAX ix 71 HAL PAF LSA ES LSA SE PY = 1
> 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 |
7
v


\
A
J
\
165
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
i
sz
-
g rg
2 F bE |
H 2 st
- : .
£ - : « #
113 : : : §
ei- © - Fa
5:5 . ¥ . -
| 4d 1 ¥ : g : §
iis : : FEPELE:
. i: sg ; e228 gg} TS
: . 232: £ 7 82.37 sf 3°
3 2:92 eset o z - s82i: 32 Lo
& 3 o O33 2 § £ _P°.z® z255. °° y
’ Beef 2% o EEPESE. = oo T GERISiion.RiTR. 8
¢ 3s $5) S23%ESiE Of ME EP BREICIELN-c cE
z 8 wis EE oo rr ITY Ceus” SaS2iss’°SC3525:C §
Z§ uff: BE 3... --9ShLiBER ZEBAcr FITIIECERCT..eiE
| 3.2  LeEl BF TRIepsiiaSiiiee goisi. ToBUTIT $ ooo ©
Tl TF Ps PSE PHD bd Sd TIT Ml Leh
gies o3tsycoiiZgE.SEiiese co. BE. b
i :  TakERsmR s3se 533
- I. - -
PO TTL I Tr Tt Fr Fit HT S43 4
HE a
RLLLEL5ERA50LERLLEE baELARREERERRERRAERR RAN AAAMAA AS A 0
SEadRad Raa ial Raa 2030000500000 40500065 5500 0h0000eEE GEBEEE
{ i
14
1 :
s i
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 “ |
i
.


\
J
\
)
N
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
- "
5. 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
4 .
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 |
7
\


\
J
\
' \
i v
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
- H
| £ 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 |
v


N
A
J
\
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.
1
|
|
|


N
A
J
\
| 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-
- sertation, Department of Computer Science,
H University of Toronto, Technical Report No. 46,
{ October 1972.
i [GAN60] F. R. Gantmacher, The Theory of Matrices,
Volume 1, Chelsea Publishing Company, New
York, 1960. d
| {GEAG7] C. W. Gear, Numerical Integration of Stiff
1 Ordinary pifferential Equations, Department of :
Computer Science, University of Illinois at
| Urbana-Champaign, Report No. 221, January 1967.
i [GE268] , "The Automatic Integration of Stiff ‘
3 | Ordinary Differential Equations”, Proceedings .
IFIP Congress - 1968, Booklet A (Mathematics),
August 1968, pp ABl-A8S5.
i [GEA71a]) , "DIFSUB for Solution of Ordinary Dif- 5 >
ferential Equations", Algorithm 407, CACM, 14, ’
I 1971, pp 185-190.
[GEA71Db] , Numerical Initial Value Problems in
Ordinary Differential Equations, Prentice-Hall,
| | Tnc., Englewood Cliffs, New Jersey, 1971. |


N
A
J
\
170
[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.
1
| ; |


N
A
U
\
| 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.
a |
;
|


»
A
J
" \
| 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 ‘
i .
)
|