- -'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