ORNL-TM-1933.txt

 

- -'1-. i.i i"

MASTE
OAK RIDGE NATIOEAS? LQQH'SRY

operated by

UNION CARBIDE CORPORATION
NUCLEAR DIVISION
for the
U.S. ATOMIC ENERGY COMMISSION

 

¥

ORNL- TM- 1933 7 = |

COPY ND-i 33

DATE - August 30, 1967

"MATEXP," A GENERAL PURPOSE DIGITAL COMPUTER PROGRAM FOR
SOLVING ORDINARY DIFFERENTIAL EQUATIONS
BY THE MATRIX EXPONENTIAL. METHOD

5. J. Ball R. K. Adams

ABSTRACT

MATEXP, a general purpose digital computer program, was
written for solving systems of ordinary differential equations
by the matrix exponential method. MATEXP has several advantages
over standard numerical integration routines. It gives virtually
exact solutions to constant-coefficient homogeneous equations
and to nonhomogeneous equations for which the forcing functions
are constant during the computation interval. The speed at which
the equations are solved and the accuracy of the solution are
essentially unaffected either by the degree of cross-coupling
of the equations or by whether or not the coefficient matrix is
nonsingular or that its eigenvalues are distinct.

The method has been extended to nonlinear equations and
equations with time-varying coefficients; this use is very
effective for engineering systems analysis problems.

NOTICE This document contains information of o preliminary nature
ond was prepared primarily for internal use at the Oak Ridge National
Laberatery. It is subject to revision or correction ond therefore does
not represent o final report.
 

 

 

LEGAL NOTICE

This report was prepored os on account of Government sponsered work. Meither the United States,

nor the Commission, nor any person acting on beholf of the Commission:

A. Mokes any warrenty or representation, expressed or implied, with respect te the accuracy,
completeness, or usefulness of rhe infermetion contained in this report, or that the use of
any information, cpparctus, mothod, or process disclosed in this report may neot infringe
privately owned rights; or

B. Assumes any liabilities with respact to the use of, or for domoges resulting from the use of
ony information, apporatus, method, or process disclosed in this report.

As used in the obove, "'person acting en behalf of the Commission' includes any employee or

contractor of the Commission, or employee of such contrector, to the extent that such employes

or confractor of the Commission, or employee of such contracter prepares, disseminates, or
provides access to, any information pursuant to his employment or controct with the Commission,

or his employment with such contractor.

 

 

 

 
‘CONTENTS

.Introduction.....l............-..C_..‘..t..l.llla..

Development of the Matrix Exponential Method

2.1 For Homogeneoué EquationsSeessseeescescceces

. 2.2 For Nonhomogeneous EqQuUatioOnSececssecsecacscas

2.3 Miscellaneous Features of the Matrix

Exponential..l.....l.l..l.........l.;...l

Description of MATEXP Program and Options

3.1 Basic Input Information........;...........

3,2 Alternative Methods of Generating the
Coefficient MatriX Accceececsccsescrccnns

3.3 Alternative Methods of Generating the
' Forcing Function Vector Zaeeeeeseeseesees
3.4 Methods for Solving Time-Varying-Parameter
| and Nonlinear Differential Equations.....

3.5 Spécial Forcing Function Subroutineés.esecese.
Sumnary and CONCIUSIONS e sssrerssnnsansesnsscanns

Appendix

5.1 Problems in_thé Evaluation of Exponéntial
Functions........................a;..;...

5.2 Detailed Description of ProgramsS.sececccesss

5.3 Fortran Listing of Programs..eccececscscsscss

Page

\O O O

11
13
13

15

16

19
20

27
28

28

30
13
1. INTRODUCTION

The matrix exponential method of solving differential equations
was first described to the authors by Prof: Henry Paynter of MIT,
who with his studentsl-3 developed this method into a practical
engineering tool. The basic technique was derived many years ago,
and even then it was an elegant method of obtaining exact solutions
for a set of constant coefficient, homogeneous differential equations.
The matrix exponential technique is ideally suited to digital
computation and is very simple to implemént, especially when compared
wifh most quadrature methods.

Ohly two persons besides Prof. denter have done extensive work

>

in this area. L. Pease” of Atomic Energy of Canada, Ltd., in-

dependently developed the method simultaneously with Paynter. The
work of Paynter and Pease formed the basis for our implementation

and, perhaps, refinement of the method, although the work of several

5-9

researchers established the rigor of the central technique.

J Suez, Automated Programming for Analog_Computers, M.S.
thesis, MIT, Aug. 1962.
2H.C.H. Lee, Some Finite Difference Models for Linear and
Nonlinear Control Studies Using Digital Computation, M.S. thesis,
MIT, Aug. 1962.

3H. M. Paynter and J. Suez, "Automatic Digital Setup and Scaling
of Analog Computers," Trans. ISA, 3, 55-64 (Jan. 1964).
hE. Artin, from O. Schreier and E. Sperner, Introduction to
Modern Algebra and Matrix Theory (1935); Translated from German,
Ch&lsea Publu CO., N-Yo’ 1951’ pp' 319‘320-

 

 

5L; Pease, DEEMS, A Fortran Program for Solving the First-=Degree
Coupled Differential Equations by Expansion in Matrix Series,
AECL-1898 (Oct. 1963, reprinted Feb. 1964).

6E. G. Keller, Mathematics of Modern Engineering, vol.Ill,
Mathematical Engineering, Wiley, N.Y., 1942, pp. 234-2L6.

TR Bellman, Introduction to Matrlx Analysis, McGraw-Hill, N.Y.
1960, pp. 165-173.

Ta
More recently, M. L. Liou of Bell Telephone Laboratories made 1mportant
contributions to the matrix exponentlal method. 10,11

Because this method can give virtually exact12 solutions to systems
of equations, it is of'considerable interest to most engineers engaged
in systems analysis, automatic control, and simulation. Also, systefis
engineers have long recognized that one essential difference between |
the analog computer and the digital computer is the:awkward (at best)
manner in which the digital macnine can performwintegration.‘ The
matrix exponential method, on}the‘other hand, requires the digital
computer to perform mainly matrix manipulations, which it can do in
a very straightforward and efficient manner.

The matrix exponential techniques have worked well for a large
general class of simulation problems which constitute the bulk of the
work in the systeme analysis>and automatic control fields. Indeed,
by use of the methods descrlbed in Sect. 3. k certain types of non-

linear equatlons can be solved as a natural exten51on of the basic

matrix exponential method.

 

8F, R. Gantmakher, Applications of the Theory of Matrices,
Interscience, N.Y., 1959, pp. 135-9 (translation of Russian
" original book: Theory of Matrices, 195k4). :

9L. A. Pipes, Applied Mathematics for Engineérs and Physicists,
24 ed., McGraw-Hill, N.Y., 1958, pp. 101-4,

Ly, 1. Liou, "A Novel Method of Evaluating Tran51ent Responses,

Proc. IEEE, 5u (1) 20-23 (Jan. 1966).

llF F. Kuo and J. F. Kaiser, eds., System Analy51s by Dlgltal
Computer, Wiley, .Y.,l966 pp. 99-129.

 

lE"Virtually exact" means that the solution can be calculated
to as great a precision as is desired, consistent with the precision
obtainable with & given computer word length. In other words, the
precision of the method is not necessarily limited by the convergence
of any approximate quadrature (integration) formula, simply because’
gquadrature is. not performed.
The matrix exponential meflhod has also been implemented and used
extensively in Fourier analyeis problems by simulating band-pass
fil‘c.e:rs.13’11L Instead of calculating correlation functions (and
subsequehtly their Fourier transforme) digital filtering can.be used
to obtain spectral density estimates and transfer functions from
noise data. Calculations using filtering technigues are of comparable
accuracy and typically more efficient than the conventional methods.
| MATEXP has also been used in a special technique to calculate the
sensitivities of the time response of a system to changes in parameter
values. 15 A descrlptlon of a subroutine which was written to
implement time response sensitivity caleulatlons is given in Sect.
5.2.3.

MATEXP has been developed and modified over a pefiod of several
~ years, and its present form reflects the considerable number of
helpful suggestions we have had from many people. We are particularly
grateful to Prof. H. M. Paynter.for first introducing us to the
method, and to Prof. T. W, Kerlin of the University of Tennessee,
and J. V. Wilson of ORNL for their help and encouragement.

2. DEVELOPMENT OF THE MATRIX EXPONENTIAL METHOD

2.1 For Homogeneous Equations

Consider the first-order scalar, linear, homogeneous differential

equation (with constant coefficient)

dx ' :
It + ax =0, - ) (1)

138. J. Ball, A Digital Filtering Technique for Efficient Fourier
Transform Calculations, ORNL-TM-1778 (July 1967).
lLLT. W. Kerlin and S, J. Ball, Experimental Dynamic Analysis of
the Molten-Salt Reactor Experiment, ORNL-TM-1647 (Oct. 1966).

 

 

Lo, oy, Kerlin, "Sensitivities by the State Variable Method,"
Simulation, 8(6), 337-345 (June 1967).

ra
whose solution is ‘ . o
x = e % (2)
o* '
An interesting characteristic of the solution is that, for any
time interval T, the value of x at the end of the interval is a
product of an exponential term e-aT_and the value of x at the beginning

of the interval, i.e.

X = € X, o (3)

This will be referred to as the "incremental solution.”
Now because & system of homogeneous linear equations of any
order can always be broken up into a set of first-order equations,

consider the following set of equations

dx : :

Lo x ta. x + a, X

dt 11 71 12 72 **** "1In "n’

dx2

TS %1 %1 *ag, Xy Foeees 8y X (4)
dx .

at "%t Tl Xp T overs By X

This array can be expressed compactly in matrix form as a first-
order, linear, homogeneous, matrix differential equation with constant
coefficients, i.e. | |

dX ' :
dt = AX 2 (5)

where X is the column vector of state variables Xi

50

33N'.... nfi
and A represents the coefficient matrix

all a]2 .o 08PN aln

a, a
A E 21 22 o e o0 98 a2n

anl an2 * &5 9 80 ann

This matrix equation has the solution
X, = X (6)

For a formal proof that Eq. (6) is the desired solution, the reader

is referred to Bellman.7 However, the following sinle proof is

somewhat less formal. First, if dX/dt = AX, then g{% = A.%% =
3 mx dt
AAX=4A" X; similarly, g;§-= a3 X, so that g—a-: A" x . (7)
dt dt
If Xt is expanded about zero in a Taylor's series,
x -x +F & | Lt ax . %
t 0 1! dt 21 th ' m! gt
t=0 t=0 t=0
With Eq. (7) substituted for the derivative,
2,2
At At
= +_ _+ * 8 9 80
S R Ty %o
or
At .
X, = ¢ X, (Q.E.D.) (8)
The "incremental solution" is
At |
Ko =€ %p (9)

A
where ¢ T, the matrix exponential, is defined analogously to the

scalar exponential as

32 3
eAT =1+ At + Légl—-+ ngl- + ..

X
4

in which I is the identity matrix

lOO .....O
Olo .88 89 O
0010 ... 0

*

O ..ll....ol

2.2 For Nonhomogeneous Equations

The matrix equation representing a system of first-order, constant
coefficient differemtial equations with nonzero forcing functions is

the nonhomogeneous equation

%=AX+Z, | (11)

where Z is the disturbanée, or forcing function,vector.

A general incremental solution of the ndnhomogeneous equation

as derived by Lioull is

C t+T R

AT A(t+'r)f =AT
= + ¢ . :
Xt+T € Xt € te _ ZT dt | (12)

An exact solution derived from Eq. (12) for the case where the

forcing function Z is constant over the interval t to t+7 is

At - AT -
Xt+T = ¢ Xt + (e -I)A th . (13)

It is important to note that the inverse of A need not be calculated

to evaluate Eq. (13) since

 

k! ?
10

2
= 7 :|:+-§':'E +£%T—L+

EAT}k-l
Kkl ?

o9 | k-1
Yy A (14)
k=1

Because this series is similar to that used to represent eAT,

the computer program can calculate the two required matrices

(c—:'}[\T-I)A-l series equals

concurrently, since the kth term of the
the (k-1)th term of the "' series times (7/k). In the MATEXP
program, the AT matrix is called the "C" matrix and the (EAT -I)A-l
matrix is called the "HP" matrix (in honor of H. Paynter).
At this point, two essential features of the matrix exponential
method are emphasized:
1. The exponential matrices can be computed by the series
approximation to nearly any,desired precision (typically,
1 part in 10° 1s specified for MATEXP calculations). Hence,
for homogeneous equations and for nonhomogeneous equations
in which the forcing functions remain constant over the
computation time interval, the solutions are virtually exact
solutions. _
2, The-solution vector can be updated successively»by a time

increment T by two matrix multiplications:

= +
XT- C XO HP ZO

XQT = C XT + HP 7.
eéc
If it i1s assumed that just one time increment value T‘iS
required, the C and HP matrices need to be evaluated only once.
An exact solution to the set of nonhomogeneous differential equations
can also be derived from Eq, (]12) for the case where the forcing |
function Z vqries linearly within the computation interval <.

In terms of the matrix exponential series approximations, the
.'.’

11

trapezoid forcing function incremental solution is

o0 _
A

AT 1.1 k-1
S zm)—) (hr)™" 2,
. k=1 .

20 () yk-1 |
T (et D)0 Ze4t ° | (15)4

k=1

C 11 : |
Liou = has also developed a recursive formula for accurate
approximations of continuous forcing functions which uses a Simpson's

rule approximation of the nonhomogeneous solution, Eq. (12), within

.the time interval =t:

At T 21 At/2 T
Xppr ™ € [Xt+6Zt]+3 € Zepcfo ¥ T Py (16)

As with the case of the step-wise varying forcing funétions, the
matrices required for Egs. (15) and (16) need to be evaluated just
once at the start. These features are not presently included in the

MATEXP code, but could readily be added as options. .

2.3 Miscellaneous Features of the Matrix Exponential

Since the matrix exponential prificiple has been a part of the
mathematical literature for many years, the matrix exponential has
had at least two other names: +the fundamental matrix, and the
transition matrix. Besides the series appfoximation method, an

9

analytical method is often used to calculate this matrix; however,
the eigenvalues of A and their eigenvectors must be calculated and
the initial condition vector must.be transformed by a matrix
comprised of the eigenvectors. It is emphasizéd-that the series
method used in MATEXP does not require that the coefficient matrix
be nonsingular (i.e., have a nonzero determinant) or that its
eigenvalues be distinct (a case where the analytical solution has
terms of the form te”" and cannot be expressed as the sum of
exponentials). The latter condition, which occurs in problems

where two time constants in a decay chain are egual, was one.of
12

the problems that Pease encountered in reactor burnup calculations -
that prompted him to develop the matrix exponential method.5

Another feature noted by Pease (but not included in MATEXP) is
that the average solution vector X could be obtained directly from
a matrix exponential type calculation,

From the mean value theorem,
T

X =%j' X, dt,
0

X can be obtained by integrating the equation for X in terms of C
and HP:

 

T T
= 1 1 [ 1
= = = = + .
X fot dt Tf cxo (HP) Zoj dt (17)
0 0
Term by term integration of the series approximations for C and
HP gives | '
. _ .
2 3
det='c I+§if—+'(§f) + A",f + ...|l=28P, (18)
o
and
; 2 |
fHPdt='12 "QI'TJ’E%JF A’f) T (19)
0 ,

The latter series, like the HP matrix calculation, could easily
be made concurrent with the other matrix exponential calculations.

The accuracy of MATEXP solutions, both in absolute terms and
compared with other methods, is difficult to estimate quantitatively
fér the general case. Even for those cases that are solved "exactly, "
the successive multiplications of the solution vector by the matrix
exponential naturally tend to accumulate errors. However, with
precise calculations of the C and HP matrices as recommended in the .
Appendix, Sect. 2.1, test cases have shown this error to be negligible
for large systems (L0 x LO), even after many thousands of updating
calculations. Lioull has developed an alternative method of évaluating
the C and HP matrices to a prescribed accuracy.

The nature of the matrix exponential method permits the use of
%

13

much larger computation time intervals 1 than would be feasible for
most numerical integration solutions. For constant-coefficient
equations and a given 1, it would be safe to assume that MATEXP would
be inherently mére accurate. As is usually the case, however, it
would be unwise to generélize about nonlinear equations. Nonlinear
solutions are discussed further in Sect. 3.4.

Eq. (20) gives a rough estimate of MATEXP so;ution times on the
IBM-7090 computer, assuming that a-negligible time is spent in the

peripheral subroutines:

Solution time(min) ® 3.0 x 1'66'(NE)2 NT , (20)
where NE is the number of equations, and NI is the number of
computation time intervals. For éxafiple, a 59 x 59 system run for
1000 time steps took 10 min, and an 8 x 8 run for 10,000 steps took
1.5 min., The solution time factor will vary from about 2 x 10-6 to
T x 10-6, depending on the amount of extra subroutine computation and
printout, and will be approximately halved for homogeneous equations.

The present "standard” version of the MATEXP program solves up
to 60th-order equations and uses about 22,000 words of‘core storage.
In a 32,000 word computer, the extra i0,000.words cafi be used for
special programming or storage, or the order of the equation’ can be
increased to about 80. Since, for larger probiéms, tape or other
slower storage devices would be required to calculate the matrix
exponential functions, the overall efficiency of the method would be
reduced. | |

Two other interesting, though perhaps purely academic, features
of the matrix exponential technique are that the solution timé
increment can be negative (allowing one to go backwards) and that the

A matrix can contain complex coefficients.
3. DESCRIPTION OF MATEXP FROGRAM AND OPTIONS

3.1 Basic Inpup Information

The MATEXP program was written with the intent that it should

be easy to use for a wide variety of differential equation problems.
1L

Unfortunately, as a program becomes more general, i.e. the more
options and special features the program has, it becomes more difficult
to explain the program and to use it for any given problem.
Consequently, any apparent awkwardness and complications in the
following discussion are due to a desire to make it general, and any
omissions are due to a desire to keep it simple.

The basic parts of the code are: +the main program, MATEXP; the
utility subroutine used for outputting, OUTPUT; and the subroutine
for calculating forcing(or disturbance)functions, DISTRB. To solve
linear, constant-coefficient differential equations that are
homogeneous (i.e. have no forcing functions) or which have only fixed
forcing functions, all the required data can be read'in>and no extra
programming is necessary. For equations of the form

%% = AX + 7,
the initial values of the X vector, the coefficient matrix A, and
the (fixed) disturbance vector 7 may be read in., Other information
required for each run is the following: |
1. number of equations,
2. initial time (or other independent variable),
3. compfitation time interval,
4., final time,
5. interval at which solution Vector X and disturbance vector z are
to be printed.

Since many elements of the coefficient matrix A are often zero,
only the nonzero elements need to be read in.. This makes it necessary
to identify each coefficient with its fow and column number. The
nonzero values of the initial condition and fixed disturbance vectors,
with theif row numbers, are read in similarly. |

Since successive runs might require no changes (or only a few)
in input data from the previous run, options are provided so that
only the altered data has to be read in. |

An option is also available whereby the last value of the X vector

from one run can be used as the starfing value of the succeeding run.
15

This option can be used if changes in the computation or printing.
interval are required in the middle of & solution or if certain
iteration or successive épproximation schemes are being used..

A complete description of the inputs and options 4is: given in

the Appendix, Sect. 5.

3.2 Alternative Methods of Generating the Coefficient Matrix A

Although the most straightforward method of inputting fifie
coeffiéien£ matrix 1s to read it in, very often it is advantageous
to have some or all of the elements calculated from system parameter
values. One option of MATEXP provides for this to be done by special
prbgramming on the first call of DISTRB.. An alternative is to use
an "algebra.table" routine developed by Kérlin and Lucius.16 This
routine calculates the matrix elements from‘input parameter values
without any specilal programming. The general expression uséd for
calculating an elemenf aij in terms of paramefiers Pk and their

exponents EkQ is

n1 Eo By . B3p Ep

E E E._ - E
11 o1 31
= 2 @& 8 + * 89 a & &
2, 5 ClPl P, P3 P CoPy P, 5 P +
or m n
a =) C, ;%fl | (21)
1] g .
1=1 k=1

A complete description of the program is given in reference 16.

Beside the fact that it is sometimes convenient to have the
coefficient matrix calculated by the computer, in some cases computer
computation is almost necessary to obtain accurate solutions. This
was the case fof one reactor dynamics calculation where the coefficients
were first carefully calculated on a 20-in. slide rule, then by the

machine. The difference in the steady-state solution for neutron

16T. W. Kerlin and J. L. Lucius, A Techniqpe'for Calculating
Frequency Response and its Sensitivity to Parameter Changes for Multi-
Variable Systems, ORNL-TM-1180 (June 1965).

 

 
16
level after a reactivity insertion was approximately a factor of 2.

3.3 Alternative Methods of Generating the Forcing Function Vector Z

When variable forcing functions are needed, a sfiecial program
must usually be written and included in DISTRB. Two‘special forcing
function subroutines have been written to simplify the programming:
DFG, for approximating arbitfafy functions; and TRIG, for approximating
variable transport lags. They are both described in Sect. 3.5.

For cases where the forcing function is a solution to an ordinary
diffefential equation, this equation can simply be added to the system
matrix, end an exact solution can be obtained. As an example, assume
that a sinuscidal forcing function is used to excite a damped spring-
mass system. The quadratic equation that describes the displacement

y of the mass with time is

2 - o
S¥ a4 vy =c sin (wt + ¢ | (22)
3t° dt C

where w is the frequency of the sinusoidal input (radians/time).

To arrange the equation in terms of first-order: derivatives, let

- 4y
xl' - dt J (23)
Xp 2 Y - (2k)
: 2 .2 B
Solving for d"y/dt" (or dxl/dt), we obtain
dx, |
T - bx, + ¢ sin (wt + ¢), (25)
and
dx
2
The equation for a pure oscillafor with frequency w is
2
d
Lo s=0- (e
17

- ds
If we let x3 = 3t and X, = ws, then
dx3'
x - —wx) (28)
dxLL
T - X3 e (29)

If the initial conditions of x_, and X) ‘are zero and -1, respectively,

3
then

sin wt , (30)

x3(t)

-cos wt . (31)

Xu(t)
. Thus éx3 could be substituted for ¢ sin (wt + @) in Eq.(25). The
required initial conditions of velocity xl(O) and displacement x2(0)
must also be specified.

The coefficient matrix for this example is

-a -b +c Ol

+1
A= 0 -
O 0 +w 0

If the sinusoidal input were introduced as a forcing functionm, it
would appear as a stair-step approximatiofi of a sine Wave, and the
accuracy of the solution would dépend on the accuracy of this
approximation. A comparison of the apfiroximate'and exact solutions
for a specific example is shown in Fig. 1. In the approximate
solution, a first-order extrapolation was uséd to approximate the
average value of the foréing function over the time interval.

In this example, the system has a natural frequéncy of 1.0
radian/sec and a damping factor of-0.25, and the driving sinusoid
has a frequency of 2.0 radians/sec. The computation interval of
0.5 sec for the afiproximatebcase gives about seveh cofipUtations
per cycle of the driving function. Figure 1 also shows the response

after a long time where the excellernt stability and accuracy of both
B ORNL DWG. 67-10215
1.0+ | ‘

 

 

  

 

S 1
Position - | N p] _ !
Xz 0 - ' ] T es [
| | | R Maximum error in
%\ ) Time (sec) approximate solution = 0.014
- Lk . ) .
™ Maximum error in initial transient ' Exact MATEXP solution
* .. approximate solution = 0.020 -
X Approximate solution, At = 0.5 sec
> : o ' :
. Fig. 1 - Comparison of Exact MATEXP and Approximate MATEXP
' . Solutions for Sinusoidal Input to Damped

- 1.0 3 o Spring-Mass System
19

solutions can be seen. This type of calculation is, historically,

very difficult to do wifih standard digital methods.l7.

3.4 Methods for Solving Time-Varying-Parameter and Nonlinear
Differential Equations

It was shown in Sect. 2 that the MATEXP method can provide exact
solutions to sets of constant-coefficient, homogeheous differential
equations and to nonhomogeneous equations for which the forcing
functions can be rebresentéd by stepwise-varying functions. oince
forcing functions are usually smoothly vérying,;the accuracy of the
solution would naturally depend on the accufacy of the stair-step
approximations.

Likewise, in the case of time-varying-parameter, or nonlinear,
equations, the variations in the coefficient matrix A can be
approximated by stepwise variations. For a variable A matrix, however,
the matrix exponentials (C and HP) would both have to be re-evaluated
at each computation interval. Although this may still be an efficient
method for low-order equations (~10 or less), it could be quite
time consuming for larger problems.

A more efficient method of solution is to modify, or "fudge;"
the forcing function vector so that it compensateé for the variation
in coefficients while the A, C, and HP matrices remain constant.

This is shown schematically in Fig. 2.

7R, A. Gaskill, "Fact and Fallacy in Digital Simulation,”
‘Simulation, 3 (5), 309 313 (Nov. 1965)
20

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Nonlinear Equations N o
Z(t )= A = £(t,X) —> X(t) (Exact)
Z(t) Ay - X(t) (Approximate)
+
6A X 96'
L

 

 

 

Fig. 2. Approximate Solution Using Fudged Forcing Functions.

Each component of the fudged forcing-function vector is calculated
by adding all the coefficient perturbation quantities in the row. TFor

example, assume one row of the matrix equation is

S dx
1
where a l’ al3, and z, are variables and 810 is a constant.
Let _ '
ay (8) = (a))g * agy
and
—_ !

Then the equation can be rewritten .

)

|
+ 2 (t) + all xl + al3 x3 .
_

= Zf(t:x)

. |
| — +
7o = (811)0 ¥ * 2p % ¥ (813)4 X5

 

Again, the foreing function zf‘would'éctfially be smoothly varying,

but in the MATEXP difference equations, 1t is approximated by a

stair-step functlon.
21

For the case where the coefficients and/or the forcifig functions
are known functions éf time, much greater accuraéy*(for a givén
computation interval T) reéults.ffqm.using'approximate.mean values,
rather than initial values, of the functions in the computation
interval. First-order approximations of the mean values can be
obtained by evaluating the time-varying forcing functions and matrix
elements at (t + 7/2) instead of at (fi). First-order extrapolations
of the mean values of the solution vector X should also be used

where coefficients are functions of X, as shown in Fig. 3.

Straight-Line
Approximation

e

 

 

 

 

 

 

 

 

 

P time

 

Fig. 3. First-Order Extrapolation of Mean Values of z and x at (t+%),

The use of an auxiliary subroutine VARCO greatly simplifies the
programming required to use first-order extrapolation calculations to
find approximate mean values of the forcing function. VARCO is
described in detail in Sect. 5.2.

The only way of guaranteeing that the solution is accurate is to
reduce the computation interval T'fintil further reductions make no

significant difference in the solution. A simple, intuitive estimation
22

of the accuracy, however, may be obtained by noting the maximum amount
of change in the solution and coeff1c1ent values within a computatlon
interval. If these changes are only a few percent of the values of
the functions at the start of the interval, then the flrst-order
approximations'will probably give very accurate answers. The'true
accuracy of the representatlon of a nonllnearlty should also be
considered when trylng to "squeeze" too much accuracy out of a
solution. - _ |

The use of fudged forcing functions for the soiutiqn of nonlinear
differential equations is'vefy effective when relatively few ef the
matrix coefficients are variable. In this case one fiight consider
the linear portion of the system of equations as being solved by an
extremely accurate analog computer, while the nonlinear portion is
simulated by a not-quite-so-accurate computer. If most of the |
matrix coefficients are variable, then the more conventlonal numerlcal
solution methods might be more practlcal than MATEXP

More detailed dlscus51ons of the theory and use of fudged forcing
functions have been found disguised in sophisticated mathematical

treatises by Wolf18 and Frazer et a1t

3.5 Special Forcing Function Subroutines

Since special programming is required in the DISTRB subroutine
to generate variable forcing functions for the differential equations,
two general purpose subroutines were written to facilitate this

pProgramming for some problems.

3.5.1 Arbitrary Function Generation - DFG

 

The arbitrary function generation-subroutine DFG provides a means
of generating approximations of singleevalued_functions of onef

variable where the arbitrary function curve is~represented by a

18 . | | |
, "A. A. Wolf, "Some Recent Advances in-the Analysis and Synthesis
of Nonlinear Systems", Am. Inst Elec. Engrs. transactlons paper
No. 61-713.

'9R. A;-Frazer, W. J. Duncan, and A. R. Cellar, Elementary
Matrices, Cambridge University Press, 1957, pp. 232-45.
23

series of linear segments (Fig. 4). The principle is identieal to
that of the diode function generator (hence DFG) used in.analog
computation.

Output
? Actual

Approximate

 

Q , S>> Input

 

Fig. 4. Subroutine DFG Representation of an Arbitrary
Function of One Variable.

DFG in its standard form arbitrarily allows for up to 8 functions
with up to 32 points (or 31 line segments) per function. Inputs
required are the ordinate and abscissa values of the line-segment
end points. If more functions 6r finer approximatibns are required,
the dimensions could be'changed eaéily. More details on the program

and a Fortran listing are given in the Appendix, Sect. 5.

 

3.5.2 Variable Transport Lag Generation - TRIG

A transport lag (also known as a pure time delay, or dead time)
acfuaiiy represents a distributed parameter system; hence, its
representation in a iumped-parameter solution will be only approximate.
The output =z from a pure delay device with an input x and a fixed

delay time T is
z(t) = x (t ~-1).

If v is variable, then the relationship between z and x is a function
of the time history of =.
The variablertime—delay problem is best illustrated by
fluid flow in a pipe where the inlet temfierature énd flow rate are
both variable. The assumptions required for a pure delay are:
1. there is no heat transfer to the pipe;
2. the fluid density is constant; -
3. plug flow exists, i.e., there is no mixing of the fluid in the

direction of flow.
2k

The technique used in TRLG is to sample the inlet temperature x
and the flow rate W at each computation time interval T, thereby
keeping an inventory on each slug of fluid in the pipe. The total
. weight of fluid in the pipe is éomputed from the initial transport
time T, and the flow rate W.: |

P, tal (1b) = W, (1b/sec) x T (sec) .
Similarly, the weight of fluid that enters during each time interval
T is W(t) x T. Since the fluid density is consfant, the weight of
fluid that leaves during that interval T 1s equal to the weight of
the inlet slug. | o o

As an example, assume that the temperature profile in the pipe
is as shown in Fig. 5 and the slug at the inlet of APO'lb is about
to enter. The slug at the Qutlet is APn at a temperature xn, wbere
APn >-APO. When APO enters, the outlet slug temperature will be
equal to X and the whole profilg will be shifted to the right
byAAPO
is then (APn - APO).

1f APO had been greater than APn, the outlet slug would have taken
as much of the upstream inventory (i.e., AP 1, &P .,

required (up to 300 samples), and the outlet slug temperature z

1b. The weight of the new slug just upstream of the exit

etc.) as

would be computed as the weighted average of the slug temperatures.
For example '
if

APO = APn + 0.5 APn_l ,

then -f--APn-xh + 0.5 APn-l Xn-l

Z = , .
, AP+ 0.5 AP |

If the maximum delay time (minimum flow rate)‘would use up too
many storage locations, the sampling would beldone every other (o:
every third, etc.) computation interval. With a variable lag, a .
minimum expected flow rate must be specified to calculate how often
to sample. | ,
The input variables suppiied by the calling program for each call
of TRIG are XT (e.g., fluid temperatures) and the flow rates W (in
Inlet

 

 

Outlet

 

 

 

 

 

 

 

 

Temperature
-
d};ff___r.._f._._.
. |

| i 0
v L : :
0 . . P i
Weight of fluid (1b) S “total

Fig. 5, Temperature Profile of .Fluid'in Pipe, -

gz
26

terms of mass/time, unity for full flow, or some percentage of full
scale). The lagged functions ZT are returned by TRIG.
On the first call of TRLG, the flag NI should be zero, and the

following input data are read in:

NLAGS = number of functions used,
TI = initial values of transport lag time for each function,
WMIN = minimum expected values of flow W. for each function.

The initial values of fluid temperatures in the pipes are set
equal to the initial values of inlet temperatures. If specific
initial temperature profiles are requifed,'they can be read in with
only a minor change being required in the program. The standard
version of TRLG provides for upAto sixvlags with up to 300 samples
per lag. If more or fewer lags or points are desired, the stétements
labeled DIMENS in the comment field can be changed accordingly.

"~ More details on TRLG and a Fortran listing are in the Appendix,
.oect. 5. .

There are two other techniques that are commonly used to represent
transport delays:

1. A series of n first-order lags, or "well-stirred tanks," with
time constants T/né o |

2. A Padé approximationfo which uses several terms of a series
approximation of e_TS_ (the Laplacian representatlon of a pure

delay) where S is the Laplacian argument .

Both the series lag and Padé methods have accuracy and flexibility
limitatiOns that would be prohibitive for certéin problems.21
Since the digital computer is quite proficient at sampling data,

- the sampled data approximation as used in the TRLG subroutine is

récommended as the most efficient and éccuraté method..

204, E. Rogers and T. W. Connolly, Analog Computation in

Engineering Design, McGraw-Hill, N.Y., 1960, pp. 419-2k.

 

 

213, q. Margolls and J. J. O'Donnell, "Rigorous Treatment of
Variable Time Delays", IEEE Trans. on Electronlc Computers, Vol.
EC-12, June 1963, pp 307-9.

 
27

L, SUMMARY AND CONCLUSIONS

The matrix exponential method has a numbef'of advantages over
the more common integration schemes for a large and significant class
of ordinary differential eQuation problems. The'speed ahd'accuracy
of MATEXP have'the pbtential of reducing computing costs for large
problems and of making more ''real-time" computations feasible for
on-line digital computation, control, and optimization calculations.

The MATEXP program has been developed over a period of several
years, mainly through use in simulation problems. There are, however,
at leasfi three other areas in which the matrix exponential method
might be effective:

1. Automatic parameter estimation - where the parameters of the
model differential équations are adjusted to optimize the
agreement betfieen theoretical énd experimental response curves.
A computer program to implement this technique is currently
under development;

2. Solution of nonlinear -algebraic equations by the method of
steepest ascents; and

3. Boundary value problens.

Other refinements that have been used with the MATEXP code
include the addition of an automatic plotting subroutine and a more
efficient output routine which prints only specified variables.
Forcing-function subroutines to solve implicit equations and
generate functions of two variables are planned as additions to the

"standard'" package.
28
5. APPENDIX

5.1 Problems in the Evaluation of Exponential Functions

The Taylor series approximation for a scalar exponential function

° ) ¢ Y. v
e’ n & =l+y+21+31+....+n1. (5.1)

This approximation also holds true when the argument y is a matrix;

hence, matrix exponential functions are amenable to digital computer
calculation, since raising a matrix to a power is a straightforward

operation. | |

It is important to note that the HP matrix calculation
HP = [exp (Ac) -'I]A'l | " (5.2)

does not require inversion of the A matrix, and can be calculated
directly from the terms of the C matrix approximation as shown
in Sect. 2.2. ,

There are several numerical problems associated with the matrix
exponential calculations. The approximations will be valid only if
1. the series will converge, - |
2. the numerical computation does not lose significance due to

overflow, roundoff, or truncation errors.

Since the evaluation of exp (Art) requires calculating:powers of the
matrix At, there is a practical limitation on the maximum value of:
the largest element in the At matrix, and experience has shown that
it is most efficient to limit this value to'about'l.O. Should the

desired T make max A, .7
1,
for the exponential calculations. The original arguments are

> 1.0, then T is halved up to 10 times

 

restored by applying the following equations as many times as

required:
29

C(Tj exp (A1)

- I (5.3)
= exp (A5) exp (A7)
HP(T).E [éxp (A7) -I] A_lt
o [exp (40 1] a7H[1 + exp (a3) -_ | (5.4)

There are also provisions in the code to keep track of the roundoff
errors in the exponential calculations. The ma.ximum values of the

largest elements in the QPT matrices Lé%%— are monitored to make sure

that they are not larger than the specified precision "P" times
lO8 (for an eight-decimal computer). When the QPT terms are summed;
the accuracy of thersummation will be approximately P, since thé
summation is carried out until the largest element in QPFT < P, If a
maximum value of a QPT element does exceed P x 108, then 7 is halved,
the exponential is calculated, and the original_T'is restored as before.
Users are caufioned that roundoff'erroré~may become significant
if restoration of the original T requires very many applications of
thé argument doubling Eqs..5.3 andv5.4. lWe know of no general rules
for estimating this limitation; however,'éhecks made on sample‘problems
indicate a "safe" boundary probafily existé'at a precision P = 10 ~ and
T halved 10 times. With a larger P and mOré halvings,ione should at
least be cautious about the results. ) |
The fidelity of the results are also questionable whenever the
ratio of the largest (absdlufe) matrix'eleméht to the sfiallest
(nonzero) element is > 10". This might be a manifestation of a very
wide range of time constants in a dynamicé problem. With a range of
~ 108, clearly the faster time constants could.be considered
"instantaneous" with respect.to the slower oqé$, and the equations

could probably be rewritten to get around-this problem.
30.

5.2 Detailed Description of Programs

Hopefully the information givefi'in this section is sufficiént to
permit the reader to use and modify MATEXP. Since we have tried
going through this typically excruciating experience with programs
from others, we have tried making things as clear as possible. In
particular, we have used many comment cards in the program listings
as a running explanation of what we are doifig; Either author would
be glad to try to help out any potential MATEXP user, and would be

happy to receive any suggesfiions for improving the program.

5.2.1 MATEXP Main Program

 

The MATEXP program consists of the main’progiam and two sub-
routines. OUTPUT and DISTRB plus any other sfibroutines called by
DISTRB. Even if DISTRB is not used, a dummy must be included.

For each case run on MATEXP, the data will include (if appropriate):
1. MATEXP Control Card,

2. Coefficient matrix (4),
3. Initial Condition Vector (XIC),
4. Any data read in by subroutine DISTRB,

5. Fixed forecing function vector (Z).

Input Data Formats - MATEXP Main Program

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

l. Control Card

Column | 1-2 6-7 11-20| 21-30 ) 31-40 | 41-50]| 51-60 | 61-62

Format | T2 13X {2 |3X| F10.0} F10.0] F10.0 | F10.0} F10.0 I2

Input | NE LL P TZERO T TMAX PLTINC | MATYES
Control Card - cont 'd

"Column | 63-64] 65-66 | 67-69 70 T1-72] 73-74 | 75-80

Format | I2 I2 I3 I1 I2 I . F6.0

Input | ICSS | JFLAG | ITMAX LASTCC TI17 | ICONTR VAR

 

 

 

 

 

 

 

 

 

 
31

NE

LL = coefficient matrix tag number

number of equations

P = precision of C and HP - recommend 10-6 or less
TZERO = zero time

T = computation time interval

TMAX = meximum time '

PLTINC
MATYES

1] = use previous A and T

I

printing time interval

coefficient matrix (A) comtrol flag

=-read new coefficients to alter A

= read entire rnew A (nonzero values)
DISTRB to calculate entire new A

= read some, DISTRB to calculate others
= DISTRB to alter some A elements
 ICSS = initial condition vector (XIC) flag

O 1 & W
il

1 = read in all new nonzero values
= read new values to alter previous vector
= use previous vector

= vector = O

v W
|

= use last value of X vector from previous run -
JFLAG = forcing function (Z) flag

1 thru 4 = same as for ICSS for constant Z

5 = call DISTRB at each time step for variable Z
ITMAX = maximum number of terms in series approximation of exp (AT)
LASTCC = nonzero for last case -
I1Z2 = row of Z if only one nonzérd, otherwise = O
ICONTR - for internal comtrol options -

0 = read new contrél card for next case

1l = go to 212 call DISTRB for new A or T

-1 = go to 215 call DISTRB for new initial conditions

VAR = maximum allowable value of largest coefficient matrix element * T
(Recommend VAR = 1.0)

ii

i
32

2. Coefficient Matrix A Formaet 4(213, E12,3) - Include if MATYES =

 

 

 

2, 3, or 5.
Column| 1-3 -6 ~ 7-18
Format I3 I3 B 12.3 ’ Repeat,
Input { Row No. Col. No.| COEFFICIENT 4 per card

 

 

 

 

 

 

 

Notes: 1. All row and column number entries on a card must
be nonzero. '

2. Insert blank card after all coefficient matrix
data is read in.

3. Data can be entered in floating point (F)
format with decimal point. :

3., Initial Condition Vector XIC Format (I2, 5(13; E12.3))- Include
if ICSS = 1 or 2 -

 

 

 

Column 1-2 3-5 6-17
Format 12 13 E 12.3 Repeat Cols. 3-17,
Input MM Row No.|} I.C. Value 5 per card

 

 

 

 

 

 

 

Notes: 1. Ail row number entries on a card must be nonzero.
2. Insert blank card after all XIC data is read in.

3. Data can be entered in F format.

L. Disturbance Vector Z Format (I2, 5(13, El2.3))- Include if
JFLAG = 1 or 2

 

 

 

Column 1-2 | 3-5 ' 6-17
Format I2 I3 : E12.3 "~ Repeat Cols. 3-17,
Input KK Row No. "Z Value 5 per card

 

 

 

 

 

 

 

Note: See notes under 3.

Two figures are included to aid in understanding the MATEXP
program. Figure 5.1 sumfiarizes the data arrangement, and Fig.
5.2 is a flow diagram.of the main program. The symbols used in
MATEXP are also listed and identified,
ORNL DWG. 67-10216

e 7 ac*

 

rMA’I‘EXP CONTROL CARD -

 

Case 2

 

 

 

Include if
JFLAG = 1 or 2

 

 

  
 
  

 

ICSS =1 or 2 7

 

 

MATYES = 2, 3, or 5

 

 

 

 

 

MONITOR
CONTROL CARDS

 

 

Fig. 5.1 MATEXP Data Arrangement

£e
34

ORNL DWG. 67-10217

 

 

 

NI:O FROM BOTTOM RIGHT
Fig. %.2¢

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

CALL DISTRB
IST CALL

 

 

 

JIFLADO 7o Yop
PTMP + P 109 Fla.2.2s
PRNT CONTROL OATA

PLTING ® PLTINC #9999

FK=0

 

 

 

Fig. 5.2a. MATEXP Block Diagram — Read or Compute A Matrix and XIC Vector.
35

ORNL DWG, 67-10218

FROM BOTTOM
Fig. 8.20

  

 

 

FIND AMAX & AMIN
‘RATIO » AMAX / AMIN
T HALVED ISTOR TIMES
(Ur T8 10) UNTIL
AMAX % T { VAR

 

 

 

 

TO STATEMENT 20
Fie.5.2C

 

 

 

 

 

 

 

 

 

          

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

PE+ O
AL+ 1.0
aoPT s
47
2
00 16 KL » |, ITMAX
KLM ». KL C(ZT) =« CITIHC(T)
ALL - T/aL ToSTATEMENT 3T
AL ® AL+ 418
TALLL » T/AL : : .
QPT * QPTHA ¥ ALL [ He2T) = HPITI+CITINNR(T) |
€*C+0QPT
- a9
o+
/1
(uFK-7)
HP ¢ HP + OPTHTALLL
47 .
' K = SR+
[ PMK-aBS (QPT(1MAX, JMAX)) | TeT#0.8
¥ o
(QPTMP-PMK) .83
+
802 ,
1STOR:
6 IF- 414
{ ( PMK-P) ‘ : ISTOR+JFK
0~ . 4,0
" :
1F e ' IF
{ PE—- 2% PMK ) 4 PRINT KLM ITMAX Db IKLM= -
0,+ ] -1) 533 ITMAX
0~ '

Fig. 5.2b. MATEXP Block Disgram — Compute C and HP Matrices.
36

FROM RIGHT SIDZ
FiG. 5.2b

20

 

TIMETZERO
PLT:0

 

 

 

6 | GO TO 54 89

JFLAG

25
Gl [ kw0

READ Z -0
271

CALL QUTPUT
1 ST CALL
(NI SET *t)

 

v

 

 

 

 

 

24

ORNL DWG. &7-1029

IF -
(IIA'IJES 9 25

+

CALL DISTRB
18T CALL

 

 

    
 
 
 
   

 

| Yo YHHPHZ(IIZ) |

 

[ YeCHX

 

 

 

SOLUTION
Xo¥

 

 

 

{

JIFLAG E)
TIME » TIME+T
PLT=PLT+T

 

 

 

 

 

 

 

 

CALL DISTRD

 

 

 

TO STATEMENT |
FiG. 8.2a

 

TO STATEMENT 212
FiG. %.20

 

 

 

PLT=0Q
KeK+)
N[O

 

 

 

 

38
FROM BOTTOM
Fi0. 5.2b
¥ - '
TIME-TMAX
0+
37

Fig. 5.2c. MATEXP Block Diagram — Compute

 

Solution Vector.
37

MATEXP MATN PROGRAM SYMBOL KEY

 

1.

2.

3.

Control Card Inputs

 

See input data format list.

Input Data
A(NE,NE) = coefficient matrix

MM = initial condition vector tag number

XIC (NE) = initial condition vector

KK = disturbance vector tag number

Z(NE) = disturbance vector

" Internal Variables

 

The following variables are listed in alphabetical order.

ADT = AMAX % T

AL = Floatlng point KIM for ALL ‘cale, KIM+1 for TALLL

ATL = T/AL with AL = KIM

AMAX = Maximum (absolute) value of element in A matrix

AMTN = Minimum (absolute) value of nonzero element in A matrix
C(NE,NE) = Coefficient matrix éxponéntial

HP( NE,NE) = Disturbance functiofi matrix exponenfial

IMAX = Row location of AMAX | |

IMIN Row location of AMIN

ISTOR = Number of times matrix exponential argument T is
halved so that AMAX x T¢VAR; later ISTOR = ISTOR + JFK

JFK = Number of times T is halved in order for matrix exponential
calculation precision to be P or better

'JJIFLAG = Flag to prevent double call of DISTRB durlng initial
time step calculation

Column location of AMAX
Column location of AMIN

il

JMAX
JMIN

K = Case number

i

KIM = Number of terms. in series approximations of exponentials

NI = Printing flag: O on initial call of OUTPUT causing printout
of A, C, and HP matrices., OUTPUT sets NI = 1 on first call.

PE = Maximum element in. (n - l)th QPT term
PMK = Maxlmum.element in nth QPT term
38

QPT(NE,NE) = Term in series approximation of C matrikx

QPTMP = Maximum perm1551ble value of element in QPT matrix.

RATTO = AMAX/AMIN. If RATIO less than 10° (for eight decimal
machine) there may be significant problems in
calculation of C and HP. o

TALLL = T/AL with AL = KIM +1

TQP(NE) = Temporary sfiorage for QPT terms

X(NE) = Solution vector B

Y(NE) = Temporary storage for X

 

5.2.2 Subroutine QUTPUT

The first time MATEXP calls OUTPUT, the coefficient matrix (A)
and the exponential matrices C and HP are printed out, alofig with the
_initial solution (X) and disturbance (Z) vectors. OUTPUT also sets
the first call flag (NI) to 1, and on subsequent calls only the X
and Z vectors are printed. A possible means of saving computing
time at the expense of storage would be to store X (and Z) values
-in-arrays for a large number of time intervals; then print the
-arrays out in blocks. Additional savings could be achleved by

printing only selected variables.

5.2.3 Subroutine DISTRB

 

Subroutine DISTRB may be called by MATEXP either to compute
matrix coefficients (A) on the first call (i.e. when flag NI = 0O)
and/or compute variable forcing-function vectors (Z). |

Other special purpose subroutines, such as VARCO, DFG, TRIG,
and any others the user may want to supply, are usually called by
DISTRB. '

Another special purpose use of DISTRB is to compute inputs
for successive MATEXP cases without requiring a control card for
each case, This is done by means of the flag ICONTR (Cols. 73-L4 on
the control card). After a case is run, the first call flag NI is
reset to O, and case number K is increased by 1l; then if ICONTR

is positive, DISTRB will be called at statement 212, where a new

 
39

coefficient matrix A or time interval T may be qalculated.' If
ICONTR is negative, DISTRB is called at statement'QlS, permitting
new initial conditions to be used.

. The program listing for DISTRB that was used in calculating the
sinusoidal forcing function for the example in Sect. 3.3 is given
in Sect. 5.3. |

Another veréion of DISTRB is used to calculate the sensitivity
of a system's time response to changes in the system's coefficient
matrix elements

ax .

da, .
1J

 

- DISTRB controls the solution of the sysfem~equations and stores
those values of the solution vector which are to be used subsequently
as forcing functions for the sensitivity calculations. To compute
the sensitivity to aij’ the jth row of the system solution£zector
is stored and is later used as a forcing function to the i row of -
the same system eq_ua'tions.l5

After solving the system equations and storing the required
elements of the response vector, the arithmetic average values of
the X's in each time interval are calculated and stored (XT).

.During each sensitivity run, DISTRB feeds the forcing function
into the system equations, and the resulting printouts of the X
vectors are the desired sensitivities.

For the sample program shown in the Fortran listing, Sect. 5.3,
the system is forced by a unit step input in row I1Z (specified on

the control card). Other control card inputs are:
JFLAG = 5

JCONIR = 1

Special input data read in by DISTRB are the row (IS) and column
(JS) numbers of the matrix elements for which sensitivities are to
be calculated, the number of time points (NTS), and the number of

sensitivity runs (NSENS), as follows:
 

1 | A 1 | 51
fs() | os(a) | (x) | 18(2) | gs(2) | (wX)f...thru Js(5) |NTT | NSENS]

 

T3 13 13 13 I3 T3

5.2.4 Subroutine VARCO

 

The VARCO (VARiable COefficient) subroutine can be used with
DISTRB to simplify the programming of problems with variable coefficient
matrix elements. In general, these elements are functions of both
time and the values of the éolution vector X. VARCO is designed to be
called by DISTRB at the start of each computation interval and to
return the mean values of time (TX), and X, (XTR), for that interval.
The mean values. of X are predicted by a first'brder-extrapolation
scheme, as shown in Fig. 3. VARCO will also cause the initial time
step to be repeated, using the first try at calculating X(T) to
estimate the mean value at g. DISTRB can then calculate the
coefficient values using TX and XTR. Use of this first-order _
extrapolation scheme results in significant improvement in accuracy

over using no extrapolation.

5.2.5 Subroutine DFG

DFG uses the principle of the analog computer's Diode Function
Generator (see Fig. 4) and uses linear interpolation to approximate
arbitrary, single-valued functions of a variable. Data for DFG is
read in the first time it is called by DISTRB (i.e., when NI = 0).
The standard program provides for up to 8 functions with up to 32
coordinates each. o N

On each successive call, DFG returns the functions ZD for
varying inputs XD, If an input XD(I) goes outside the specified
limits, the output is a straight-line apbrokimation of ZD(I) based
on the slope of the function at the boundary, and an error message
"DFG(I) RANGE EXCEEDED" is printed.

The inputs read in by DFG are:

NDFGS Number of functions used

NPTS(8) Number of points in approximation for each function
41

XP(32,8) Independent variable points ..
7ZP(32,8) Dependent variable points

The input format is as follows:

"-Card No. 1 (I2, 8X, 8I3)

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Column 1-2 11-13 '
= Repeat Cols, 11-13
Format I2 8X I3 ' . '
\ 7 more times for
Variable | NDFGS NPTS(1) NPTS(2) to (7)
Card No. 2, 3....etc. (8EL0.3) o
~ Column 1-10 11-20 21-30 31-40 Repeat as required
Format E10.3 E10.3 E10.3 E10.3 | for DFG(1); Max.
Variable | X2(1,1) | z8(L 1) | x2(2,1) | zp(2,1 | O mumbers per card
NOTES: 1. When all data for DFG(l) has been-ehtefed, start

DFG(2) data on new card; etc.

2. Enter independent variable points XP in order,
progressing from most negative to most positive

values.

3. F Format entries (with decimal point) may be used.

5.2.6 Subroutine TRLG

 

TRLG (TRansport LaG) is described in some detail in Sect. 3.5.

.~ The input functions XT (e.g. fluid temperature) and the mass flowrates

W (in terms of either mass/time, unity for full flow, or some

percentage of full scale) are supplied by the calling program DISTRB,

and the lagged functions ZT are returned by TRLG.

On the first call

of TRLG (when NI = O), the following input data is read in:
NLAGS Number of functions used

TI(6) Initial value of transport lag time for each function

WMIN(6) Minimum expected value of mass flow W for each function

The program is set up assuming that subroutine VARCO is also

called by DISTRB.

VARCO has a restart feature which repeats the

initial time step calculation; thus the TRLG functions will not be

updated on the second call. If VARCO is not used, this second call

 
L2

omission may be deleted by removing statement 33 in the TRLG program.
The input format for TRLIG is:

Card No. 1 (I2)

 

Column 1-2

 

Format 12

 

 

 

 

Variable NLAGS

 

Card No. 2 (6E10.3)

 

Column 1-10 Repeat 5 moref|

 

Format ¥10.3 times for

 

 

 

 

 

Variable | TI(1) TI(2) - (6)

 

Card No. 3 (6E10.3)

 

Column 1-10 Repeat 5 more

 

Format | £E10.3 times for

 

 

Variable WMIN(1) WMIN(2)7'~(6)

 

 

 

 

 
-43-

5¢3 -~ FORTRAN LISTING OF PROGRAMS

$IBFTC MAIN DECK

fi(‘l(‘)(‘b(‘)fififififififlflfiF\(\(‘s(\fi(‘lflfifififi(‘l(‘lflffififlf‘l(‘)fi(‘)dfifl(‘)(‘\fifififi(‘l(‘lflflflfi'fififlfi

PROGRAM MATEXP FOR THE 7090 - FORTRAN 4

THfS PROGRAM CALCULATES THE SOLUTION OF A MATRIX OF FIRST
ORDERs SIMULTANEOUS DIFFERENTIAL EQUATIONS W/ CONSTANT COEFFICIENTS

OF THE FORM DX/DT # AX + Ze«
THE METHOD IS PAYNTER—S MATRIX EXPONENTIAL METHOD

THE SOLUTION IS GIVEN FOR INCREMENTS OF THE INDEPENDENT
VARIABLE (T) FROM TZERO THROUGH TMAX

COMPUTES MATRICES C # EXP(A#T) AND
HP # (C-I1)*A INVERSE
SOLUTION X(N#T) # CH*X{(N—|)*¥T)+HP*Z ((N—1)*T)
SERIES CALCULATION OF C AND HP MONITORED TC
ASSURE SPECIFIED SIGNIFICANCE.
IF T IS REDUCED FOR C AND HP CALCSe>s
CRIGINAL ARGUEMENTS ARE RESTORED BY -
CL2¥TH#C(THy*C (T
HP (2#T)#HP(T)+C(T)#*HP(T)

OUTPUT FROM THE PROGRAM IS PRINTED AT INTERVALS PLTINC.
THE PROGRAM USES SUBROUTINES DISTRB AND OUTPUT

INPUT FOR THE PROGRAM CONSISTS OF
- " ONE CONTROL CARD :
THE COEFFICIENT MATRIX A (UP TO 60 X 60) o DIM
THE INITIAL CONDITION VECTOR X :
A FIXED DISTURBANCE VECTOR <Z

A VARYING Z CAN BE GENERATED BY DISTRB
VARIABLE COEFFICIENT EQUATIONS MAY BE SOLVED BY APPROPRIATE
FUDGING OF THE DISTURBANCE FUNCTION SUBROUTINE.

CONTROL CARD INPUT INFORMATION

NE#NO. OF EQUATIONS (12)
LL#COEFFe MATRIX TAG NCe (I12) -
P#PRECISION OF C AND HP (Fi0e0) -~ RECOMMEND |e«0E-6 OR LESS
TZERO#ZERO TIME (Fl0.0) ,
T#CCMPUTATION TIME INTERVAL (F|10s0)
TMAX#MAXIMUM TIME (FlQ0e0)
PLTINC#PRIMTING TIME INTERVAL (FiDe0)
MATYES#COEFFe. MATRIX (A) CONTROL FLAG (I12)
| #USE PREVIOUS A AND T
2#READ NEW COEFF.S TO ALTER A
3#READ ENTIRE NEW A (NON-ZERO VALUES)
4#DISTRB To CALC. EMTIRE NEW A
E#READ SOMEs DISTRB TO CALCe OTHERS
6#DISTRB TO ALTER SOME A ELEMENTS
ICSS#INITIAL CONDITION VECTOR (XIC) FLAG (I2)
| #/READ IM ALL NEW NON-ZERO VALUES
2#READ NEW VALUES TO ALTER PREVIOUS VECTOR
3#USE PREVIOUS VECTOR
LEVECTOR#D
YOO OO CTYO OO OO DYy Oy Y O

™y

Oy OOy

FH@

OO0

o

B#USE LAST VALUE OF X VECTOR FROM PREVIOUS RUN
JFLAGH#FORCING FUNCTION (2) FLAG (12}
| THRU 4#SAME AS FOR 1CSS FOR CONSTANT Z
S#CALL DISTREB AT EACH TIME STEP FOR VARIABLE Z
I TMAX # MAX. NOo OF TERMS IN SERIES APPROX.
OF EXP(AT). (I3) |
LASTCC # NON—-ZERQ FCR LAST CASE (1)
11Z # ROW NOe OF 2 IF ONLY ONE NON-ZEROs
OTHERWISE #0O- (12)
"ICONTR — FOR INTERNAL CONTROL CPTIONS (12)
O#¥READ NEW CONTROL CARD FOR NEXT CASE
| #GO TO 212 CALL DISTRB FCR NEW A OR T
~1#GO TO 215 CALL DISTRB FOR NEW I+Ce=S
VAR # MAX. ALLOWARLE VALUE OF LARGEST CQEFFe MATRIX ELEMENT * T
(RECOMMEND VAR#1.0) (F6.0)

DIMENSION A(éD;éD),C(GD’éD);HP(éflgéfll,OPT(60960)§
IX{60)sY(60)+2(60)sXICL60),TQP(60)

COMMON CsHPsA»sQPTsXsZsY s ITMAX sKKsLL s MM,
| JUFLAGsXICyNT s TIME s TMAX s TZERCSNEsTQP T
211ZsI1CONTRsPLTINCIMATYESSI1CSSsJFLAGHPLT

K#CASE NUMBER
NI#0 ON |-ST PASS. SET TO | ON 1-5T CALL OF OUTPUT.
K#|
NI #0
| READ (5,100) NEsLLsPsTZEROsTs TMAX sPLTINCsMATYES s 1CSSs

IJFLAGs ITMAX sLASTCCs 1125 I1CONTR VAR
100 FORMAT(2(12+3X) s5F 10e0s312s13s1142125F6eC)

COEFFICIENT MATRIX INPUT
GO TO (3+9992s2+2+3)sMATYES
DO 9n J#!1sNE
90 A(IsJ)#OWD
IF(MATYES-4)9953599
99 DO 9| I1#1+1379
MATRIX ELEMENTS 5(ROW,s COLUMNSs VALUE)
ALL I AND J ENTRIES ON CARD MUST BE NON-ZERO.
A BLANK CARD IS REQUIRED AFTER ALL ELEMENTS ARE READ IN.
READ (5s101) IleJl oDl 9129329029139 J3sD3914sJ4,D4
|0t FORMAT (4(213sE1243)) ' ,
IF(I1)2+3492
92 AlI1,J11#D!
ACT2,9J2)V%D2
A¢q3’JBfi#D3
91 A(14,Jޤ#o4

INITIAL CONDITION VECTOR XIC INPUT
3 GO TO(4s 12096954861 91CSS
4-DO 93 1#14NE
93 XIC(I}#0e0
120 DO 94 I1#1,15 | '
ALL ROW (I) ENTRIES MUST BE NON=-ZERO
A BLANK CARD 1S REQUIRED AFTER ALL ELEMENTS ARE READ IN.
READ (5%,+95) MM 111 sD1[s1124D125113,D1351143D14s1155D15
95 FORMAT(I245(134E12e32}) '

 

DIMEN
DIMENS
Yy () (Y O

N

96

94

81
82
214
212
213

5.

IF (1111636396
XIC(II1)Y#DI I
XIC{II2)#D|2
XIC{II3}#D13
XIC(II4)1#D 14
XICIII5)#D15

MM#0
DO 7 I#1sNE
XIC(I)#0.0

IF(ICSS~5)814214481I
DO 82 I#14NE
X(CT)#XIC(I)
IF(MATYES=3)21352135212
CALL DISTRB
JJFLAG#D | o _ -
QPTMP # MAXe. PERMISSIBLE ELEMENT OF QPT FOR 8 DECIMAL COMPUTER
MATRIX CALCe LOSES SIGNIFICANCE IF LARGEST '
ELEMENT IN SERIES APPROX. MATRIX QPT 1S
GREATER THAN P#*|.CES
QPTMP#P* | +E8 -

WRITE (6+211) KeNEsPoTo
IPLTINCyMATYES s ICSSsJFLAGYICONTR s ITMAX, I|Z’VAR9QPTMP

211 gFORMAT ( I2ZHIMATEXP CASELI3/17H NOe. OF EQUATIONS,

113/720H SPECIFIED PRECISIONsFI248/6H TIME o

28HINTERVAL,F 1848/ 15H PLOT INCREMENTF|7e8//

316H CONTROL FLAGS =/1H 35Xs6HMATYESs14/1H

45X s 4HICSSs 16/ 1H $5Xs5HUFLAGsI5/1H s5Xs6HICONTR 14/
534HOMAXe TERMS IN EXPONENTIAL APPROXaesI5/

806

402

40

407
408

613H SINGLE Z ROWsI14/20H MAXe ALLCWABLE A#DT9sF9.3/
727H MAXe ALLOWABLE QPT ELEMENTsFlle3)

PLTINCHPLTINC*#0+9999

JFK#D

IF(MATYES—=1)20+20+806

SCAN MATRIX FOR MAXe AND MINe NON-ZERO ELEMENTS.
IMAX# |

JMAX# | B

AMAX#ABRS (Al 1ls1))

DO 401 I#14NE

DO 40! J#I»NE

IF(AMAX=ABS (A(I1sJ7))402s401 401
AMAX#ABS (A(IsJ))

IMAX#I :

JMAX#J

CONTINUE

IMIN#IMAX

JMIN#IMAX

AMIN#AMAX

DO 4N9 1#I| sNE

DO 409 J#1»NE

IF(A(IsJ)) 40794099407 —
IF(ABS (A(I1sJ))~AMIN) 408,409,409
AMIN#ABS (A(IsJ))

IMIN#I

JMIN#J
Oy YO

Yy (v

‘)

()

(1

409

413
403
404

405

410

46~

CONTINUE
RATIO#AMAX/AMIN
AMIN # MINIMUM NON-ZERO ELEMENT
ISTOR#D
ADT#AMAX*®*T
DO 403 I#!1,1!
IF{VAR=ADT) 413+404s404
ISTOR#ISTOR+I
ADTHADT¥D 5
T#ADT /AMAX
COMPUTATION INTERVAL T IS HALVED ISTOR
TIMES (1g#MAXe) SO MAXe ELEMENT IN A®*T
IS LESS THAN VARe ,
WRITE (6s405) IMAX s IMAXsA{IMAX s IMAX) 9sADT o T
| IMINesJMINSZA({IMINsJMIN)LZRATIO
FORMAT (3 |HOMAX<COEFFe MATRIX ELEMENT # A(sI123s1HssI1293H) #5»
| El5e4/13H MAXe A#DT # oF 128 42Xs l4HWITH DELTA T #,Fi5.8/

230HOMINIMUM NON-ZERO ELEMENT # A(sI12sIHss129s3H) #sE15e4/
318H RATIO AMAX/AMIN #sEI5.4)

IF(ISTOR-10)8s410s410
WRITE {(6s411)

411 gFORMAT (34HDA*DT STILL GREATER THAN ALLOWABLE,

48

49

| 1.9H AFTER 10 HALVINGS.)

GO TO 37 _

CALCULATION OF MATRIX EXPONENTIALS C AND HP
DO 9 I#IsNE

DO 9 J#I1 sNE

ClIsJ)#Do

DO D I#1sNE
C(IsI)#1a

SKIP HP CALCS. FOR HOMOGENEOUS EQUATIONS
IF (JFLAG=4)48,51,48

DO 49 1#1sNE

DO 49 J#1sNE

HP (I +J)#0e

DO 50 I#! 4NE’

50

51

HP (14 1)#T

PE#0.0

DO 11 I#!,NE.

DO |1 J#1sNE

QPT(IesJ)HCI(I o J)}

FORM THE MATRIX EXPONENTIALS CHEXP(A*#T) AND HP#((C—I1)%A INVERSE)
ALZ! oD

DO 16 KL#IsITMAX

KLM#KL

CALL#T /AL

AL#AL+! 0
TALLL#T/AL

 

e e ——

~e?

Lo
— Ll e
'.")-I-T"

DO 18 I#!4NE

C
C
DC |3 J#1sNE
TQP(J)#DeD
DC 13 KX#1 4NE
13 TQP(J)YETQAP(JI+QPT(IsKX)*¥A(KXsJ)
c o AR
DO 18 J#I sNE
18 QPT(IsJI#TQP(J)*ALL
C .
C QPT#MATRIX TERM IN SERIES APPROXe #( (A*T)%*%#K)/K FACTORIAL
DO 44 I#1sNE
DO 44 J#| oNE
44 C(IsJIHC(II 4y I+QPT(IsJ)
C
[F (JUFLAG-4)45 447445
C

45 TF(ITMAX=KL)47 9474145
145 DO 46 I#I| sNE
DO 46 J#I1sNE
46 HP(I s JYHHP (I s J)+QPT (s J)*TALLL

FIND MAX ABS ELEMENT IN QPT AND CALL IT PMK

Y YOy M

LARGEST QPT ELEMENT USUALLY IN ROW IMAXs COLUMN "UMAX
47 PMK#ABS (QPT(IMAXsJMAX))
IF(QPTMP~-PMK) 83,83s502
502 IF(PMK—=P) 40640616 ' ,
C SCAN OTHER QPT ELEMENTS ONLY WHEN QPT(IMAX, JMAX) IS LESS THAN P
406 DO |4 I#I1sNE -
DO 14 J#I oNE
|4 PMK#AMAX | (PMKsABS (QPT(IsJ)))
IF(PMK=P)1T7s17s16

C
C PRESENT MAXe QPT ELEMENT SHOULD BE LESS THAN
C HALF PREVIOUS MAXe TO INSURE CONVERGENCE
|7 IF(PE=2e%¥PMK) 16921321 ~
_ |6 PE#PMK
C
21 WRITE (6+200) ~ KLM
c |
200 FORMAT (44HQONCe OF TERMS IN SERIES APPROXe OF MATEXP # 412)
C : _ e
IF(ITMAX=1)20920+538
538 IF(KLM=1TMAX) 4144+83,83
C
JFKHEJFK+ 1
IF(JFK~T7)303+304+304
304 WRITE (6+305) PMK
305 QFORMAT(32HD7 TRIES AT HALVING T NeGes PMK#sF1246)
GC TO 37 |
303 WRITE (6+210) KLMsPMK o T

210 FORMAT(2!HOMAXe ELEMENT IN TERM,I13,8HOF QPT #sE1l+3/
| 35H TRY HALVED TIME INTERVAL DELTA T #,Fi15.8]
GO TO 8 |
-48-

414 1STOR#ISTOR+JIFK

C ORIGINAL ARGUMENTS OF C AND HP MATRICES RESTORED IF ISTOR GREATER THAN O ,
IF(ISTOR) 20205416 '
416 WRITE (64+415) ISTOR

415 FORMAT(26HOTOTAL NCe OF T HALVINGS #,413)
DO 417 KR#|sISTOR -
IF(JFLAG=4) 41944189419 :
C SKIP HP CALCS. FOR HOMOGENEOUS EQUATIONS
419 DO 420 I#I1sNE '
DO 421 J#I| 4NE
TQP(J)Y#0e0
DO 421 KX#I1 sNE
421 TQP(JIATAP(J)+HP (1 sKX)*C (KX J)
DO 420 J#1 4NE
420 HP (I s N #TQP(J)+HP( I+ J)

418 DO 430 I#1,4NE
DO 430 J#!14NE
430 QPT(IsJ)#00
DC 431 I#I1 4NE
DO 431 J#! 4NE
DO 421 KX#| sNE
431 QPTU(I s JI#QPT( L 9J)+CITIoKX)¥C(KXsJ)
DO 432 1#I| 4NE
DO 432 J#| 4NE
432 C(LsJIHQPT(I4J)
417 TH2.0%T

ClIlsJ) IS THE MATRIX EXPONENTIAL CHEXP(A%*T)
AND HP(IsJ) IS THE ((C~-I1)Y*#A INVERSE) MATRIX
NOwW WE READ (OR CALL SUBRCUTINE FOR) DISTURBANCE VECTOR

P OYOOY N

20 TIME#TZERO
PLT#O.
GO TC (264121 927925955),,JFLAG
55 IF(MATYES-3)2155215427
215 CALL DISTRB
11Z#112Z
GO To 27

26 DO 97 1#I1,NE
97 Z(1)#0.0
121 DO 98 I#1415 »
C ALL ROW (1) ENTRIES MUST BE NCN-ZERO
C A BLANK CARD 1S REQUIRED AFTER ALL ELEMENTS ARE READ INe
READ (54595} TKKeI1219D211224sD224123sD2351249D24451254D25
IF(I21§27+27+78 - '
78 Z(121)#D21
Z(122)V#D22
Z(123)#D23
Z124)#D24
98 Z2(125)#C25

25 KK#D
DO 28 I#!sNE
28 Z(1)#Q. |

C ON |-ST7 CALL OF OQUTPUT NI SET TO |
27 CALL OUTPUT -

 
O OO

OO0 M

53

56
30
29

702
703

700
32

52
31

ONE
NOW

33
35

37
34

40

-49-

COMES THE EQUATICN SOLUTION BASED ON
XINT)EMEX(NT=1)+((M=1)A INVe)¥Z(NT~1)

IF (JFLAG=4)29454456

DO 53 I#Is4NE -

YUIVHC T s 1y #X U]

DO 53 J#2sNE

YOI)AY I )Y +C (I s J) %X (J)
IF{112)52452,702
IF(JJFLAG)I30s2923D

CALL DISTRB

IF{I11Z2)700+700s54

ONLY ONE Z-TERM CALC, IF I1Z 1S GREATER THAN ZERO
DO 703 I1#I1sNE

YOI #Y (DI +HP (T T1Z)I%Z (1 1Z)

GO To 52

DO 32 I1#!IsNE
YUIYHCITI o 1) #X ) +HP LT 1) *Z (1)

DO 32 J#2sNE :
Y(I)AY(I)+C(Tad) X (J)+HP (T s ) ¥Z (J)
DO 31 I1#I sNE

X{Iy#Y (1)

TIME INCREMENT OF THE SOLUTION HAS JUST BEEN FOUND
PLOT AND PRINT IF PLTINC INTERVAL HAS ELAPSED

JJIJFLAG#|

TIME#TIME4T

PLT#PLT+T
IF{PLT=PLTINC) 35, 33 33
CALL OUTPUT

PLT#Q

IF({TIME-TMAX 2437937
IF(LASTCC)4Ds 34540
K#K+ |

NI#0

PLT#NeN
IF(ICONTR)2I59I92I2
STOP

END

 
50~
$IBFTC QUT DECK
SUBROUTINE QUTPUT
C
C _
DIMENSICN A(60+s60)9sC(60960)sHP(60960) sQPT(E60+60) s DIMENS
IX(60)sY(60)s2(60)sXICI6G)sTQP(60) ‘ _ - DIMENS
c _
COMMON CoHPsAsQP ToX9Z Y s ITMAX sKKsLL s MM,
| JUFLAG O XICHNI s TIMEs TMAX s TZEROSNESTQP, T
211Z5sICONTRSPLTINCIMATYESSICSS s JFLAGSPLT
C
IF(NI)29s 192
I NI#I
NC#10 , | | (
DO 11 NCM#I+51410

WRITE(Es200) LLs ((A(IsJ)sJENCMsNC) »I#1sNE)
200 FORMAT (2HQAsI2/(1H SIPIBE!1+2))

IF(NE=NC) 10s!0s1! '
I NCANC+10

10 NC#ID
DC 21 NCM#!1 51,410 : ' '
- WRITE(S59201) ((CUI9sJ) s JENCMeNC) 9 I#I14NE)
201 FORMAT (2HOC/(IH S IPIDEI142))
[FINE-NC) 20920s21
21  NC#NC+I10 '

20  NC#I10 -
DO 31 NCME1s51,10
WRITE(69202) ((HP(IsJ) s JENCMsNC) s I# 1 sNE)
202 FORMAT (3HOHP/(1H sIPIOEI143)) |
IF(NE=NC) 252531
31 NCHNC+10

2 WRITE(65203) TIMEs (X(I)sI#IsNE) -
203 FORMAT(4H T #,IPE(Qe2s4H X #s /(IH $5Xs 10E1143))
IF(JFLAGsNEL5) GO TO 30

WRITE(65204) (Z(1)sI#1sNE)
204 FORMAT(6HOZ # »IPIOEI143/(1H 3»5XsICEI143))
30 RETURN
END

 
~51-

$IBFTC SsuBz DECK

aNaNE

SUBROUTINE DISTRB
DISTRB FOR REPCRT EXAMPLE

DIMENSION A(60+60)9C(60s60)sHP(£60+60) sQPT(60s60)
IX(60)sY(60)sZ2(60)sXICI(eD0)sTQP(6D)

COMMON CoHP sAsQPToXsZ oY s ITMAX sKKoLL sMMy

| JUFLAGsXICoNI o TIMEs TMAX s TZEROSNE s TQP 4T,
211 ZsICONTRSPLTINCsMATYES s ICSS 3 JFLAGSPLT

TXH#TIME+DeS5#T

Z{IY#SIN (2.0%*TX)
RETURN
END

$IBFTC DSENS DECK

C

SUBROUTINE DISTRB

DISTRB FCR TIME RESPONSE SENSITIVITY CALCS,
DIMENSION A(60+60)9sC(60+60)sHP(60s60)sQPT(60+60) s
IX(60)sY(60)sZ(60)sXIC(60)sTQP (60}

COMMON CsHPsAsQPTsX9ZsY s ITMAX 3KK s LL y MM,

| JJFLAGsXICoNI s TIME s TMAX s TZEROSNESTQP 5T
21 1ZsICONTR4yPLTINCYMATYESsICSSsJFLAGSPLT
DIMENSION IR(5)sIS(15)9JS(15)sIQ(30)sXT(5,1000)
I XSEN(15530) sXPSI(30) |
IF(NI)lsls2

IF(ICONTR+2)5 4443

IF(ICONTR+2)7 646

INITIAL INPUTS AND CALCS.

READ (5 100)(IS(I)sJS(I)sI#1395)sNTIsNSENS
FORMAT(6(21344X))

NDT#1 '

ICONTR#=2

NTIMO#NTI—1

DO 8 I#1sNE

Z(1)#0.0

DURING SOLUTION OF SYSTEM EQUATIONS

DC 20 I#!sNSENS

ICO#JS(1) -

XT(IsNDT)#X(ICO)

NDTH#NDT+ |

GO TO 30

JUST AFTER SYSTEM SOLUTION IS COMPLETED
IST#N

ICONTR#-3

DO 21 I#1sNSENS

DO 21 J#I 4NTIMO

XTI o J)#DaB* (XTI J)+XT(I sJ+11))

XT # AVG VALUES OF SENSITIVITY EQN INPUTS
WRITE(64+102)
FORMAT (3HOXT/(1H 4I10EI]e3))

AFTER COMPLETING EACH SENSITIVITY RUN -
ISTHIST+} | - .
IF(IST-NSENS)31s31932 '

((XT(IoJd) s JHEI sNTI) s I#I1 +NSENS)

DIMENS
DIMENS

29880105
29880107
298840108

29880113
29880115
29880117

29880123

29880201

298802172

29880203
29880205

29880209
2988021 1
29880213
2988U2 1 4

29880215
29880217
32

31

101

4 |

30

- -52-

GO TO NEXT CASE

[CONTR#D

PLTINC#TMAX

TMAX#De0O

NI#I

GO To 30

I 1Z#ISCIST)

COLe I1Z OF HP MATRIX MULT. BY 2z
WRITE(6s101) : :  IS({IST)sJSI(IST)
FORMAT(IBHUSENSITIVITY T0O A(9139|H9’I39IH))
TIME#TZERO

NDT#I

DO 41 I#I sNE

X(1)#0e0
Z(1)#0.0
JJIFLAG#O

"DURING EACH SENSITIVITY RUN -

Z(IlZ)#XT(IST;NDT)
NDT#H#NDT+ |

RETURN

END

29880219
29880221

29880301
29880303

29880305

29880309

29880315 °

29880317
" $IBFTC SUBV DECK

OO O OOy MM

|2
30

SUBROUTINE VARCO(XTR,TX)

FOR USE WITH DISTRB AND MATEXP FCR
VARIABLE Z-Se GIVES {-ST ORDER EXTRAP.

FOR AVGe X AND TIMEs PLUS RESTART
ON 1-ST INTERVAL. DISTRB FORM #

CALCe MATRIX COEFFe=S}
CALL VARCO(XTRsTX)

ETC.

IF NI#0

CALCe Z-S USING XTR(I)=S AND TX (TfME);

DIMENSION A(60,60)C(60+60)sHP (60 6D),OPT(6u,60),
IX(60)sY(60)sZ(60)sXICI60),TQP(60)

COMMON CsHP sAsQP T oX3Z s Ys ITMAX sKKsLL 4 MM,
F JUFLAG o XICeNI g TIMEs TMAX s TZEROSNEsTQP 4T

DIMENSION XTR(60)sXL(60)

IF(NI)1 sl e2

FIRST ENTRY

NV # |

TX#TZERO+O 5% T

DO 10 I#I| 4NE
XTROIVH#XIC(I)

GO Tn 30

IF(NV)3s3s4

SECOND ENTRY

NV #0

TIME#TZERO

PLT#0eDN

DO |1 I#IsNE
XLII)#XICI(I)
XTROI)#0e 5% {XLIT)+X (1))
X{I)#XIC(I)

Go To 30

ENTRIES AFTER SECOND
TX#TIME+005*T

DO 12 I#I|sNE
XTROII#X(I)+0e 5% (X(I)=XL(I))
XL{T)y#X(I)

RETURN

END

21 1ZyICONTRyPLTINCIMATYES s ICSSsJFLAG,PLT

25880101
29880)03
29880105

- 29880107

29880109
2988011 |
29880113
29880115
29880117
DIMENS .

DIMENS

29884118

29880120
29880121
29880122
29880124
29880202
29880204
293880206
29880208
29880210
29880212
29880214

29880216
29880218
29880220
29880222
29880224
29880301
29880303
29880305
29880307
29880309
2988031 |
29880313
,_5)4__

$IBFTC FGEN CECK

N OO OOy O Oy O M

100

-
101
86

SUBROUTINE DFG(XD,s2D)

EQUIVALENT TO 8 DFG—-S WITH UP TO 32
POINTS EACHes CALLED BY DISTRB.

INPUTS ARE
NDFGS NOe OF DFG-S5 USED
NPTS NOe OF POINTS IN EACH DFG
XP INDEPENDENT VARIABLE DFG POINTS
ZP DEPENDENT VARIABLE DFG POINTS

XD IS THE INPUT VARTABLE AND ZD THE QUTPUT

CIMENSION A(60+60)sC(60+60)sHP(60s60) sQPT (605960}
IX(6D0Y Y (60)sZ2160)sXICL60)sTQP(60)

COMMON CoHP sAsQPT9X9Z oY s ITMAX 9KKsLL s MM,

| JJFLAG e XICoyNI 9 TIME s TMAXSsTZEROSNEsTQP s T
211 Z9sJCONTRIPLTINCIMATYES s ICSSyJFLAGYPLT

DIMENSICN XP(32+8)9ZP(3298)9SL(32+8)sNPTS(8)
|JP(8)sZ2D{(8)+XD(8)

IFINIY 1921
FIRST CALL COMP.

READ (5,100)
FORMAT(1258X5813)
DO 86 I#!4NDFGS
NP#NPTS (1)

READ (5s101)
FORMAT (8E 1043}
WRITE (65,200)

NDFGSsNPTS -

(XP(JsI)sZP(Js1) sJ#IsNP)

Ia(XPUJeIYsZP(Jsl) s J#HI 4NP)

2000FORMAT (4HQOCFGs 13,1 7H XP AND ZP IRNPUTS/

e

|18
D
13

| (ITHO 94 (2E12e494X)) )

DO 3 I#!1 sNDFGS

MENPTS (1)

DO 3 J#I M '

SLIJs IV#(ZP(J+1 91 )=ZP(Js 1))/ (XPLJ+1s1)-XP(Js1))

DO 5 I#!| sNDFGS

DO &4 J#2932
IF(XD(I)=XP({Js1))59544
CONTINUE

JP(I)#J

CALCSe MADE EACH TIME

DO 6 I#! sNDFGS

J#JIP (1)
IF(XD{I)I=XP(JsI))1Ds!llsl2
IFI(XDUI)=XP{J=1s1))13414415
J#J—| -
IF(J=1)16s16s10

J#2

GO To 19

ZD(I)HZP(J=1y1)

GO To 6

JH#I+ |

IF(NPTS(I)=J) 17518518

29880105 .

29880106
2988007
29880108
29880109
129880112
29880113
25880110
298801 14
29880115
29880116
DIMENS
DIMENS
25880117
258806118
29880119
29880121 .
29880122
29880123

29880124 -

29880125
29880201

29880202

29880204
29880205

29880207
29880208
29880209
2988U210
2988021 |
25880212
29880213

2988021414 .

29880215
29880216

29880218
29880219
29880220
-55-

JHNPTS(1)

GO TO 19
ZDUIVHZP (U 1)
GO TO 6

WRITE (6s102) 1

FORMAT (4HODFGs I3 16H RANGE EXCEEDED.)

ZDUIYVH#ZP(J=1 s 1) +SLIJ=1 gT)#(XD(IV=XP(J=1,1)).
JP(I) STORES VALUE OF XD LOCATION

TO USE AS FIRST TRY NEXT TIME.
JP LI #J - ,

RE TURN
END

29880222
29880223

29880224
29880225

29880301
29880302
298803303
29880304
$IBFTC TRLAG

OOV OO OO0 OO O v O

N OON

N

Oy Oy

Oy Oy ) Oy

23

24

22

ZD

®%#¥t¥% NOTE -

33
31

26

56~

~ DECK
SUBROUTINE TRLG(XTsWsZT)
VARIABLE TRANSPORT LAG GENERATOR = FORTRAN IV

USES UP TO 300 POINT APPROXIMATION
UP TO 6 VARIABLES.

FOR
USES INVENTORY CALC,.

INPUTS FOR EACH LAG (TOTAL # NLAGS)
e INPUT FUNCTION XT{(1I)
2e MASS FLOWRATE W(I)
3. INITIAL VALUE OF LAG TIME TI(I) . ,
4o MINIMUM EXPECTED VALUE OF MASS FLOW:WMIN(I)'

OUTPUTS ARE LAGGED FUNCTIONS ZT(I)

DIMENSION A(60+60)»C(60s60)sHP(60s60) sQPT(60+60)

IX(60)sY (60),Z2(60)sXICLED)»TQP(60D)

COMMON CoHP sAsQPTsX9ZsY s ITMAX yKKsLL sMM,

[ JJFLAG 9 X1ICsNI s TIMEs TMAXsTZEROSNESTQP 5T
211ZsICONTRsPLTINCyMATYESs1CSSyJFLAGSPLT

DIMENSION XT(6)9W(6)9TI(6)9WMIN(6)’ZT(6)9XS(BDD’6)9

IPS{3006) sKT(6) s JT(6) s XIMP(6E) 9 JMP(6) s NIMP (6)

NI # |-ST CALL FLAG (# 0 ON
T # COMPUTATION TIME INTERVAL

| -ST CALL)

IF(NI)20s21 20

FIRST CALL COMP.
READ(5,100) NLAGSsTI 4 WMIN
FORMAT(I2/(6E1De3))
WRITE(6s101) TIsWMIN

FORMAT (26HQTRLG INPUTS - TI AND WMIN/(IHDs6EI1845))

DO 22 1#!sNLAGS
XIMP(I)#1e0
XSUIsT)#XT(I)
PSCIsI)#W(II*TICI)
XNSPH#PS (1 s1)/ (WMIN(I)*T)
DO 23 M#1,10 |
PI#XJIMP (L) %XNSP
IF(300e0-P11)23424424
XJIMP I )#XIMP(1)+1.0 -

JMPUIYHIFIXI(XJIMP (1))

CKT(I)#2

JTCI)#]

CNJMP (1) #1

NV#-1

CALCS. MADE EACH TIME

NVAENV+ | | |
IF A RESTART FEATURE IS USED (WHERE THE INITIAL TIME
STEP CALCULATION IS REPEATED),
OMIT THE TRLG CALCe THIS
REMOVING STATEMENT 33
IF(NV)31932,31

DO 17 I#1sNLAGS
IF(NJMP(I)—JMP(I))26,27,27
NIMP (T #NJUMP (1) +1

| -ST CALL OMISSION MAY BE DELETED BY

 

THE FLAG NV AND STATEMENT 33 WILL

25880105
29880106

29880108
29880109
29880110
298801 1§ |
29880112

298801 14

DIMENS
DIMENS

DIMENS
DIMENS

29880121 -
29880123

DIMENS
DIMENS
29880202
29880203
29880204
29880206
DIMENS
29880209

29880212
29880213

29880214

29880216.

29880218
29880219
27

@ O

=57~

GO TO 17

NJMP (1) #1

KEKT (1)

JHEJTIT)

XSIKsT)#XT(I)

PS(Ks I )EXIMP (L) *¥W (I ) %®T

J#NCe OF ELEMENT AT EXITe K#NO. AT ENTRANCE
IF(PS(Js1)=-PS(KsI)) | s2s3 |
ZTUIYH#XS(Js 1)

IF(J~300)1697s7

JT (1) #]

GO TO 30

JTCIV#I+]

GO TO 30

COLLT#XS(Js1)
COLLP#PS(Jy1)
DO I5M#1+300
IF(J-300)85955

J#0

J#J+]
PQ#COLLP+PS(Js1)

IF(PQ-PS{KsI)) 1191213
CCLLTH(COLLT*COLLP+XS(JsII®PS(JsI1}))/PQ

COLLP#COLLP+PS(Js1)
ZT(1)#(COLLT*COLLP+XS(J,1)%¥PS(Js1)1/PQ

IF(J=-300) 14516416
JT(IY#1

GO To 30
JTII#J+

GO To 30

PS(Js1)#PQ-PS(KsI)
ZT(I)#(COLLT*COLLP+XS(Js 1) ¥PS(Js1) )/ (COLLP+PS(Js 1))
JT LI #J) | '

GO To 30

ZTUIYVH#XS (D)
PS{Js ) #PS(Js 1) -PS{K,I)

[IF({K=300)44+545
KT(I)#!

GO TO 1|7

KT (D) #K+]
CONTINUE

RE TURN
END

29880220
29880221
29880222
29880223
29880224

29880301
29880302
29880303
DIMENS.

25880305
29880306
29880307
29880308
29880309
298830310
2988031 |
DIMENS

DIMENS

29880316

29880319
29880320

DIMENS

29880401

29880402
29880403

29880404
29880405

29880407
29880408
29880409
298806410

29880412
25880413

DIMENS

29880416
29880417
25880418
29880419
 

 

 

 

aaete

   

 

s
I

 

 

 
1-30.
31.
32-61.
62 .
63.
6L,
65 .
66.
67.
68.

70.
72%

81.

ok.

95.
96.

OWoHEoEaoaggarHInIi ool raEqgd RN oy B nE D

-*

*

TH<4GpqUbbQQUqg<<HoogEzzgngddagqEdgderm g ggE R

59

INTERNAL DISTRIBUTION

. Adams

Anthony (K-25)

. Blalock
. Borkowski

Bullock
Burke
Clark

. Dandl
. Danforth

Davidson
Ditto

. Duggins

Engel
Epler

. Farris

Fowler (Y-12)
Fowler

Fry

Hagen
Harrill

. Hyland:
. Jamieson

Jordan

. Kaye

Kedl -

. Kress
. Kryter
. Lawson

125-
Lucius

Mankin, Jr. (K-25)128-

. Martin
. McClain (XK-25)
. McLain

Metz
Miskell (Y-12)

. Moore
. Mossman

135-

97.

8.

99.
100,
101.
102.
103.
104,
105.
106,

. 107.

108.
109,
110,
111.
112.
113.
11k,

115.

116.

“117.

118.
119.
120.
121.
122,
123.
124,
126.
127.
132,
133.
134.
149,

150.

2RO namsE

L.

D EqgOHGgEqgNHOEZ Q"D OO QO

<HtHEHEHORDWOQKHKE DN

M' M.

ORNL-TM-1933

. Nestor

Nowlin

. O'Brien

Olsen

. Partain
. Patten
. Peele

. Perry

Pitkanen
Prince

. Roux
. Sadowski

Sartory
Smith
Smith

Squires

Stone

. Tallackson
. Tobias
- Trinko
. Tucker
. Walker

Watson

. Whatley
. Whitt

Wichner
Wilson
Yarosh

Central Research Library

Document Reference.Section

Laboratory Records Department

Laboratory Records, ORNL R.C.

ORNL Patent Office

Division of Technical Infor-
mation Ext.

Research and Development
Division ORO
151.
152.
153.

15k,
155,
156.
157,
158.
159,
160.
161.
162,
163.
164,
165.
166.
167 .
168.
169.
170.
171.
172.

J
J
R.
E
S

60

ORNL-TM-1933
EXTERNAL DISTRIBUTION

S. G. Bloom, Battelle Memorial Institute, Columbus, Ohio

Neal Carter, Battelle Northwest Labs., Richland, Washington

E. H. Cocke-Yarborough, Electronics Division, AERE,
Harwell, England |

S. J. Gage, University of Texas

R. P. Gardner, North Carolina State University
B. E. Gilliland, Clemson University

S. H. Hanauer, University of Tennessee

T. H. M. Hung, Deere & Company, Moline, Illinois
J. M. Jansen, Jr., University of Illinois

T. W. Kerlin, University of Tennessee

‘L. B. Koppel, Purdue University

J. B. H. Kuper, Brookhaven National Laboratory

Rufus Oldenburger, Purdue University

H. M. Paynter, Massachusetts Institute of Technology
W. Prados, University of Tennessee '

. C. Robinson, University of Tennessee
¥, Saxe, North Carolina State University
Siddall, AECL, Chalk River, Ontario

. E., Stephenson, University of Arkansas

Otis Updike, University of Virginia

T. J. Williams, Purdue University

W. C. Wright, University of Tennessee

 

:1n