- Email: [email protected]
Subroutines for integration of systems of first order ODE's
Computer Physics Communications North-Holland. Amsterdam
48 (1988) 293-304
293
SUBROUTINES FOR INTEGRATION OF SYSTEMS OF FIRST ORDER ODE’s G. DELIC
and S.M. MALHERBE
Centre for Nonlinear Studies and Department Johannesburg, South Africa Received
of Physics, University of the Witwatersrand,
WITS 2050,
1 July 1987
PROGRAM SUMMARY Title of program:
RKPC
Peripherals used: card reader,
Catalogue number: ABBH
No. of cards in combined program and test deck: 617
Program obtainable from: CPC Program Library, versity of Belfast, N. Ireland (see application issue) Computer: IBM 3083524; Installation: Witwatersrand Computer Centre Operating system: VM/SP Programming language LANGLVL (77)
line printer
HP0
Queen’s Uniform in this
University
of
the
Nature of the problem Numerical integration of the vector equation Method of solution Application of second or predictor-corrector
4.2
used: FORTRAN
Keywords: Runge-Kutta, predictor-corrector, ordinary differential equations
VS LEVEL
1.4.1,
High speed storage required: 20272 bytes for OPT(O) compile No. of bits in a byte: 8
OOlO-4655/88/$03.50 0 Elsevier Science Publishers (North-Holland Physics Publishing Division)
first
dy/dx
order
= f( x, y).
to sixth order, fixed step, Runge-Kutta integration formulae.
Restrictions on the complexity of the problem Vector dimension set for 20 components. Typical running times Dependent on required range and accuracy: is in the range 0.14 to 0.18 s.
B.V.
for the test cases it
294
G. D&c, S.M. Malherbe
/ Integration
of systems
offirst order
ODE’s
LONG WRITE-UP 1. Introduction
as
This write-up documents a set of subroutines which apply Runge-Kutta (RK) or predictor-corrector (PC) methods for the numerical integration of systems of first order ordinary differential equations (ODE’s). It is not uncommon to find variable-order, variable step size routines either as separate packages [l] or as part of a mathematical software library such as NAGLIB [2]. However such general purpose routines tend to be lengthy and can require more storage than is absolutely necessary for specific applications. Furthermore, the changes in order or step size are often internal to the package and beyond the control of the user. For the specific application discussed in ref. [3] these considerations lead to the design of a simple ODE package which allows for the selection of fixed-order, fixed-step size integration formulae of either RK or PC type taken from the text of Lapidus and Seinfeld [4]. The subroutine package can be either applied as it stands or modified for specific applications. For example, other integration formulae may be substituted, or (in the PC case) adaptive step size algorithms tested. In the application discussed in ref. [3] the package proved useful in monitoring accuracy and stability in the integration of a system of nonlinear ODES. As part of the package systems of non-stiff and stiff ODE’s (systems I, II, III and VI of ref. [4]) are included as test cases as discussed in section 5. Some discrepancies are apparent between the results of the present package and those tabulated in ref. [4] even though the same integration formulae are used. The remainder of this section develops some detail of the general background to subsequent sections. Consider a system of N, first order ODE’s yII=fr(x,
Yl?...,
Xv,)
(1)
with initial condition Y,(xo) =Yo,,
i= l,...,
N,.
(2)
When written in vector form eqs. (l), (2) are given
y’=f(x,
Y>
(3)
and Yo =Y(xo>.
(4)
Define the sequence { x, } of equally spaced values of the independent variable x as x,=xg+nh,
(5)
where h is a constant step size and n is the integer sequence n = 0, 1, 2,. . . , N. A sequence of approximations { y,, } to { y( x,,)} can be generated by a linear multistep integration formula of the general type k
Y,+l=
c
k
(YiYn+i-r + h c
i=l
PiY,l+i-i+
T(x,
A),
1=0
(6) where the 2k + 1 parameters (~i,. . . , cxkr &,, . . . , fik are constants independent of h. If k = 1 then eq. (4) gives a single step method. If & = 0 then the formula is explicit and if &, # 0 it is implicit and requires a value for Y,‘+I = (dY/dx)],=,n+,.
(7)
Two important classes of methods for obtaining integration formulae such as eq. (6) are known as Runge-Kutta (RK) and predictor-corrector (PC). The RK method generates y’, f(x, y) at the points x, and also at values of x, y between successive x,, y, whereas the PC method generates Y’, f(x, Y) only at the points x,. Each method has associated with it an error T(x, h) which is the error produced at step (n + 1) if y, and y,’ on the right hand side of eq. (6) are known without error. In T(x, h) x is in the interval [X,-k, x,+1 1. Sources of error are truncation error and round-off error. Truncation error occurs because the rhs of eq. (6) is not an exact representation of Y, + i due to the use of a finite interpolating polynomial for y. Round-off error is due to finite precision arithmetic on the computer. The calculation of yn+, from inexact values for y,, y,’ propagates error through subsequent steps. This source
ofsystemsoffirst order
G. D&c, SM. Malherbe / Integration
of error is studied under the question of stability of a given integration formula. Detailed error and stability analyses for the integration formulae used here are to be found in ref. [4]. The second to sixth order formulae applied here used only explicit single step forms of eq. (6) which correspond to &, = 0 and ai = 0, i > 1. This choice restricts the PC formulae (as described in section 3) to the Adams-Bashforth predictors and the Adams-Moulton correctors. In the numerical experiments reported in ref. [4] these PC formulae were found to perform well. The corresponding selection of RK single step formulae is described in section 2. In section 2 and 3 the vector notation introduced above is suppressed to simplify the exposition. 2. Runge-Kutta
method
The RK method evaluates f(x, y) within the interval (x,, y,) to (x,+i, y,+r) in a general single step formula of the type =y,+
Yn+l
e
(8)
Wjrj,
j=1
where wi are weight coefficients and u is the number of evaluations of f(x, y) with j-l
q=hf
X,+cjh,
y,,+
C
Y,>.
00)
Different sets of parameters { wj}, {c,}, {a,,,,}, determine different RK integration formulae and each set specifies (x, y) values at which f(x, y) is evaluated. Each set of parameters is conveniently represented in the form of a Butcher Table [4] c
A
t
WT
Explicit RK formulae (which are the only type considered here) require only r,,,, m = 1,. . . , j - 1 to calculate 5. Thus the explicit form of eq. (11) would be 0 ’ c2
a21
c3 . . .
a31 . . .
a32
c,
au1
au2
. . .
auu--l
WI
w2
...
W”_l
(12) W”
RK formulae are classified by the pair (p, u) where p is the order of the RK formula. Second to sixth order RK formulae as taken from Lapidus and Seinfeld [4] are listed in tables 1 to 7. Tables 4 and 6 are alternative fourth and fifth order formulae which allow direct comparison with results tabulated in ref. [4]. The tables have been multiplied by a suitable constant B so that all the parameters wj, c,, a,, appear as integers. The specific formulae chosen are Improved Euler (second order), Nystrom (third order), Classic and Kutta (fourth order), Butcher (fifth and sixth order). For these formulae u(p) is the number of terms in the summation of eq. (8) for an order p RK formula and has the values u(P)=P,
2~~4, u(6) = 7,
u(5) = 6,
(13)
aj,,,r,,,
for j = 1, 2,. . . , u, where by definition c, = 0 and
=hfb,,
295
m=l
i
‘1
ODE’s
where A is a lower triangular matrix containing the multiplying factors aj,, c is the vector of coefficients c, and wT is the transposed vector containing the weights wj.
3. Predictor-corrector
method
In the PC method a predictor formula of the type in eq. (6) is used to obtain an approximate value j,+r for y,+r. A corrector formula of the same type then uses this value to obtain an improved estimate for y,+r. The basic steps in the method are:
Table 1 Runge-Kutta 02
second order (2, 2), B = 2
I
2
296
G. Delic, SM.
Table 2 Runge-Kutta
Malherbe
/ Infegrafion
1. Prediction 2. Evaluation
third order (3,3), B = 24
0 16 16
Table 3 Runge-Kutta 0 8 16 24
Table 4 Runge-Kutta
Y -,‘+, =fk+*,
16 0
16
6
9
3. Correction 4. Evaluation
9
fourth order (4,4), B = 24 8 -8 24
24 -24
24
3
9
9
3
fourth order (4,4), B = 24
Table 5 Runge-Kutta 0 1260 2520 5040 7560 10080 I
12 0 0
12 0
24
4
8
8
4
fifth order (5,6), B = 10080
1260 0 5040 1890 - 7200
2520 -10080 0 5760
10080 0 17280
5670 -17280
11520
784
0
3584
1344
3584
784
Table 6 Runge -. Kutta fifth order (5,6), B = 10080 0 2520 2520 5040 7560 10080
1260 -5040 0 2880
10080 0 17280
5670 -17280
11520
784
0
3584
1344
3584
Table 7 Runge-Kutta
484
3564
3564
-140%
-1408
of Y,+ i. of Y,‘, 1.
Yn+l
-rn+,.
(15)
784
26;_7680 0
(14)
6. The estimate of eq. (15) may be used in an adaptive PC algorithm where the step size is adjusted during the integration.
Sixth Order (6,7), B = 5280
~~;;jj~:~~~
Y,+,>.
1. The form of eq. (6) used for the corrector may be different from that used as the predictor and would therefore.have a different truncation error. 2. Convergence properties of the iteration determine the efficiency of the corrector step. is 3. The number of f(x, y) function evaluations determined by the number of iterations. 4. Starting values for y, and Y,’ must be provided with the number of such values being determined by the order of the PC formula. error en+i 5. Error estimates for the truncation are easily obtained from c n+1=
2520 1260 0 1890 - 4320
of a value JJ,,,~ for Yn+,. of the derivative from the ODE
The method is referred to as the PE(CE)s method by Lapidus and Seinfeld [4] since steps (3) and (4) are iterated S times. In a alternative method, namely PEC, step (4) is omitted but the function value using the predicted value j,,+i is used in the next cycle. Steps two and three may be iterated to give a P(EC)s method. Thus with S = 1 the PEC method gives Y,+i and jji+, while the PECE method gives Y,,,, and y,‘,,. Some features which are specific to the PC method are as follows,
0 12 24
of systems offirst order ODE’s
484
The PC formulae used here are of the Adams form with CX,= 1, 0~;= 0, i > 1 in eq. (6) and the predictor is the Bashforth form with & = 0 while the corrector is the Moulton form with & Z 0. The PC formulae together with their respective truncation errors are listed in tables 8 and 9. The tabulated values are the numbers dkP, where the divisors d, are given in the second column in each case. In the case of the predictor formulae the order p = k and the expression for the (local)
G. Delic, S.M. Malherbe / Integration of systems Table 8 Adams-Bashforth d,
d,&
d&z
d,&
2 3 4 5 6
2 12 24 720 1440
3 23 55 1901 4217
-1 -16 -59 - 2714 - 1923
5 31 2616 9982
Table 9 Adams-Moulton
correctors
k
d,
d.&
d,P,
d&
1 2 3 4 5 6
2 12 24 720 1440 60480
1 5 9 251 475 19087
1 8 19 646 1427 65112
-1 -5 - 264 - 798 -46461
T(x,
error assumes \
dk&
d&
-9 -1274 - 7298
d,&
1 106 482 37504
251 2887
dk&
-19 - 173 -20211
the form
h) = chp+lp+l)(s).
d&s
27 6312
d&
T(x, h)/(hk+lyCk+‘))
- 475
5/12 9/24 251/720 475/1440 19087/60480
dk&
T(x, b)/(hk+2y’k+2’)
-863
-l/12 -l/24 - 19/720 - 27/1440 - 863/60480 - 1375/120960
defined (16)
The term y (pfl) denotes the derivative of order p + 1 evaluated at some point { on the interval [X,-k, x,+1 ] and an expression similar to eq. (16) holds for the corrector formulae. Values of the constant C are listed in the last column of tables 8 and 9.
COEF
PRECOR
4. Description of the program The package consists of the driving routine MAIN, BLOCK DATA and subroutines PRECOR, RUNGEK, COEFF, ACTUAL, FCN, COMPR, FWDTRN, BCKTRN and OUTPUT. A detailed description of the routine functions is given in the source listing. A brief description follows. MAIN
297
predictors
k
truncation
offirst order ODE’s
Reads input data and initial values for the test case(s). It calls subroutine COEFF to calculate the coefficients for the RK (see section 2) and PC (see Section 3) methods. MAIN invokes either the PC or the RK method according to the user
RUNGEK
parameter
METHOD
= 1
(PC), 2 (RK). Assigns values from BLOCK DATA and calculates the coefficients for the RK method (see section 2, eqs. (12) (13) and tables 1 to 7) and the PC method (see section 3, tables 8 and 9). This implements the PC method described in section 3. It calls subroutine RUNGEK to obtain the required K-1 starting values (where p = K is the order of the method). There are two options selected by setting the parameter OPTION. For OPTION = 1 the corrector is of the same order as the predictor while with OPTION = 2 the corrector order is one less than the predictor. The corrector is iterated S times where S is set using the input parameter ITERNO. This routine implements the RK method described in section 2 for eqs. (8)-(10). It has two parameters, T the initial x value for the integration and TFINAL which is the x
G. D&c, S. M. Malherbe
298
ACTUAL
FCN COMPR
FWDTRN
BCKTRN
OUTPUT
/ Integration
value at which the integration terminates. This is a user defined routine which allows the calculated values from the integration to be compared to the actual values of the function at a given x value when they are known. Evaluates the function f(x, y) of eq. (3). Shifts the array passed to it down by one index allowing the calculated and function values to be stored as arrays of dimension N, by 6. For the stiff system test case this routine uses the forward transform matrix FWD to forward transform vector YIN (see section 5). For the stiff system test case this routine uses the back transform matrix BT to back transform vector A00 (see section 5). Prints out the calculated value, the actual value, the absolute (ERTYPE = 1) or the relative error (ERTYPE = 2) and the difference between the corrected and predicted values (for the PC method only), at step multiples of the parameter BLOCSZ. OUTPUT also prints out the total number of function evaluations.
of systems offirst order ODE’s
The system of eqs. (17), (18) (see also chapter 6 of ref. [4]) has the analytical solution _yi(x) = e-O.iX + eso”, y2(x)
= ee50X,
ys(x)
= eesoX + e-i20x,
and is a suitable model for equations which result from the spectral transform method of ref. 131. In this method the matrix A is written in diagonal form A = YAZT
ZTY = YZT = I.
(23)
and
(24)
The test cases consisted of three non-stiff systems and one stiff system. These examples were systems I, II, III and VI of table 2.3 in ref. [4]. The stiff case is system VI given as
Y
-0.1 0
A= i and Y;=pJJl.
I(Z) = 11y”’
0
- 49.9 -50 70
0 0 -120
I
(25)
with the forward transform defined as y”’ = zry
(26)
and the back transform defined as y = Yy(=).
where
(22)
For the matrix A of (18) a simple calculation shows that
From the method of ref. [3] eq. (17) becomes
(17)
(21)
where A is the diagonal matrix of eigenvalues h, = -0.1. h, = -50, A, = -120 and Y, Z have columns which are eigenvectors of A, AT respectively with normalisation
5. Test cases
Y1 = AY,
(20)
(18)
(27)
Then in place of eq. (17) the PC or RK method is applied to eq. (25) with the initial vector
yo”’ = zryo
(19)
(28)
and at any subsequent value of x the solution y is obtained from the back transform of eq. (27).
G. Delic, S.M. Malherbe / Integration of systems
Non-stiff systems initial values are
I, II,
III
and
corresponding
.Y’= -y,
y,=l,
(29)
y’=
y,=l,
(30)
y,=o,
(31)
+y,
y’=l-y, with respective
analytical
solutions
y(x)
= eeX,
(32)
y(x)
= e+X,
(33)
y(x)
= 1 - eex.
(34)
These non-stiff systems have been included because of differences in accuracy in comparing the present package with results given in ref. [4] despite the use of the same integration formulae. These comparisons implement the RK fourth and fifth order formulae of tables 4 and 6. As an example the RK fourth order formula in the present package reproduced the errors of table 2.8 in ref. [4] but with an increasing discrepancy with decreasing step size which reduced from to five to two significant figures for h in the range 0.25 to 0.01. The remainder of this section describes the four test run outputs. The package is set up to run for system VI (test run output 2 and 4) and needs only minor modifications to subroutines FCN and ACTUAL to run systems I, II, III (test run output 1 and 3). To compare RK and PC efficiency the number of calls to FCN is printed at the end of each run.
Test run output Test run I
This shows results for systems I, II, III corresponding to I = 1, 2, 3 for the RK Classic fourth order formula of table 4. Results of table 2.8 in ref. [4] are reproduced for h = 0.25 with system I. While not shown here, results of table 2.5 in ref. [4], were also reproduced for both systems I and II.
offirst order
ODE’s
299
Test run 2 This shows results for the stiff case of system VI for the fifth order RK formula of table 6. The three eigenmodes correspond to I = 1, 2, 3 respectively. Results tabulated in tables 2.9 of ref. [4] are reproduced for eigenmodes 2 and 3 but not for eigenmode 1. Calculations with all other step sizes listed in Table 2.9 of ref. [4] showed the same feature. Further calculations for h = 0.01 with the fourth order formula showed the results to be one to two orders of magnitude better than those listed in table 6.3 of ref. [4]. Test run 3
This shows the results for systems I, II, III corresponding to I = 1, 2, 3 for the fourth order PC method with S = 3. The last column E is the difference between corrected and predicted values in eq. (15) as computed in the last iteration. The relative error compares favourably with RK fourth order but costs fewer function calls. The two additional iterations improved the result for all three systems but more so for systems I and III. Further calculations for system III with this step size showed better performance than listed in table 4.6 of ref. [4] for p = 3 with S = 1. Furthermore the same case with S = 3 showed more than an order of magnitude improvement whereas in ref. [4] it does not. Test run 4
This shows results for the stiff case of system VI for the fifth order PC method with S = 3. The accuracy is comparable to test case 2 with a fourth order RK formula but with approximately half the number of function evaluations. Further calculations with a fourth order formula and h = 0.01 showed the solution at x = 0.4 to be unstable for all three eigenmodes with S = 1. However, a comparison with table 6.3 of ref. [4] showed stable solutions for S = 2 with similar or better accuracy for the first two eigenmodes and the third unstable while all three eigenmodes were stable at x = 0.4 with S = 6.
G. D&c, S.M. Malherbe / Integration
300
Acknowledgements This research was performed
under contract
the Plasma Physics Division of the Atomic Corporation of South Africa Limited.
to
Energy
J.W. Daniel, E. Denman and for Boundary-Value Problems
Ordinary Differential Equations (Lecture Notes in Computer Science, Nr. 76, Springer-Verlag, Berlin, 1979). [2] Numerical Algorithms Group Library Manual (Mark ll),
NumericalAlgorithmsGroup, Oxford, 1984. [3] G. Delic, J. Math. Phys. 28 (1987) 39.
.141_ L. LaDidus
and J.H. Seinfeld, Ordin& Differential Equations York, 1971).
References [l] B. Childs, M. Scott, Nelson, eds., Codes
of systems offirst order ODE’s
P. in
Numerical (Academic
Solution of Press, New
G. Delic, SM. Malherbe / Integration of system offirst order ODE’s
TEST RUN INPUT SYSTEMS 0.0 3
I,Il,lI 10.0 ;0
A::
l:o
0.0 1.0 0.0 0.0 0.0 1.0
:*: l:o
Test
I 4 2
input
SYSTEM
VI 0.005 5 2 : 5773:0269189625787 0:577350269189625787 0.577350269189625707 -1.0 1.73205080756887719 -1.0 -50.0 1.0
Test
0.5
run
SYSTEMS
input l.II.III 16.0.
0.0
4 2
i.0
0.0 0.0
~~00
1.0
Test
0: 0 1 .o 0.0 0.0 1.0
SYSTEM 0.0 3 1.0 0.0 0.0 1.0 0.0
0.0 -0.1 2.0
run
input
25
1
::: 1.0 :.: l:o -120 0 2,o
3
3 1.0
::: 0.0 1.0
2
2
0.0 1.0 3 0.0 0.0 1.0 0.0 0.0 10.1 2.0
8
A.: 0:o 1.0 0.0 0.0
g.0" l:o run
0.250 2 0.0 0.0
0.250 1 0.0 0.0 1.0 0.0 0.0 1.0
8
2
::: 4
VI 0.005 5 1 &5773;0269189625787 0.577350269109625787 0.577350269109625707 -1.0 1.73205080756807719 -1.0 0.5
-50.0 1.0
25
1 0.0 0.0 1.0 ::o" 1.0 -120.0 2.0
301
302
G. D&c, SM.
ofsystemsoffirst order
Malherbe / Integration
ODE’s
TEST RUN OUTPUT FIRST
ORDER
ODE
INTEGRATION USING
INTEGRATION
RANGE
(
0.0
.
10.0000)
STEP
SIZE
H =
CALCULATED 0.135346141957132571D+OO 0.7388665273572B5885D+Ol 0.864653858042867401D+OO
: 3
I
CALCULATED O.l83185781426802750D-01 0.545923745249014694D+O2 0.98168142185731970OD+OO
I 1 7 3
CALCULATED 0.247934887775202935~-02 0.4033647818540229650+03 0.997520651122247928D+OO
I
CALCULATED 0.335570305169483366D-03 0.298032735626711087D+04 0.999664429694830492D+OO
run
output
FIRST
ORDER
RANGE
(
0.0
,
0.5000)
STEP
SIZE
1 2 3
1 2 3
CALCULATED 0.975313638694454754D+OO 0.372666612242002666D-05 0.3726666216109402660-05
OF
H =
CALCULATED
2 3
I
CALCULATED 0.951229424514601477D+OO O.l38880403879931175D-10 O.l38880403879931261D-10
FUNCTION
= TEST
480 RUN
FOR
SYSTEM
ORDER
K =
BLOCSZ
5
=
25
0.1250 ABSOLUTE ERROR -.335421834751770120D-08 -.335421852651947988D-08 -.353910147704843830D-08
0.2500 ABSOLUTE ERROR -.129499883039230212D-10 -.129503412879631305D-10 -.129504544344124665D-10
0.3750 ABSOLUTE
-.369704267200177128D-13 -.375000924030652613D-13 -.3750009245517763990-13
0.9631944249149548550+00 0.719413303032561335D-08 0.719413303035423793D-08
ERROR
0.5000 ABSOLUTE ERROR 0.541233724504763813D-15 -.965230285046578347D-16 -.965230285046578347D-16
ACTUAL 0.951229424514602018D+00 O.l38879438649646128D-10 O.l38879438649646215D-10 NUMBER
:
RUNGE-KUTTA
ACTUAL
TOTAL
CALLS
PACKAGE
ACTUAL 0.975313638681504766D+OO 0.37266531720787387OD-05 0.372665326565496825D-05
x= 1 2 3
RELATIVE ERROR -.401242858251871595D-03 0.264434327820700567D-03 O.l82172246732899583D-07
ACTUAL 0.989508254630109177D+00 O.l93045413622772166D-02 O.l93076003854822348D-02
x= 0.963194424914991826D+OO 0.719417053041801641D-08 0.7194170530446693110-08
RELATIVE ERROR -.320981409001717987D-03 0.211553056889359650D-03 0.107713400826058441D-06
0.00500
x=
:
a
RELATIVE ERROR -.240726399068044142D-03 O.l58668988789710483D-03 0.598183835206303743D-06
INTEGRATION
x= CALCULATED 0.989508257984327524D+OD O.l93045749044624818D-02 O.l93076357764970053D-02
I
=
10.0000
USING INTEGRATION
BLOCSZ
8.0000
NUMBER ODE
I,II,III
6.0000
ACTUAL 0.453999297624848508D-04 0.220264657948067161D+05 0.999954600070237523D+OO
2
SYSTEMS
RELATIVE ERROR -.160477827934384288D-03 O.l05782123373425274D-03 0.299409266619492296D-05
ACTUAL 0.335462627902511799D-03 0.2980957987041728300+04 0.999664537372097498D+OO
CALCULATED 0.454181461600671844D-04 0.220206412411300007D+05 0.99995458185383991lD+OD
FOR
RELATIVE ERROR -.802356950840200623D-04 0,5289246i149274090470-04 O.l25583018589638564D-04
ACTUAL 0.247875217666635838D-02 0.403428793492735llOD+03 0.997521247823333651D+OO
TOTAL Test
K = 4
ACTUAL 0.183156388887341804D-01 0.545981500331442362D+02 0.981684361111265821D+OO
x= I
RUN
4.0000
x=
: 3
ORDER
ACTUAL O.l35335283236612688D+OO 0.738905609893065018D+01 0.864664716763387312D+OO
x=
: 3
TEST
2.0000
x=
: 3
:
RUNGE-KUTTA
0.25000
x= I
PACKAGE
OF
FUNCTION
CALLS
=
1818
VI
run
3
(
0.0 ,
10.0000)
CALCULATED 0.453459079295524943D-04 0.220429713008490489D+05 0.999954654092070419D+00
I 1
_.^^,
:
FIRST
.^_ ,..,
CALCULATED 0.3351519066869318090-03 0.2982695652975190400+04 0.9996648480933130330+00
I 1 2 3
.^^_,
CALCULATED 0.2477109968343606050-02 0.403596830792699052D+03 0.9975228900316563670+00
I 1 2 3
^_.^.
CALCULATED O.l83083362285561945D-01 0.5461180783343751700+02 0.9816916637714437760+00
: 3
I
0.864682965655215713D+00
CALCULATED 0.1353170343447842730+00 0.738967535732674730D+Ol
RANGE
output
1 2 3
I
INTEGRATION
Test
X=
2.0000
H = 0.2500000
x= 4.0000
ACTUAL 0.1353352832366126880+00 0.73890~6098930650180+01 0.864664716763387312D+00
SIZE
6.0000
8.0000
TOTAL
NUMBER
OF
BLOCSZ
RUN
=
FOR
=
390
RELATIVE ERROR O.l18991005525731378D-02 -.749348815016180896D-03 -.540242855947055639D-07
RELATIVE ERROR 0.926246889326484614D-03 -.5829219804558662060-03 -.3108254858694123830-06
RELATIVE ERROR 0.662514122311699621D-03 -.4165228230467581960-03 -.164628906532041695D-05
RELATIVE ERROR 0.3987117360387415100-03 -.250151338186178657D-03 -.743890854050894981D-05
8
SYSTEMS
RELATIVE ERROR O.l34842085463476737D-03 -.8380751043245300410-04 -.2110516536017628500-04
K = 4
: TEST
CALLS
ORDER
FUNCTION
10.0000
ACTUAL 0.4539992976248485000-04 0.220264657948067161D+05
x =
ACTUAL 0.335462627902511799D-03 0.2980957987041728300+04 0.999664537372097498D+OO
x=
ACTUAL 0.247875217666635838D-02 0.40342879349273511OD+O3 0.997521247823333651D+OO
x=
PACKAGE
PREDICTOR-CORRECTOR
INTEGRATION
USING
ODE
ACTUAL O.l83156388887341804D-01 0.5459815003314423620+02 0.9816843611112658210+00
STEP
ORDER
S =
3
OPTION
=
-.265890461:764337640-09 0.39450334594221203bD-01 0.265890462425311114D-09
-.196519816;72285863D-08 0.533813431479757128D-02 0.1965198170861270910-08
-.145247927:49049502D-07 0.7223177764277325600-03 O.l4524?927157976520D-07
E -.1073528346877405950-06 0.977388239711274309D-04 O.l0735283467?332254D-06
-.793070793:56220882D-06 0.1322498312905651120-04 0.7930707933700986700-06
ITERATIONS
1,11,111
2
run
4
CALCULATED 0.95122942451461D206D+OO 0,138967652607537578D-10 O.l42639648789881117D-10
: 3
I 1 2 3
I
CALCULATED 0.963194424918334111D+DO 0.719751281556389742D-08 0.715908123567748933D-08
0.5000)
I
,
CALCULATED 0.975313639816041064D+00 0.37277877087345061OD-05 0.373181012025953838D-05
0.0
I 1 2 3
(
: 3
RANGE
output.
CALCULATED 0.989508523169545434D+OO O.l93072267566415565D-02 O.l93063564496288317D-02
INTEGRATION
rest FIRST
SIZE 0.1250
H = 0.0050000 x=
0.2500
0.3750
TOTAL
NUMBER
OF
ORDER
FUNCTION
0.5000
ACTUAL 0.951229424514602018D+OO O.l38879438649646128D-10 O.l38879438649646215D-10
x=
ACTUAL 0.963194424914954855D+OO 0.719413303032561335D-08 0.719413303035423793D-08
x=
ACTUAL 0.975313638681504766D+OO 0.372665317207873870D-05 0.372665326565496825D-05
x=
PACKAGE
5
TEST
BLOCSZ
RUN
=
FOR
25
SYSTEM
=
960
ABSOLUTE ERROR -.818789480661052949D-14 -.882139578914496007D-14 -.376021014023490212D-12
ABSOLUTE ERROR -.337925520899062803D-11 -.337978523828407489D-11 0.350517946767486024D-10
ABSOLUTE ERROR -.113453629824178392D-08 -.113453665576740294D-08 -.515685460457013427D-08
ABSOLUTE ERROR -.268539436257264619D-06 -.2685394364339895730-06 O.l24393585340318302D-06
K =
:
CALLS
PREDICTOR-CORRECTOR
INTEGRATION
USING
ODE
ACTUAL 0.989508254630109177D+OO D.l93045413622772166D-02 O.l93076003854822348D-02
STEP
ORDER
VI
S =
3
OPTION
E O.l8470606144154825lD-16 O.l84706061441548251D-16 -.488409904565982497D-12
0.956642944:09092232D-14 0.956642944809092232D-14 0.51129102523574545OD-10
0.495471408:55925636D-11 0.495471408355925636D-11 -.534528062523821371D-08
0.256618117:66739735D-08 0.256618117066739735D-08 0.563836206335483095D-06
ITERATIONS
=
2
48 (1988) 293-304
293
SUBROUTINES FOR INTEGRATION OF SYSTEMS OF FIRST ORDER ODE’s G. DELIC
and S.M. MALHERBE
Centre for Nonlinear Studies and Department Johannesburg, South Africa Received
of Physics, University of the Witwatersrand,
WITS 2050,
1 July 1987
PROGRAM SUMMARY Title of program:
RKPC
Peripherals used: card reader,
Catalogue number: ABBH
No. of cards in combined program and test deck: 617
Program obtainable from: CPC Program Library, versity of Belfast, N. Ireland (see application issue) Computer: IBM 3083524; Installation: Witwatersrand Computer Centre Operating system: VM/SP Programming language LANGLVL (77)
line printer
HP0
Queen’s Uniform in this
University
of
the
Nature of the problem Numerical integration of the vector equation Method of solution Application of second or predictor-corrector
4.2
used: FORTRAN
Keywords: Runge-Kutta, predictor-corrector, ordinary differential equations
VS LEVEL
1.4.1,
High speed storage required: 20272 bytes for OPT(O) compile No. of bits in a byte: 8
OOlO-4655/88/$03.50 0 Elsevier Science Publishers (North-Holland Physics Publishing Division)
first
dy/dx
order
= f( x, y).
to sixth order, fixed step, Runge-Kutta integration formulae.
Restrictions on the complexity of the problem Vector dimension set for 20 components. Typical running times Dependent on required range and accuracy: is in the range 0.14 to 0.18 s.
B.V.
for the test cases it
294
G. D&c, S.M. Malherbe
/ Integration
of systems
offirst order
ODE’s
LONG WRITE-UP 1. Introduction
as
This write-up documents a set of subroutines which apply Runge-Kutta (RK) or predictor-corrector (PC) methods for the numerical integration of systems of first order ordinary differential equations (ODE’s). It is not uncommon to find variable-order, variable step size routines either as separate packages [l] or as part of a mathematical software library such as NAGLIB [2]. However such general purpose routines tend to be lengthy and can require more storage than is absolutely necessary for specific applications. Furthermore, the changes in order or step size are often internal to the package and beyond the control of the user. For the specific application discussed in ref. [3] these considerations lead to the design of a simple ODE package which allows for the selection of fixed-order, fixed-step size integration formulae of either RK or PC type taken from the text of Lapidus and Seinfeld [4]. The subroutine package can be either applied as it stands or modified for specific applications. For example, other integration formulae may be substituted, or (in the PC case) adaptive step size algorithms tested. In the application discussed in ref. [3] the package proved useful in monitoring accuracy and stability in the integration of a system of nonlinear ODES. As part of the package systems of non-stiff and stiff ODE’s (systems I, II, III and VI of ref. [4]) are included as test cases as discussed in section 5. Some discrepancies are apparent between the results of the present package and those tabulated in ref. [4] even though the same integration formulae are used. The remainder of this section develops some detail of the general background to subsequent sections. Consider a system of N, first order ODE’s yII=fr(x,
Yl?...,
Xv,)
(1)
with initial condition Y,(xo) =Yo,,
i= l,...,
N,.
(2)
When written in vector form eqs. (l), (2) are given
y’=f(x,
Y>
(3)
and Yo =Y(xo>.
(4)
Define the sequence { x, } of equally spaced values of the independent variable x as x,=xg+nh,
(5)
where h is a constant step size and n is the integer sequence n = 0, 1, 2,. . . , N. A sequence of approximations { y,, } to { y( x,,)} can be generated by a linear multistep integration formula of the general type k
Y,+l=
c
k
(YiYn+i-r + h c
i=l
PiY,l+i-i+
T(x,
A),
1=0
(6) where the 2k + 1 parameters (~i,. . . , cxkr &,, . . . , fik are constants independent of h. If k = 1 then eq. (4) gives a single step method. If & = 0 then the formula is explicit and if &, # 0 it is implicit and requires a value for Y,‘+I = (dY/dx)],=,n+,.
(7)
Two important classes of methods for obtaining integration formulae such as eq. (6) are known as Runge-Kutta (RK) and predictor-corrector (PC). The RK method generates y’, f(x, y) at the points x, and also at values of x, y between successive x,, y, whereas the PC method generates Y’, f(x, Y) only at the points x,. Each method has associated with it an error T(x, h) which is the error produced at step (n + 1) if y, and y,’ on the right hand side of eq. (6) are known without error. In T(x, h) x is in the interval [X,-k, x,+1 1. Sources of error are truncation error and round-off error. Truncation error occurs because the rhs of eq. (6) is not an exact representation of Y, + i due to the use of a finite interpolating polynomial for y. Round-off error is due to finite precision arithmetic on the computer. The calculation of yn+, from inexact values for y,, y,’ propagates error through subsequent steps. This source
ofsystemsoffirst order
G. D&c, SM. Malherbe / Integration
of error is studied under the question of stability of a given integration formula. Detailed error and stability analyses for the integration formulae used here are to be found in ref. [4]. The second to sixth order formulae applied here used only explicit single step forms of eq. (6) which correspond to &, = 0 and ai = 0, i > 1. This choice restricts the PC formulae (as described in section 3) to the Adams-Bashforth predictors and the Adams-Moulton correctors. In the numerical experiments reported in ref. [4] these PC formulae were found to perform well. The corresponding selection of RK single step formulae is described in section 2. In section 2 and 3 the vector notation introduced above is suppressed to simplify the exposition. 2. Runge-Kutta
method
The RK method evaluates f(x, y) within the interval (x,, y,) to (x,+i, y,+r) in a general single step formula of the type =y,+
Yn+l
e
(8)
Wjrj,
j=1
where wi are weight coefficients and u is the number of evaluations of f(x, y) with j-l
q=hf
X,+cjh,
y,,+
C
Y,>.
00)
Different sets of parameters { wj}, {c,}, {a,,,,}, determine different RK integration formulae and each set specifies (x, y) values at which f(x, y) is evaluated. Each set of parameters is conveniently represented in the form of a Butcher Table [4] c
A
t
WT
Explicit RK formulae (which are the only type considered here) require only r,,,, m = 1,. . . , j - 1 to calculate 5. Thus the explicit form of eq. (11) would be 0 ’ c2
a21
c3 . . .
a31 . . .
a32
c,
au1
au2
. . .
auu--l
WI
w2
...
W”_l
(12) W”
RK formulae are classified by the pair (p, u) where p is the order of the RK formula. Second to sixth order RK formulae as taken from Lapidus and Seinfeld [4] are listed in tables 1 to 7. Tables 4 and 6 are alternative fourth and fifth order formulae which allow direct comparison with results tabulated in ref. [4]. The tables have been multiplied by a suitable constant B so that all the parameters wj, c,, a,, appear as integers. The specific formulae chosen are Improved Euler (second order), Nystrom (third order), Classic and Kutta (fourth order), Butcher (fifth and sixth order). For these formulae u(p) is the number of terms in the summation of eq. (8) for an order p RK formula and has the values u(P)=P,
2~~4, u(6) = 7,
u(5) = 6,
(13)
aj,,,r,,,
for j = 1, 2,. . . , u, where by definition c, = 0 and
=hfb,,
295
m=l
i
‘1
ODE’s
where A is a lower triangular matrix containing the multiplying factors aj,, c is the vector of coefficients c, and wT is the transposed vector containing the weights wj.
3. Predictor-corrector
method
In the PC method a predictor formula of the type in eq. (6) is used to obtain an approximate value j,+r for y,+r. A corrector formula of the same type then uses this value to obtain an improved estimate for y,+r. The basic steps in the method are:
Table 1 Runge-Kutta 02
second order (2, 2), B = 2
I
2
296
G. Delic, SM.
Table 2 Runge-Kutta
Malherbe
/ Infegrafion
1. Prediction 2. Evaluation
third order (3,3), B = 24
0 16 16
Table 3 Runge-Kutta 0 8 16 24
Table 4 Runge-Kutta
Y -,‘+, =fk+*,
16 0
16
6
9
3. Correction 4. Evaluation
9
fourth order (4,4), B = 24 8 -8 24
24 -24
24
3
9
9
3
fourth order (4,4), B = 24
Table 5 Runge-Kutta 0 1260 2520 5040 7560 10080 I
12 0 0
12 0
24
4
8
8
4
fifth order (5,6), B = 10080
1260 0 5040 1890 - 7200
2520 -10080 0 5760
10080 0 17280
5670 -17280
11520
784
0
3584
1344
3584
784
Table 6 Runge -. Kutta fifth order (5,6), B = 10080 0 2520 2520 5040 7560 10080
1260 -5040 0 2880
10080 0 17280
5670 -17280
11520
784
0
3584
1344
3584
Table 7 Runge-Kutta
484
3564
3564
-140%
-1408
of Y,+ i. of Y,‘, 1.
Yn+l
-rn+,.
(15)
784
26;_7680 0
(14)
6. The estimate of eq. (15) may be used in an adaptive PC algorithm where the step size is adjusted during the integration.
Sixth Order (6,7), B = 5280
~~;;jj~:~~~
Y,+,>.
1. The form of eq. (6) used for the corrector may be different from that used as the predictor and would therefore.have a different truncation error. 2. Convergence properties of the iteration determine the efficiency of the corrector step. is 3. The number of f(x, y) function evaluations determined by the number of iterations. 4. Starting values for y, and Y,’ must be provided with the number of such values being determined by the order of the PC formula. error en+i 5. Error estimates for the truncation are easily obtained from c n+1=
2520 1260 0 1890 - 4320
of a value JJ,,,~ for Yn+,. of the derivative from the ODE
The method is referred to as the PE(CE)s method by Lapidus and Seinfeld [4] since steps (3) and (4) are iterated S times. In a alternative method, namely PEC, step (4) is omitted but the function value using the predicted value j,,+i is used in the next cycle. Steps two and three may be iterated to give a P(EC)s method. Thus with S = 1 the PEC method gives Y,+i and jji+, while the PECE method gives Y,,,, and y,‘,,. Some features which are specific to the PC method are as follows,
0 12 24
of systems offirst order ODE’s
484
The PC formulae used here are of the Adams form with CX,= 1, 0~;= 0, i > 1 in eq. (6) and the predictor is the Bashforth form with & = 0 while the corrector is the Moulton form with & Z 0. The PC formulae together with their respective truncation errors are listed in tables 8 and 9. The tabulated values are the numbers dkP, where the divisors d, are given in the second column in each case. In the case of the predictor formulae the order p = k and the expression for the (local)
G. Delic, S.M. Malherbe / Integration of systems Table 8 Adams-Bashforth d,
d,&
d&z
d,&
2 3 4 5 6
2 12 24 720 1440
3 23 55 1901 4217
-1 -16 -59 - 2714 - 1923
5 31 2616 9982
Table 9 Adams-Moulton
correctors
k
d,
d.&
d,P,
d&
1 2 3 4 5 6
2 12 24 720 1440 60480
1 5 9 251 475 19087
1 8 19 646 1427 65112
-1 -5 - 264 - 798 -46461
T(x,
error assumes \
dk&
d&
-9 -1274 - 7298
d,&
1 106 482 37504
251 2887
dk&
-19 - 173 -20211
the form
h) = chp+lp+l)(s).
d&s
27 6312
d&
T(x, h)/(hk+lyCk+‘))
- 475
5/12 9/24 251/720 475/1440 19087/60480
dk&
T(x, b)/(hk+2y’k+2’)
-863
-l/12 -l/24 - 19/720 - 27/1440 - 863/60480 - 1375/120960
defined (16)
The term y (pfl) denotes the derivative of order p + 1 evaluated at some point { on the interval [X,-k, x,+1 ] and an expression similar to eq. (16) holds for the corrector formulae. Values of the constant C are listed in the last column of tables 8 and 9.
COEF
PRECOR
4. Description of the program The package consists of the driving routine MAIN, BLOCK DATA and subroutines PRECOR, RUNGEK, COEFF, ACTUAL, FCN, COMPR, FWDTRN, BCKTRN and OUTPUT. A detailed description of the routine functions is given in the source listing. A brief description follows. MAIN
297
predictors
k
truncation
offirst order ODE’s
Reads input data and initial values for the test case(s). It calls subroutine COEFF to calculate the coefficients for the RK (see section 2) and PC (see Section 3) methods. MAIN invokes either the PC or the RK method according to the user
RUNGEK
parameter
METHOD
= 1
(PC), 2 (RK). Assigns values from BLOCK DATA and calculates the coefficients for the RK method (see section 2, eqs. (12) (13) and tables 1 to 7) and the PC method (see section 3, tables 8 and 9). This implements the PC method described in section 3. It calls subroutine RUNGEK to obtain the required K-1 starting values (where p = K is the order of the method). There are two options selected by setting the parameter OPTION. For OPTION = 1 the corrector is of the same order as the predictor while with OPTION = 2 the corrector order is one less than the predictor. The corrector is iterated S times where S is set using the input parameter ITERNO. This routine implements the RK method described in section 2 for eqs. (8)-(10). It has two parameters, T the initial x value for the integration and TFINAL which is the x
G. D&c, S. M. Malherbe
298
ACTUAL
FCN COMPR
FWDTRN
BCKTRN
OUTPUT
/ Integration
value at which the integration terminates. This is a user defined routine which allows the calculated values from the integration to be compared to the actual values of the function at a given x value when they are known. Evaluates the function f(x, y) of eq. (3). Shifts the array passed to it down by one index allowing the calculated and function values to be stored as arrays of dimension N, by 6. For the stiff system test case this routine uses the forward transform matrix FWD to forward transform vector YIN (see section 5). For the stiff system test case this routine uses the back transform matrix BT to back transform vector A00 (see section 5). Prints out the calculated value, the actual value, the absolute (ERTYPE = 1) or the relative error (ERTYPE = 2) and the difference between the corrected and predicted values (for the PC method only), at step multiples of the parameter BLOCSZ. OUTPUT also prints out the total number of function evaluations.
of systems offirst order ODE’s
The system of eqs. (17), (18) (see also chapter 6 of ref. [4]) has the analytical solution _yi(x) = e-O.iX + eso”, y2(x)
= ee50X,
ys(x)
= eesoX + e-i20x,
and is a suitable model for equations which result from the spectral transform method of ref. 131. In this method the matrix A is written in diagonal form A = YAZT
ZTY = YZT = I.
(23)
and
(24)
The test cases consisted of three non-stiff systems and one stiff system. These examples were systems I, II, III and VI of table 2.3 in ref. [4]. The stiff case is system VI given as
Y
-0.1 0
A= i and Y;=pJJl.
I(Z) = 11y”’
0
- 49.9 -50 70
0 0 -120
I
(25)
with the forward transform defined as y”’ = zry
(26)
and the back transform defined as y = Yy(=).
where
(22)
For the matrix A of (18) a simple calculation shows that
From the method of ref. [3] eq. (17) becomes
(17)
(21)
where A is the diagonal matrix of eigenvalues h, = -0.1. h, = -50, A, = -120 and Y, Z have columns which are eigenvectors of A, AT respectively with normalisation
5. Test cases
Y1 = AY,
(20)
(18)
(27)
Then in place of eq. (17) the PC or RK method is applied to eq. (25) with the initial vector
yo”’ = zryo
(19)
(28)
and at any subsequent value of x the solution y is obtained from the back transform of eq. (27).
G. Delic, S.M. Malherbe / Integration of systems
Non-stiff systems initial values are
I, II,
III
and
corresponding
.Y’= -y,
y,=l,
(29)
y’=
y,=l,
(30)
y,=o,
(31)
+y,
y’=l-y, with respective
analytical
solutions
y(x)
= eeX,
(32)
y(x)
= e+X,
(33)
y(x)
= 1 - eex.
(34)
These non-stiff systems have been included because of differences in accuracy in comparing the present package with results given in ref. [4] despite the use of the same integration formulae. These comparisons implement the RK fourth and fifth order formulae of tables 4 and 6. As an example the RK fourth order formula in the present package reproduced the errors of table 2.8 in ref. [4] but with an increasing discrepancy with decreasing step size which reduced from to five to two significant figures for h in the range 0.25 to 0.01. The remainder of this section describes the four test run outputs. The package is set up to run for system VI (test run output 2 and 4) and needs only minor modifications to subroutines FCN and ACTUAL to run systems I, II, III (test run output 1 and 3). To compare RK and PC efficiency the number of calls to FCN is printed at the end of each run.
Test run output Test run I
This shows results for systems I, II, III corresponding to I = 1, 2, 3 for the RK Classic fourth order formula of table 4. Results of table 2.8 in ref. [4] are reproduced for h = 0.25 with system I. While not shown here, results of table 2.5 in ref. [4], were also reproduced for both systems I and II.
offirst order
ODE’s
299
Test run 2 This shows results for the stiff case of system VI for the fifth order RK formula of table 6. The three eigenmodes correspond to I = 1, 2, 3 respectively. Results tabulated in tables 2.9 of ref. [4] are reproduced for eigenmodes 2 and 3 but not for eigenmode 1. Calculations with all other step sizes listed in Table 2.9 of ref. [4] showed the same feature. Further calculations for h = 0.01 with the fourth order formula showed the results to be one to two orders of magnitude better than those listed in table 6.3 of ref. [4]. Test run 3
This shows the results for systems I, II, III corresponding to I = 1, 2, 3 for the fourth order PC method with S = 3. The last column E is the difference between corrected and predicted values in eq. (15) as computed in the last iteration. The relative error compares favourably with RK fourth order but costs fewer function calls. The two additional iterations improved the result for all three systems but more so for systems I and III. Further calculations for system III with this step size showed better performance than listed in table 4.6 of ref. [4] for p = 3 with S = 1. Furthermore the same case with S = 3 showed more than an order of magnitude improvement whereas in ref. [4] it does not. Test run 4
This shows results for the stiff case of system VI for the fifth order PC method with S = 3. The accuracy is comparable to test case 2 with a fourth order RK formula but with approximately half the number of function evaluations. Further calculations with a fourth order formula and h = 0.01 showed the solution at x = 0.4 to be unstable for all three eigenmodes with S = 1. However, a comparison with table 6.3 of ref. [4] showed stable solutions for S = 2 with similar or better accuracy for the first two eigenmodes and the third unstable while all three eigenmodes were stable at x = 0.4 with S = 6.
G. D&c, S.M. Malherbe / Integration
300
Acknowledgements This research was performed
under contract
the Plasma Physics Division of the Atomic Corporation of South Africa Limited.
to
Energy
J.W. Daniel, E. Denman and for Boundary-Value Problems
Ordinary Differential Equations (Lecture Notes in Computer Science, Nr. 76, Springer-Verlag, Berlin, 1979). [2] Numerical Algorithms Group Library Manual (Mark ll),
NumericalAlgorithmsGroup, Oxford, 1984. [3] G. Delic, J. Math. Phys. 28 (1987) 39.
.141_ L. LaDidus
and J.H. Seinfeld, Ordin& Differential Equations York, 1971).
References [l] B. Childs, M. Scott, Nelson, eds., Codes
of systems offirst order ODE’s
P. in
Numerical (Academic
Solution of Press, New
G. Delic, SM. Malherbe / Integration of system offirst order ODE’s
TEST RUN INPUT SYSTEMS 0.0 3
I,Il,lI 10.0 ;0
A::
l:o
0.0 1.0 0.0 0.0 0.0 1.0
:*: l:o
Test
I 4 2
input
SYSTEM
VI 0.005 5 2 : 5773:0269189625787 0:577350269189625787 0.577350269189625707 -1.0 1.73205080756887719 -1.0 -50.0 1.0
Test
0.5
run
SYSTEMS
input l.II.III 16.0.
0.0
4 2
i.0
0.0 0.0
~~00
1.0
Test
0: 0 1 .o 0.0 0.0 1.0
SYSTEM 0.0 3 1.0 0.0 0.0 1.0 0.0
0.0 -0.1 2.0
run
input
25
1
::: 1.0 :.: l:o -120 0 2,o
3
3 1.0
::: 0.0 1.0
2
2
0.0 1.0 3 0.0 0.0 1.0 0.0 0.0 10.1 2.0
8
A.: 0:o 1.0 0.0 0.0
g.0" l:o run
0.250 2 0.0 0.0
0.250 1 0.0 0.0 1.0 0.0 0.0 1.0
8
2
::: 4
VI 0.005 5 1 &5773;0269189625787 0.577350269109625787 0.577350269109625707 -1.0 1.73205080756807719 -1.0 0.5
-50.0 1.0
25
1 0.0 0.0 1.0 ::o" 1.0 -120.0 2.0
301
302
G. D&c, SM.
ofsystemsoffirst order
Malherbe / Integration
ODE’s
TEST RUN OUTPUT FIRST
ORDER
ODE
INTEGRATION USING
INTEGRATION
RANGE
(
0.0
.
10.0000)
STEP
SIZE
H =
CALCULATED 0.135346141957132571D+OO 0.7388665273572B5885D+Ol 0.864653858042867401D+OO
: 3
I
CALCULATED O.l83185781426802750D-01 0.545923745249014694D+O2 0.98168142185731970OD+OO
I 1 7 3
CALCULATED 0.247934887775202935~-02 0.4033647818540229650+03 0.997520651122247928D+OO
I
CALCULATED 0.335570305169483366D-03 0.298032735626711087D+04 0.999664429694830492D+OO
run
output
FIRST
ORDER
RANGE
(
0.0
,
0.5000)
STEP
SIZE
1 2 3
1 2 3
CALCULATED 0.975313638694454754D+OO 0.372666612242002666D-05 0.3726666216109402660-05
OF
H =
CALCULATED
2 3
I
CALCULATED 0.951229424514601477D+OO O.l38880403879931175D-10 O.l38880403879931261D-10
FUNCTION
= TEST
480 RUN
FOR
SYSTEM
ORDER
K =
BLOCSZ
5
=
25
0.1250 ABSOLUTE ERROR -.335421834751770120D-08 -.335421852651947988D-08 -.353910147704843830D-08
0.2500 ABSOLUTE ERROR -.129499883039230212D-10 -.129503412879631305D-10 -.129504544344124665D-10
0.3750 ABSOLUTE
-.369704267200177128D-13 -.375000924030652613D-13 -.3750009245517763990-13
0.9631944249149548550+00 0.719413303032561335D-08 0.719413303035423793D-08
ERROR
0.5000 ABSOLUTE ERROR 0.541233724504763813D-15 -.965230285046578347D-16 -.965230285046578347D-16
ACTUAL 0.951229424514602018D+00 O.l38879438649646128D-10 O.l38879438649646215D-10 NUMBER
:
RUNGE-KUTTA
ACTUAL
TOTAL
CALLS
PACKAGE
ACTUAL 0.975313638681504766D+OO 0.37266531720787387OD-05 0.372665326565496825D-05
x= 1 2 3
RELATIVE ERROR -.401242858251871595D-03 0.264434327820700567D-03 O.l82172246732899583D-07
ACTUAL 0.989508254630109177D+00 O.l93045413622772166D-02 O.l93076003854822348D-02
x= 0.963194424914991826D+OO 0.719417053041801641D-08 0.7194170530446693110-08
RELATIVE ERROR -.320981409001717987D-03 0.211553056889359650D-03 0.107713400826058441D-06
0.00500
x=
:
a
RELATIVE ERROR -.240726399068044142D-03 O.l58668988789710483D-03 0.598183835206303743D-06
INTEGRATION
x= CALCULATED 0.989508257984327524D+OD O.l93045749044624818D-02 O.l93076357764970053D-02
I
=
10.0000
USING INTEGRATION
BLOCSZ
8.0000
NUMBER ODE
I,II,III
6.0000
ACTUAL 0.453999297624848508D-04 0.220264657948067161D+05 0.999954600070237523D+OO
2
SYSTEMS
RELATIVE ERROR -.160477827934384288D-03 O.l05782123373425274D-03 0.299409266619492296D-05
ACTUAL 0.335462627902511799D-03 0.2980957987041728300+04 0.999664537372097498D+OO
CALCULATED 0.454181461600671844D-04 0.220206412411300007D+05 0.99995458185383991lD+OD
FOR
RELATIVE ERROR -.802356950840200623D-04 0,5289246i149274090470-04 O.l25583018589638564D-04
ACTUAL 0.247875217666635838D-02 0.403428793492735llOD+03 0.997521247823333651D+OO
TOTAL Test
K = 4
ACTUAL 0.183156388887341804D-01 0.545981500331442362D+02 0.981684361111265821D+OO
x= I
RUN
4.0000
x=
: 3
ORDER
ACTUAL O.l35335283236612688D+OO 0.738905609893065018D+01 0.864664716763387312D+OO
x=
: 3
TEST
2.0000
x=
: 3
:
RUNGE-KUTTA
0.25000
x= I
PACKAGE
OF
FUNCTION
CALLS
=
1818
VI
run
3
(
0.0 ,
10.0000)
CALCULATED 0.453459079295524943D-04 0.220429713008490489D+05 0.999954654092070419D+00
I 1
_.^^,
:
FIRST
.^_ ,..,
CALCULATED 0.3351519066869318090-03 0.2982695652975190400+04 0.9996648480933130330+00
I 1 2 3
.^^_,
CALCULATED 0.2477109968343606050-02 0.403596830792699052D+03 0.9975228900316563670+00
I 1 2 3
^_.^.
CALCULATED O.l83083362285561945D-01 0.5461180783343751700+02 0.9816916637714437760+00
: 3
I
0.864682965655215713D+00
CALCULATED 0.1353170343447842730+00 0.738967535732674730D+Ol
RANGE
output
1 2 3
I
INTEGRATION
Test
X=
2.0000
H = 0.2500000
x= 4.0000
ACTUAL 0.1353352832366126880+00 0.73890~6098930650180+01 0.864664716763387312D+00
SIZE
6.0000
8.0000
TOTAL
NUMBER
OF
BLOCSZ
RUN
=
FOR
=
390
RELATIVE ERROR O.l18991005525731378D-02 -.749348815016180896D-03 -.540242855947055639D-07
RELATIVE ERROR 0.926246889326484614D-03 -.5829219804558662060-03 -.3108254858694123830-06
RELATIVE ERROR 0.662514122311699621D-03 -.4165228230467581960-03 -.164628906532041695D-05
RELATIVE ERROR 0.3987117360387415100-03 -.250151338186178657D-03 -.743890854050894981D-05
8
SYSTEMS
RELATIVE ERROR O.l34842085463476737D-03 -.8380751043245300410-04 -.2110516536017628500-04
K = 4
: TEST
CALLS
ORDER
FUNCTION
10.0000
ACTUAL 0.4539992976248485000-04 0.220264657948067161D+05
x =
ACTUAL 0.335462627902511799D-03 0.2980957987041728300+04 0.999664537372097498D+OO
x=
ACTUAL 0.247875217666635838D-02 0.40342879349273511OD+O3 0.997521247823333651D+OO
x=
PACKAGE
PREDICTOR-CORRECTOR
INTEGRATION
USING
ODE
ACTUAL O.l83156388887341804D-01 0.5459815003314423620+02 0.9816843611112658210+00
STEP
ORDER
S =
3
OPTION
=
-.265890461:764337640-09 0.39450334594221203bD-01 0.265890462425311114D-09
-.196519816;72285863D-08 0.533813431479757128D-02 0.1965198170861270910-08
-.145247927:49049502D-07 0.7223177764277325600-03 O.l4524?927157976520D-07
E -.1073528346877405950-06 0.977388239711274309D-04 O.l0735283467?332254D-06
-.793070793:56220882D-06 0.1322498312905651120-04 0.7930707933700986700-06
ITERATIONS
1,11,111
2
run
4
CALCULATED 0.95122942451461D206D+OO 0,138967652607537578D-10 O.l42639648789881117D-10
: 3
I 1 2 3
I
CALCULATED 0.963194424918334111D+DO 0.719751281556389742D-08 0.715908123567748933D-08
0.5000)
I
,
CALCULATED 0.975313639816041064D+00 0.37277877087345061OD-05 0.373181012025953838D-05
0.0
I 1 2 3
(
: 3
RANGE
output.
CALCULATED 0.989508523169545434D+OO O.l93072267566415565D-02 O.l93063564496288317D-02
INTEGRATION
rest FIRST
SIZE 0.1250
H = 0.0050000 x=
0.2500
0.3750
TOTAL
NUMBER
OF
ORDER
FUNCTION
0.5000
ACTUAL 0.951229424514602018D+OO O.l38879438649646128D-10 O.l38879438649646215D-10
x=
ACTUAL 0.963194424914954855D+OO 0.719413303032561335D-08 0.719413303035423793D-08
x=
ACTUAL 0.975313638681504766D+OO 0.372665317207873870D-05 0.372665326565496825D-05
x=
PACKAGE
5
TEST
BLOCSZ
RUN
=
FOR
25
SYSTEM
=
960
ABSOLUTE ERROR -.818789480661052949D-14 -.882139578914496007D-14 -.376021014023490212D-12
ABSOLUTE ERROR -.337925520899062803D-11 -.337978523828407489D-11 0.350517946767486024D-10
ABSOLUTE ERROR -.113453629824178392D-08 -.113453665576740294D-08 -.515685460457013427D-08
ABSOLUTE ERROR -.268539436257264619D-06 -.2685394364339895730-06 O.l24393585340318302D-06
K =
:
CALLS
PREDICTOR-CORRECTOR
INTEGRATION
USING
ODE
ACTUAL 0.989508254630109177D+OO D.l93045413622772166D-02 O.l93076003854822348D-02
STEP
ORDER
VI
S =
3
OPTION
E O.l8470606144154825lD-16 O.l84706061441548251D-16 -.488409904565982497D-12
0.956642944:09092232D-14 0.956642944809092232D-14 0.51129102523574545OD-10
0.495471408:55925636D-11 0.495471408355925636D-11 -.534528062523821371D-08
0.256618117:66739735D-08 0.256618117066739735D-08 0.563836206335483095D-06
ITERATIONS
=
2