Robert H. Schuler Cornegie-Mellon University Pittsburgh, 15213
I
(
~altulator-plotterRoutines for Numerical Integration of Kinetic Expressions
The kinetic expressions which describe other than the most simple of chemical problems usually involve differential equations which cannot he integrated to give an analytical relation between the variables. With presently available computers these cases can, of course, he treated by numerical methods but for many students and researchers there is undoubtedly a certain reticence toward such an approach because of the required programming, particularly for cases where only a limited amount of information is desired. During the past several years this author has treated a variety of kinetic problems with a numerical integration routine built around a programmable calculator-plotter combination. The main program allows convenient entry of the appropriate differential equation as a subroutine with the kinetic parameters heing given in this subroutine as general constants which can be specified numerically when calling the main program. Because the kinetic expression can be written into the calculator in what is basically its algebraic form, the equation can he programmed with relative ease and for cases of moderate complexity such programming usually does not hecome excessively involved or time consuming. We will illustrate several examples below. More importantly, in actual operation the results are plotted as they are calculated so that the student gets the feeling that he is following the simulated reaction as it progresses and is ahle to interact with the calculations as the more important kinetic features become apparent. It has been found that in practice such direct interaction provides insight into the critical features that can affect the course of a given reaction and that one can frequently proceed iteratively to fit the kinetics to a given problem or experimental situation. Such an approach is of obvious value to the practicing kineticist who is faced with problems which, as indicated above, are frequently mathematically complex. The pedagogic value has been demonstrated by students a t both the graduate and advanced undergraduate levels who, it has been shown, can readily master the relatively simple programming requirements of the calculator even though they have had little or no prior experience along these lines. We describe below first the requirements of equipment and an appropriate master program. These descriptions are then followed by several examples of the writing of kinetic subroutines. While attention has been focused on kinetic examples it should be pointed out that the integration routine used is general for any integral of the form of
where f'(x,y) can be expressed in its algebraic form. With appropriate modifications it is, therefore, potentially useful in other applications. What is done here is simply to put a standard approach using numerical integration methods into a format which minimizes the thinking required to carry out the actual mechanics of the integrations. The student is thus allowed to be concerned priSupported in part by the U S . Atomic Energy Commission. 166
/ Journal of Chemical Education
marily with the differential equation and, of course, the results. Hardware The initial calculator-plotter configuration used by this writer was acquired some four years ago and consists of a programmahle Hewlett-Packard 9100A calculator with a 9101A extended memory, a 9125A plotter and 9120A printer as peripheral equipment. The printer is very useful hut not essential. The instruction set consists of 64 instructions coded by two digit octal numbers. They are keyed in their algehraic form hut referred to in the display and program listings by their octal names. The currently availablk equivalent equipment with a somewhat extended program language is a Hewlett-Packard 9810A calculatorhaGing 1012-memory steps together with a 9862A plotter. An alphanumeric printer is optional. A more advanced instrument is the Hewlett-Packard 9820A calculator which possesses many features that, together with a 9862A plotter, make it a more desirable package. The instruction set of this latter instrument is considerably larger than that of the 9100A and a number of editing functions make its programming relatively easy. Alphanumeric display and printout are included and listings are in the programming language itself so that they are readily readable. The basic 9820A calculator comes with 173 memory registers which is adequate for a somewhat abbreviated fonn of the program discussed here hut for more general use it is recommended that it be acquired with 429 memory registers. Basically one needs a programmahle calculator or minicomputer which is capable of heing interfaced with an X-Y recorder or oscilloscopic display. The instrument should he ahle to store a t least 11 numerical quantities needed in the integration routines plus any additional quantities needed for initial conditions and parameters in the actual examples. An absolute minimum of 20 numerical storage registers is recommended but one should preferably have more. In addition the calculator should he capable of storing the required program information which, depending on the calculator language, will usually amount to at least several hundred program steps. (The author's program for the Hewlett-Packard 9100A calculator is a rather basic program and requires 23 registers (322 program steps); that written for the 9820A instrument is considerably more elaborate and requires 198 registers (-1500 program steps).) Programmable statements which ~ e r m i tloopine and hranchine are rewired. Calculations should, ofAco;rse, he carrieduout with sufficient enough accuracy (i.e. 8-10 significant figures) that rounding errors are unimportant, butmost desk cal'ulators possess this feature. Ease of programming and program editing is extremely important. Costwise the basic calculator can he expected to run -$5,000 with an additional ~ $ 3 , 0 0 0needed for the plotter and necessary interface. One can, by printing out the results, avoid or postpone the initial cost of a plotter but this form of output becomes relatively time consuming and sacrifices much of the insight available in the graphical display of the results of the calculations. In general this mode of operation cannot be recommended. This author has found that the cur-
rently available Hewlett-Packard 9800 calculator series satisfies the above requirements admirably. The Master Program' Each instructor will have his own ideas on the requirements of an appropriate master program so that this writer can only discuss what he feels are the more important aspects that should be considered and suggest a general approach. The software supplied by the manufacturer2 is usually written in such a general form that it is considerably less convenient than a program developed for a specific application. The master program should, of course, be written in such a way that the differential equation can be conveniently entered as a subroutine. For ease of operation the program should require a minimum of numerical entries. The essential entries are the graphical scales, the quantities defining the initial chemical composition of the system, and the parameters involved in the kinetic expression. Provision should be made to skip entries not needed or not required to he changed between calculations. It is convenient also to provide a means of altering the integration increment so that one can readily test the accuracy of the approximations used in the numerical apnmach bv carrvina . - out the calculations with a smaller or larger increment. (Rounding errors are usually unimportant hut significant errors can be introduced by the finite increments used in the summations, especially when the derivative varies rapidly with the independent variable.) Various embellishments on a minimum program can be attempted and may, in fact, be required by certain applications but one should avoid requiring excessive entries because these only distract the user from the benefits of the interactive mode of operation. Printout of unnecessary data similarly should be avoided since it is time consuming and a considerable distraction. The following illustrate some of the operating characteristics of the programs developed by this writer. The program for the 9100A calculator requires as its initial entry a definition of the X and Y scales of the plot. The magnitude of these scales in the desired coordinate units per inch of plot are entered into the X and Y registers of the calculator. If the differential equation has been appropriately entered and the initial conditions and parameters specified, then the integration can be set into motion with no further entries by pressing the instruction key CONTINUE. The master program calculates and stores the X and Y plot scaling coefficients and the magnitude of the increment in X required to plot ten points/ in., stores the initial values of X and Y in auxiliary registers for possible requirements of the differential equation, plots the initial value of Y at X = 0, and then calls the integration subroutine which in turn calls the differential equation as a subroutine. The integration routine uses a third-order Runge-Kutta method to calculate AY over the AY interval AX and outputs a new value of Yt+1 = Y, a t X,+1 = X, AX which is scaled and plotted by the master program. This process is then repeated until the program is terminated. In order to enter or change a parameter or initial value the calculator flag3 can he set when the scale entries are made, in which case CONTINUE will divert the program to an auxiliary suhroutine which allows one to rapidly examine the current values of the parameters and make any desired changes in the registers storing the various quantities. The first entry to this latter subroutine is the number of parameters to be examined. If the flag (which has been cleared) is set again a t this point the calculation increment and initial values can be changed. When modification of the parameters has been completed the program is returned to the master Droeram and the inteeration ~roceeds.The master program is labelled Po a i d reqiires eight storage registers and the integration subroutine is labelled Plo and requires
+
+
seven storage registers and can be called directly at any time to enter or change the parameters. The differential equation subroutine must he labelled Pz (since this is how it is named in PI) and if entered last can readily he changed without reentering the other programs. Programs Po, PI, and PIO are permanently stored on magnetic cards and program P z can either be stored on and entered from a magnetic card or from the keyboard. The above approach has been adapted and extended considerably to give a rather elaborate but more versatile master program for the 9820A calculator. In this case one can use the alphanumeric display capability to identify the required input and the large number of available flags to provide options for controlling the mode of operation and entering various types of data. The program, when it is first run after being loaded from a magnetic card, will cause the display immediately to show WRITE EQUATION. Pressing the keyboard edit command RECALL causes the line number of the differential equation subroutine to be displayed followed by "DY/DXH; SFG 9.4 At this point the user can either key in the differential equation or insert it from a magnetic card. Statements entered from the keyboard to describe the differential equation will be printed out automatically for a permanent record. With an appropriate equation in memory a start of the program will cause OPTIONS to be displayed at which point the desired options can be specified in terms of their numerical names shown in the table. The program sets the flags corresponding to these numbers and the flags are tested at the appropriate points in the program. After the desired options have been entered, pressing the continuation key RUN PROGRAM causes the options to be executed and the integration and plotting to begin. Termination is automatic at the right-hand limit of the plotter. Use of options in the way described here permits the user to skip unwanted entries. For example if one wishes to examine the effect caused by a change in a particular rate constant one obviously should not have to reenter the coordinate scales, the initial values, or other rate constants. In this case one specifies option 2 which calls the parameter subroutine which in turn allows a change to be made in the rate constant (or constants) of interest. The program is then automatically returned to the master program and the integration commences. No other entries will be required and all other quantities will retain the values used in the previous calculation. Minimizing the entries in this way encourages a very desirable interaction between the user and the calculations. The Integration Routines The master routine for the 9820A calculator, unless otherwise instructed, uses a Runge-Kutta subroutine for calculating the incremental values of Y. This procedure is 'Appropriate master programs for the Hewlett-Packard 9100A and 9820A configurations described in the text are available from the author as CMU Radiation Research Laboratories Special Repart 18. 2Hewlett-Packard Program 09100-70006 is, for example, a genera1 integration routine for use with the 9100A calculator. 3A flag is a special register which contains a zero or a one and can be set independently from the numerical entry. It can be tested during the running of the program to control branching depending on whether it is clear (0) ar set (1). The Hewlett-Packard 9810A calculator has a single flag and the 9820A calculator has 16 flagsof which four (flags0, 13, 14, and 15) serve special functions. 'In the language of the 9820A program, lines can be addressed either directly or by literal fields enclosed in quotes. The aubroutine for the differential equation has been named "DYJDX." When the oroeram is first loaded this subroutine contains a dummy .ta;emenr ahmh sets flap. 9 (SF&?> when the program is firzt run and causes the program t o stop at the immedietcl! pre. ceding h e and dcsplay the insrrurtmn \\'RITE EQL'ATIOS. Volume 52. Number 3.March 1975
/
167
somewhat time consuming in that it calls the differential equation four times (see below) so that if the derivative does not change significantly over the time increment a simple summation of the form
may he adequate5 and may be carried out by specifying option 6. The simplicity of the 9820A language is illustrated by the fact that the suhroutine for a summation of the form of eqn. (1) (which here is called "INTI") requires only two lines, i.e.
(GSB designates a subroutine call and RET a return to the calling program.) The first line increments X by AX12 to calculate the derivative a t the midpoint of the interval (AX is stored in C) and then calls "DY/DX" which computes the derivative and stores it in register Z. The second line increments Y by d,X(d Y/dX) and X by a n additional AX12 and then returns operation to the main program which draws a line to the coordinate indicated by the appropriately scaled incremented values of X and Y. The Runge-Kutta suhroutine ("INTG"), which is used if option 6 is not specified, is quite straightforward and using auxiliary registers R60 through R63 for storage of intermediate values of Y can he written in five lines
CZ
"TNTG"; Y R61; Ci2 X
-CZ
CZ
R60
+
-
+
+
-+ -+
R62; R62i2 R63; C/2 X
+
R6O; GSB "DYIDX X; R61/2 R60 Y; GSB "DYIDX R60 Y; GSB "DYIDX" X, R63 R60 + Y; GSB "DYIDX 2R63 CZ)/6 Y;RET
+ (R51 + 2R62 +
+
-
It will be noted that in this latter case the differential equation is called four times so that if any auxiliary quantity is incremented during the computation of the derivative this fact must be taken into account. The best procedure is to insert an appropriate subroutine or subroutine call between the first two lines of the above so that the same programming of the differential equation can he used for hoth the simple summation and Runge-Kutta methods of calculation. The ease with which the 9820A programs can be edited becomes apparent here. The Runge-Kutta routine is valuable for equations of the form Y' = P(X,Y) where Y varies significantly over the increment AX. It offers no advantage over the simple integration if the derivative is not a function of Y hut this is not the usual situation in kinetic examples. Writing the Differential Equation One must, of course, translate the appropriate kinetic expression into the language of the calculator being used. The user must first become familiar with the details of the particular language involved and the following is given only to illustrate the approach with the Hewlett-Packard 9820A instrument. Let us take as a simple example the differential equation describing the disappearance of an intermediate present a t concentration [C] a t time t when first- and second-order processes are hoth important,G i.e.
The language of the 9820A instrument very much resembles algebraic notation so that programming is extremely straightforward. In this case the variables t and [C] in 168
/
Journal of Chemicai Education
Figure 1. Results obtained from integration of the kinetic expression describing the simultaneous disappearance of a reactant by first- and second-order processes (eqn. ( 2 ) ) . Curve 1AJ for k , = 1.0 s - ' , k2 = 0: Curve ( 8 ) tar k , = 0. k 2 = 1.0 M-' s - ' and Curve (CJ for k , = 0.5 ssl, k 2 = 0.5 M-' 5 ' . Calculations are made at 0.02-s intervals. These results are essentially identical to those obtained by plotting the integral 01 eqn. 12) (see footnote 6 ) directly. With the R~nge-Kuttaroutine deviations are observed only if the calculation interval exceeds 0.5 s whereas the simple summation produces comparable deviations tor calculations at 0.02-s increments.
eqn. (2) are simply redefined as X and Y and stored in the resisters of the same name and the two constants are stored in registers R1 and R2. The subroutine can then be written in its algebraic format and should appear as "DYIDX;
-
RIY
-
R2YY
-
Z
The simolicitv of ooeration in this case is auparent. The language of the 9 1 0 0 ~is more primitive and;kquires considerablv more attention to detail.7 In both of these cases the language is unique to the particular calculator but other instruments may employ more conventional languages such as Basic. 5The overall speed for the Hewlett-Packard 9100A calculator9125A plotter combination is largely limited by the time required for plotting (-1 s/point) so that a Runge-Kutta integration procedure, which allows one to calculate a smaller number of points for a given accuracy is usually preferred. Alternatively one could calculate points more frequently and plot only a fraction of the points (i.e. every 5th or 10th point) but the Runge-Kutta approach is more efficient. 61n this case the equation can be integrated to give [C] = flt), i.e. [C] = [Clo~-""ll+ (k~/kdCo(1- e - R " ) ] l . ?In the master program written for the 9100A instrument the current value of [C] is stored in register E and the parameters k~ and kz can readily be stored in registers 201 and 202 by subroutine Plo. Programming of eqn. (2) involves recalling [C] from register E, multiplying the quantity found there by the parameter stored in register 201, squaring [C] and multiplying it by the parameter stored in register 202, summing the two terms, negating the sum, and leaving the result in register Y. The following keystrokes can be used to do this E; i ; X; f ; 2; 0; 1; FORMAT; r ; X, 1 ; s ; ' ; 2; 0;2; FORMAT; n;X; 1 ; +; 1 ;CHANGE SIGN; t ;FORMAT; END. This subroutine is then stored as program Pz. While this instruction sequence looks formidable it is quite straightforward once one has even an elementary feeling for the calculator language and inspection will show that it quite directly carries out the algebraic operations indicated on the right-hand side of eqn. (2). (The instruction FORMAT, a, far example, calls into the X register the value stored in the register named in the immediately previous numerical entry-register 201 in the first instance-and the instruction 7 duplicates this value in the Y register. The terminating instructing FORMAT; END acts as a subroutine return.) As mentioned above it has been found that students when confronted with a specific problem can pick up this language quite rapidly.
Several plots of the integral of eqn. (2) obtained by these numerical methods are illustrated in Figure 1. Most examples are, of course, more complicated and one is frequently confronted with a problem in which additional information is required to define the quantities in the differential equation in terms of either the dependent or independent variahles. The only requirement is that these definitions can be put into a mathematical form that can be adequately handled by the calculator. In certain instances a quantity within the differential equation must itself be stated in the form of an integral, in which case an appropriate subroutine must be built into.the procedure to increment the quantity properly. In many kinetic examples the rate of reaction depends on the product and reactant concentrations in such a way that in the differential equation the dependent and independent variables (usually time or extent of reaction in the latter case) are not separable and the integral cannot be obtained in closed form. A very gwd practical example of such a situation is the seemingly simple case describing the production of Iz when HI reacts with radicals produced at a constant rate. This example has previously been treated by Mani and Hanrahans who used an IBM 709 computer together with numerical methods related to those described here. In this case radicals produced at a rate P Reactants
P
R.
Figure 2. Production of l l z j by reaction of radicals with unit concentration Hi. Rewits Obtained by numerical integration of eqn. ( 7 ) as described in the text. Values ol k s / k r are far ( A ) 0.01. ( 8 ) 0.1. I C ) 0.5. ( D l 1 . (El 1.9. and I F ) 10. Abscissa is given in terms of the equivalents of total reaction. Calculations are made at intervals corresponding to 0.02 units of reaction. As is indicated by the solid curve the caicuiation for k 8 / k a = 0.02 (Curve A ) becomes invalid for reactions greater than 1.0 units because the residual [ H I ] in the caiculation becomes negative. A vaiid calculation can be made by using a smaller interval ar indicated by the dashed curve (calculated at intervais 0.002 units) or by branching as described in the text. Curve E describes the iodine production curve in the 1.9. actualexperiment (seereference (9))showing that k 6 / k 4
-
(3)
are consumed by rapid reaction with HI
so that the initial rate for iodine production will be 'hP. However, as iodine builds up a complication soon sets in because the radicals also react with the product iodine with a rate constant similar to that for eqn. (4).
The differential equation descrihing this kinetic scheme is
Figure 3 . Production of I , far unit concentration HI soiutians initially containing added I , for k b / k n = 1.9. Coordinates are as per Figure 2. Note that the curve tor [ I 2 ] = 0.5 has very nearly zero initial dope as expected lor the approximately equal importance of reactions (4) and ( 6 ) under these Conditions. One can readiiy see the pronounced effect of 1 , on the initial slope but the very srnaii effect on the slope near to the completion of the reaction.
where the stoichiometric equation
can be used to relate [HI] to the variables [Iz] and t. Even though both the chemical system and the differential equation appear to be extremely simple it turns out that, except for the special case where ka = ks, the variahles in eqn. (7) are not separable and the integration cannot he performed to give an algebraic expression for [I21 as a function of time.9 For use with the 9820A master routine the differential equation can be written as
-
-
8Mani. I.. and Hanrahan. R. d., J. Phys. Chem.. 70, 2233
(1966). %ee Perner, D., and Schuler, R. H., J Phys. Chem., 70, 2224 (19661. While [Is] cannot be given directly as a function of time, for the ease where the iodine concentration is initially zero, a change of variables allows eqn. (8) to be integrated to give both [I.] and t as a function of the new variable. These variables can, therefore, be related to each other though not explicitly in algebraic farm. The numerical approach illustrated here is more generally useful and one can easily treat cases where iodine is initially present as is illustrated in Figure 3.
(with kz/kl stored in R1, [HI10 in R2 and P in R3; FzJo will automatically be stored in R42). The first statement computes HI as [HIIo + [I210 - [I] - Pt/2 and temporarily stores the result in register Z. The second computes ikz/ kl) [12]/[HI] and the third uses this result in eqn. (7). The three statements can be written on separate lines as indicated or, if desired, all can he written on one line. Graphical output based on this program is illustrated for various values of kz/kl in Figure 2 and for various initial values of in Figure 3. These results are very similar to those of Mani and Hanrahan.8 When the RungeKutta subroutine is used and values calculated in 0.1-in. increments, each curve of 150 points takes approximately 100 s to calculate and plot. The calculations are accurate except for cases where kz/kl is very small so that at t = [HIIo/P one approaches the discontinuity indicated by the apex of the triangle in Figure 2. At this point the residual [HI] calculated in line 1 above may he insufficient for the Volume 52, Number 3, March 1975
/
169
requirements of the next increment so that the calculated value of [HI] becomes negative. As is shown in the figure the suhsequent calculations, of course, then become invalid. By inserting the statement
between lines one and two of the above program a test for a zero or negative value of [HI] will be made and, if this condition is found, the calculation will branch to the differential equation
which describes the suhsequent kinetics since a t this point effectively no HI is left to react. By adding the statement IF R3X > ZR2
+ R41);GTO " E N D
as the second line of the subroutine a test will he made to determine when Pt becomes greater than 2([HI]o + [Izlo), i.e. when all of the HI and Iz have been consumed, a t which point the program will automatically terminate. These additions indicate the power availahle for control of the program flow and the ease with which this power can be applied. In adding tests such as these one must, of course, sacrifice a certain speed of calculation hut the calculations are, in general, so rapid that the added time is usually not serious. For all except very complex situations 150 points along the curve (15 in. a t 10 points/in.) are calculated and plotted with the 9820A calculator plotter in under 2 min. The 9100 calculator which uses the slower 9125A plottertakes -2 min. longer. One can, of course, go on with many more examples but the above illustrate the general approach and some of the more salient features of the programming. The power is obvious. We are coming to the point in time where calculation methods such as those described here are becoming commonplace and students should be aware of the characteristics of the available methods. The approach described here stresses simplicity of programming and maximization of interaction between the user and the calculations. One cannot stress too strongly the importance of this latter noint since it is ouicklv found that the total immersion in the calculations produced by an interactive mode of operation emphasizes the relative importance of the different
170
/ Journal ot Chemical Education
Calla subroutine "EQUA." h Pressinx edit command RECALL recalls the differential equation and prints its address and the current erpreJsion. If one desire= one can then rewrite the expression with the revision also being printed. The plot oriein ia initialized as (0.0) and X and Y scales are set to 1
issettoo. Y(01 i ~ s t o r e d i n R l 4 2 ) . ~ Allowa initial valves of X and Y to be entered. X(0) is stored in R(41) and Y(01 in R1421.r Allows X and Y scales to be entered. Origin is set s t (0.0). If X scale is set eaual to 0 on first pegs then "on-zero origin can be entered.d Instrncte program to use a simple summation in the integration. If this option is not specified a Runge-Kutta subroutine will be
tion. Allowa a ehsnee of inteeration increment,
or the mastel
rd calculator written for a ~ e w l e t t - ~ a c k a9820A
with 9862A plotter. Options 10 through 12 are open and can be added at the discretion of the user. Options can be entered in any order but are tnted in t h e order 1, 5, 9, 4, 3, 2, 8, 7, 6. Keying RUN PROGRAM with no entry causes program to exit from option selection loop. Integration and plotting automatically commences after the data required by all selected options have
option 4 takes precdence and 3 is not =muted. d ~thisoption f is exercised the d o t origin isset to (0.0) and X and Y scales to 1unit/in. Other plot parsmetera must be reentered. mmber of mathematical .For a differential equation requiring operatiom the ~ ~ ~ ~ e - ~ takes ~ t t a 50 s to caleulate summation specified by option 6 requires snd plot 150 points. he only 15 S. Eseh of these timea is increased by approximately 1%for each mathematical operation in the equation and by about 10% for each conditional test.
kinetic quantities much more strongly than is possible when one has large volumes of computer output to examine a t some time when one is not so vitally interested in the calculations. For this purpose the computing device is best dedicated to a single user and presently available pmgrammable calculators are both adequate and very suitable. They are in a cost range that is not too unreasonable for most institutions. The master programs written by the author, appropriately annotated as to the details of operation, are available upon request.'
