EXPFIT4 — a FORTRAN program for the numerical solution of systems of nonlinear second-order initial-value problems

ELSEVIER

Computer Physics CommunicmicJJs i00 (1997) 71-80

...........

EXPF1T4 - a FORTRAN program for the numerical solution of systems of nonlinear second-order initial-value problems L.Gr. Ixaru ~, H. De Meyer 2, G. Vanden Berghe 3, M. Van Daele Department of Applied Mathematics and Computer Science. Universit¢it Gent. grijgslaan 281-$9, B-9000 Gtnfo I~lgi~m Received 8 January 1996

Abstract in this paper we present a FOI~rRAN p~gram which solves the initiai-valne problem associaeed with nons+Jff systems of the form y " = f ( x , y ) . The program is based on a family of exponential-fitted four-step met.hods, The ~ c~ particular|y suited to solve second-order initial-value problems arising from physics. PACS: 02.30.Hq; 03.20.+i, 95.10.Ce Keywords: Newtonian physics; Celestial mechanics; Linear and nonlinear second-order differential equations; InRial-valt~¢~ m s ; Exponenthi fittit~g; Multi~ep methods; Stepsize control

PROGRAM SUMMARY 71tie of program: EXPFIT4 Catalogue identifier: ADEZ Program obtainable from: CI:'C Program Libra'y, Que+.'n'sUniversity of Belfast, N. Ireland Licensing provisions: hone Compaters: 386/486/Penlium-based PCs, T800 transputer (MSDOS): Sun Span: (Unix) Operating systems under which the program has been tested: MSDOS, Unix

FORTRAN (version 5.0). SPARCompiler F77, 3L FORTRAN [ I ] ) No. of bytes in distributed program, including test data, eg'.: 140616 Distribution format: ASCII Keywords: Newtonian physics, celeslial mechanic~, ~ ~ linear second-order differengal equa~on,+, ,~2,~ial-vah~ exponential fining, mulfistep metl,,~ls, stcpsiz¢ Nature of physical problem All physical problems which [end themselves to a d c s e ~ by a nonstiff system of linear and/or n
Programmi,+g language used: FORTRAN77 (compilers: MSJ Pennar;¢at t~dd~ss: Institute of Physics and Nuclear Engineering, P.O. Box MG-6, Bin:hadst. Romania. = Research Dip~ctor at the National Fund for Scicmific Research (N.F.W.O.). Belgium. 3 Corresponding aatb,3x;e-mail: [email protected]. 0010-4655/971517.00 Copyfigh: ~ 1997 Llsevie~ ~iet,c~: B.V. All tights reserved. PI! S0010-4655 ( 96 ) 00146-4

72

L.Gr. lxaru et al./Computer Physics Communications i00 fl997) 71-80

for systems of less than I0 equations and if no more than 1000 steps are required.

A [~'vkm~y de~eh~'~ed ~xth-~der four-step exponential-fitted me~t~od [2] is used ;n p~redictor-correc:or mode. The efficiency of ~ : method t~hes on the variability of the stepsize and on the of ~ k~fl mmea~hm eJ~'.

References l l[ Parallel Fortran (version 2.1.3) User Guide. 3L Lid.. UK

T>~'at rg~u~ing time Depem:ls on Lhe sir.e of the probh:m: of lhe ordct" of a few ~conds

t21 L.Gr. i~aru. G. Vanden Berghe. H. De Meyer and M. 'fan Deele. Comput. Phys. Comman. I00 (1997) 56.

(1990).

LONG WRITE-UP L Introduction Many physical problems, such as e.g. pendulum problems, anharmonic oscillator problems, celestial mechanh.~ problems, etc., are such that a system of noniiaear second-order differential equations of the type

J/~=f(x,y),

x E [a,b~,

y ( a ) = y o,

y(a)=y~,

(!.1)

wt'~reby y and f denote real vectors with n components, must be eff~tively solved. Most of the numerical codes for h a ~ l i n g second-order initial-value problems, which are available in existing program libraries, have been designed for general, unspecialised use. In the physical context, howt~.ver, one quite considers problems for which it is known that the solution components exhibit an oscillatory behaviour. For such a class of problems, h is expected that it is beneficial to develop specialized methods. In a previous paper I l L the present authors have constructed an algorithm on the basis of a family of exponential-fitted four-step methods and it has been shown that, indeed, this algorithm performs very well on classical physics problems of the kind mentioned above. 1"i7r~¢program EXPF!"T4 presented in this paper is the FORTRAN implementation of that algorithm. 2. Method of solution In our code F.XPFIT4, the numerical integration of the initial value problem ( 1. I ) for a system of linear or nonlinear second-order differential equations proceeds in a steplike manner by the application of exponentialfitted mul6step formulae. At each pew step an estimation of the error is obtained, from which, by comparison with a gi,~n tolerance, it is decided whether the stepsiz¢ should he maintained or adapted. The code furthermore includes specific procedures to start the integration and to compute the numerical solution in a set of ptcdefined

poims. In what follows, certa/n essential parts of the algorithm will be emphasized. The order of description is, however, not the same as 1.he order in which these parts are executed in the program. A discussion of the flow is postponed until the next section. ( i ) IN the generic case of adva~cing the num,~erical solution one step, five-point (or equivalently, four-step) exponential-fitted methods are used which are of sixth order in the constant stepsize h. They belong to ehe one-parameter class of methods

S(q: z):Y~.2 + al (Yk+l + Y~-I) + a2Yk +Yk-2 = h2 [ ~ (A,2 + A - 2 ) -" b, (A+, + A - , ) + ~ A ] •

(2.1)

The weights az,a2, bo, bl and b2 are calculated upon the coedition that (2.1) integrates exactly the four scalar differential equations y" = ,a~y, i = i,2, 3,4 wP.h given ¢a/2's, either (i) all real or (ii) two real and

L.Gr. Ixaru et al./Computer Physics Communicati~nL~100 (1997) 7;-80

73

one complex conjugate pair or (iii) two complex conjugate pairs° By this c o n d i t ~ , the w ~ g h ~ 2 2 2 2 and one ~ explicitly upon the four-component vector z = h 2 [tou,ta2,ta3,ta4] free ~ q. We adopted the convention that S(0;z) is the method for which a2(0; z) = 0 and S ( | ; z ) is the ~ f~r which at ( 1; z) = 0. (2) For given values of z and q, the problem of finding the weights for S(0; z) and S( l; z) reduces for bo~h metheds to the problem of solving a system of linear equations of the type 4

5-~Ui(Zj)ti='.U(Zj),

j= 1,2,3,4,

(2.2)

/~1 with respect to ti, i = 1,2,3,4, whereby ui, i = ! , 2 , 3 , 4 and v are given functions of one ( c o m ~ c x ) variable and zj,j = l, 2, 3,4 are ~he elements of the parameter-vector z introduced abcrc¢. ~ two or more of the z's are close to each other, the system (2.2) becomes almost singular. In [2} we have established a regularization procedure enabling to generate in all cases the numerical solmh~ with a high degree of accuracy. (3) The calculation of the approximation Y~+2 proceeds by the application of two metluxIs S(q;z) w i ~ the same z but different values of q. A first approximation is obtait~i by a predictol P, which is tee e x ~ t fommla of the type (2.1), but which in general is nonstable. Hence, we further use a c ¢ ~ c t o r C, again of the type (2.1) to generate a second appro×in~ion ofyk+ 2. It follows from stability considerations the q value of the correcter must be chosen in the interva| - 1 < qc < O. In the program, we opted for t ~ choice qc = O, because it is the one which gives rise to the smallest proportionality factor in the e x ~ for the local truncation error. The predictor ~pproximation is used as the starting value in an i~'afion process carried out on the cotr~tor. It consist:; of one Newton-Raphson type iteration if at least one of the z's differs ftom zero (i.e., when the exponential fitting is active) and of one fixed-point i t ~ if -all z's are zero. In the former case, use is made of the Jacobian of the given system ( l . l ) . (4) The z~ values for equivalently of the ~a2 values) from which the predictor and correcter are bui|L a ~ sei~cted on the basis of the eigenvalue spectrum of the Jacobian of the given system ( t . | ) . If not all of the eigenvaiues a~e real and of the same sign, all z's ate set equ'~ to ~ , o and in that case e x ~ . ~ fitting is not activated. If, on the contrary, all eigenvalues are real and of the same sign a~d ~2m the largest eigenvalue in absolute value, the z's are chosen as fcl|o~'s: 2 [ !,0.7,0.3,0]. z = h 2w m (2.3) (5) Once the appt'~ximation to y generated from the con'ector formula is available m seven consecutive points, an estimation ci' the local truncation error (ire) ;s obhaiaed from the expression (3 !q~ - 159)h 2 . -(6, ! t-- " - 2 z m F ( 4 ) + l . 2 1 z 2 F C 2 ) - . 2 1 z 3 f o ) [ ~ ° ,

Ire= I - - - 6 . 9 4 8 0

(2.4)

2 ~ and F(6) , j¢(4) and F(2) ar~ certain linear combinations of the f ' s at those poi.-.L~ whereby zm = tamh"

(s~ [1;). (6) The algorithm includes the possibiiityto change the stepsize on the basis of the comparison of the local truncation error estimation against a given local tolerance. Actually, three ~ypes of stepsize medifica~ioas are foreseen: doubling or halving the CmTCnt stepsize or multiplying the current stepsize by a factor 3/2. The decision whether the step:~ize sheu'd be modified, and if so, what type of modificafiou ~ b~ carried out, proceeds as follows. From the estimation Ftc of ihe local truncation error and the given |oca~ tolerance tot the deviation factor fh := (i.e/tol)-x/s is computed. If f~ < 0.9, the stepsiz¢ is i ~ y halved. Furthermore, if at four consecutive steps 0.9 < f~ < 0.95, the stepsi~ is halved and if, a/so at four consecutive steps, 1.5 < fh < 2, the stepsize is increased by a fa~or 3/2. Finally, if f~ > 2 at f ~ consecutive steps, the stepsize is do~L led. In all othe.~ cases, tb~ st~size is left unchanged.

74

L. Gr. lxaru et aL /Computer Physics Comraunications 100 ~1997) 71-80

(7) Doubling the s ~ s i z e goes without special efforts. On the other hand, when the .~tepsize is halved or muhipJicd by a factor 3/2, speci'aI ,:are has to be taken to generate approximations of y at certain in,.--,~a'~.~d_iaq,• points before the integration can be continued with the new stepsize. The~e approximations are furnished by a suitable linear combination of four-step exponential-fitted methods (algorithm N in [ 1 ]). It turns out that three ad~it,.'cnal function evaluations are required by these types of steFsize modification. (8) The .Iacobian of the system is reevaluated and new values of the z's are introduced after 100 conslant steps and also each time the stepsize is modified. (9) Before the four-step method (2.1) can be applied in a systematic way, approximations (at least of the eighth order in the stepsize) of the solution of the initial value problem (1.1) at tour consecutive equklistant p o i ~ s :co = a and x j = xo + j h , j = i ,2,3, are r,~quired. ~ s d y , an ~proximadon of the solution at xj is obtained from the application of a Runge-KuttaN.,/s~6m one-step formula of sufficiently high order. In |t~,: program, we have opted for the DOPRIN mc~md (see, e.g., [ 3 ] ) with ~he default tolerance 10 -I'~. Next, a p p r o x i ~ i o n s at the points x2 and x3 are obtained by the combination of several two-step and four-step (expo~6'nfial-fitted) methods which are applied with a stepsize being half of the initial stepsiz¢ (see [ | ] ). Once these are known, the general formula (2.1) is applied three times to compute y and f the points x i = :co + j h , j = 4,5,6. From the dam at these seven consecutive equidistant points, the loc~ ernncalion error is estimated and the deviation factor f~ is computed. In the case that fh < 0.95, tbe ShYing procedure is entirely rep~'..ated with the new stepsize equal to 0.9f~ times the initial stepsize. (10) Usually, it is required to genera'.e the numerical solution of (I.1) in a given set of points. ~ince, in general, these points do not coincide with those uzed in the integration process, an interpolation algorithm of degree seven is designed to produce the required approximations.

3. ~

c~scripCon

The program consists of a main routine F.,XPFIT4 and a number of su:~routines, all written in FORTRAN77. ~ F r l ' 4 acts as an internal driver-routine which essentially dictates the order in which algorithmic steps are to be performed. ~ user should however provide the program with three more routines, i.e. * a main program which acts as an external driver and is concerned with input/output operations and settles the values of user-adaptable parameters, such as t.he to|e,'ance, a guess value of the initial stepsize, etc.; * a sabroutine, n~.,,~ [~-kq, which describes the right-hand sides of the system ( !. I); * a subroutine, named SOL, containing the information about the. points at which the numerical solution is

required. Cett~n essential chaxacterisfics of the algorithm, such as the number of iterations prescribed for NewtonRaphson itor'~io~ and the value of q, are fixed by means of program parameters. They are all defined in the subroutine DLP~'IT4 and should not be changed by the user, unless he wants to carry out experiments with modified versions of the presented algorithm. In order to clarify the used conventions, a brief description of these algorithm parameters is given. N r I ~ Number of iterations in the corrector to advance the solution along each new step. This is fixed as NITER ~ 2 if T0L _> 9.10 -3 and as NITI~ = 1. other.vise. I T , D Number of iterations in the corrector when genecating the ar.iditional starting values for stepsize halving or i~'e.asing by a factor 3/2. ~ used value is U~40D ,, 1. Q Value o f q in ~he scheme S ( q ; z ) . It has been set simply as Q = O. The user can consider also sinaller values. However, one should be aware that for Q <_ - ! severe instabilities may appear.

L.Gr. Ixaru et a~ /Computer Physics Communications t00 (1997) 71-80

75

I0PT An input/output parameter with two possible values, 0 and I. IOFr = O (TOFI'. = 1) impJics ~ tb¢ classicel method S(q: O) (the exponential-fitted method S(q; z ) ) is used, with functional ~ . T g c w t o n - ~ ) iteration for the corrector. The initial value is I0PT = t and it will be preserved as long a~ s ~ . ~ q e will decide that the exponential-fitted version is appropriate. If, howes~, at some point along the i ~ FITMET will decide that the classical method is the apwopriate method, its value is modified to IOPT = O. In this case, the solution will be advanced by S(q;0) but checks of w.hett'~r the exponential fitting ~ be activated again will be done at each NSTgP which is a multiple of 100. inv~lw,~i subroutines c~n be gro,Jped as follows: ( 1) A subroutine which allows timing tim vector z in terms of the given problem: F I ~ (which, at its tin'n, calls JACOB and GF~Z). (2) Subroutines which, once z is known, furnish the weights of the involved two- and four-stq~ m ¢ ~ . ~ : F.,FN3, EFM5 (which call CGEBAS, GDF3, GDF5, GEWCand CWGHT)and CORNET. (3) Subroutines to start the solution: STAIITand DOPlgIN. (4) Subroutines to advance the solution one step: ADV which, at its turn, calls CORRECT. (5) Subroutin~ [or the ste~ize modification: PRELMOD(calling STOM, WHALVE)and HIDSOL. (6) Auxiliary, mostly matrix operation subroutines. A short description of all subroutines is given below. ADV Allows computing whe wlues at x of y and f when the values at x - j h , j = 1,2, 3, 4 are known. COF_.BAS Calculates, foc given complex z (double complex variable g), the values of ~ ( z ) , f l o ( z ) , f l l ( z ) , . . . . ~ ( z ) . The ~j's, detin~d in [4], are proportional to successive derivatives of ~(z). These data are used to construct and regularize the algebraic linear systems which yield the weights of the involved two- and four-step metlmds. COPY Furnishes a copy B of a given rectangular matrix A. CORRECT Iterates for the corrector in terms of the input parameters IOPT and NITL'~. b"t~lgT Computes, for given q, the weights of S ( q ; z ) in terms of the weights of S(0;z) and S ( l ; z ) . ~GtrI' Helps regularizing the algebraic linear system of equations for the method weights. El~r~3 Calculates the weights of the two.step method~; these are necessary for the starting procedure. Twocomponent vectors P~1., P~0 and 1tB1 give the values of the (real) weights at,bo and b~ in Eq. (4.1) of [ I ] as follows: if or, input, INDEg ffi 1, ILA.I(M), RBO(M) and ItBI(M) give the weights of Pe for M = 1 and of P~ fer M = 2; if INDEX ffi 2, the weights of C2 are given at M = 1. Remark I. For INDEX = 1, only the first two double complex type elements g(1) and 2(2) ~a'¢ used. They

must be either both real or a pair of complex conjugate numbers. For DJDF.g ffi 2, three values Z ( 1 ) , Z(2) and Z(3) are used. They must be either all real, o:' one real and a pair of complex conjugate vales. In these cases the r~ultant weights arc r~al. EFM5 Calculates the weights of S(0;z) and S ( I , z ) for given complex vector z (double complex variable g of dimension 4). The resultant weights are found in 1~0(:[~), Pd,l ( I Q ) , PA2(I~), RB0(T~), RBI(I~) and RB2(I~), with IQ either 0 or 1. Remark 2. The l~burcomponents of the input vector Z must he eitbe~ all :'¢~1, two real and a pair of complex

conjugate values, or two pairs of complex conjugate values. Olllv i,a these cases the resulta~-~t weighLs are real. EIG Acts as an interf~.ce routine for the calculation of the eigenvalues of the J~obian raatrix by ~:alling in sequence appropriate routines taken from the Numerical Recipes L i b ~ [ 5 ].

L.Gr. lxaru er al. /Computer Physics Communications 100 (1997) 71-80

E g ~ r l ' 4 Driver romine which contains the definition of algorithm parameters and controls the order of the FTrIMI~ Basically allows fitting the parang~exs of the method in terms of the form o f f . To this aim R commands the ~ h ~ of the .;acobian, the calculation of its eigenvalues and, on this basis, the generation of the doub~ c o m p ~ vector g with fi~u" elements. C~culafes the coefficients of the linear system for the two-step method weights and also their first cigh~. derivatives. Does tim same as G'DF3 for ~he co¢.ffici¢n~s of the linear system for the four-step method weights. G~¢Z G e n e r a ~ the values of zi, i = 1,2, 3,4 in the restrictive case when the eigenvalues of the .Iacobian are real and of tim santo sign. Decides on w ~ the linear system should he regularized and, ,~fso, performs the. regulanzation. I~P Comtructs ~ ¢ in~'polation polynomial of seventh dcgrc¢ through four given equidistant points such that at these points the polynomial and its second derivative coincide with y ( x ) and f ( x , y), respectively. 3AC~ Consm~tg a numerical approximation of the Jacobian matrix in terms of the right-hand side of the given syslcm of differential equations. LI/F2 C ~ s , on the basis of the Jacobian of the given system, the LU-decompos~ nuarix to b¢ used in the Newton-Raphson iteration process for solving the corrector equation with respect to Yk+2 (sea Eq. (2.1)). MID~L M/dpoim exl~nential-fitted interpolation subroutine. I ~ Keeps track of tim values of f at a succession of steps. It actually allows having the values of f stored at the last thirteen points. This information is further used in F,XPFIT for both calculating the error and ~ n g the starting dam when the stepsize is halved, increased by a factor 3/2 or doubled, as dictated by" the error cc~.'.,"o|. Compmes tim mauix-product A. B of two given rectangular matrices A and B. Congruous the matrix A + t . B for given scalar t and two given rectanguiar matrices A and B. Congress the w~,'ights of tim N algorithm. The weights AHP(I), BliP(I) (for tim predictor) and AHC(I), BHC(1) (for the corrector), I=0,1,2,3,4, are conununicated to MIDSOL by COMMON/COMID/. b'~AlgT C a l c u l a ~ the additional s~rting valu~ Y2 and y~ by the five-stage procedure describe/J in [ ! ]. A subrou~n¢ called by P~I/:LI~D to produce the individual compormnts of the algorithm N. L~FI~ Gen~'at~ a rectangular matrix of given dimensions with unit diagonal elements and all other elements ~o. WHAL~R~Solves the linear system which results from the five conditions imposed on algorithm N. C~nerat~s a rectangular zero-matrix of given dimensions. The ~ uses fur'd~exthe Numerical Recipes [5] subroutines BAL~C, E&MH~, HQILLUDC~ and Lb~KSB. complex afithn~ic is required, the last two subreutines are replaced by CLUDC~ and CLUBKSB,r~p¢cfively. Also, ~ program uses the subroutine DOPRIN (and the assoc~ate, i dummy subroutine 30LIYl'2) described ~n [3]. The las~ mentioned codes are published and belong to the pubhc domain software. 3. L The driver program and other routines provided by the user

The following ¢x~anpl¢ of a driver program and of the two other sL,broutines to he written by the user contains all infom~ion required to write own applications. _*2 C C C

:mI$ I$ THE;DgZVlN~;PRI~PJ~I'1"0TESTEXP!:IT4. THE TESTCASEI5 THE ~;~gt~ P[~UL~ PROBLEM.

DI~$1~

Y(255) ,F(R~)

L.Gr. lxartt et al./Computer Physics Communications i00 (1997) 71-80 C O ~ N / C O ~ O L / X H I N , ~ , IOONTO,ICOr~

CO~NISU~I~

, ~

,NFAZL

C

C

THE PARAMETERS XNIN, X~iX, INPUT TOLERANCE (GIVTOL), NU~BEP.OF

C

F~ATIONS (NF~), THE VALUES AT XNIN OF THE SOLUTION (Y(1 ),

C

Y(2) . . . . Y(NEQ)) AND OF ITS FIRST DERIVATIVE (F(1),F(2) . . . . . F(NEQ))

C C C C

AI~D A GUESS VALUE (HIN) PO~ THE INITIAL STEPSIZE FOLLOW. NOTE THAT THE DIMENSION OF VECT0mS Y AND F SHOULD BE 5 TIMES NP, WHERENP HAS THE VALUE 51, AS DECLARED IN SBE.EXPFZT4. XMIN-O. DC XMAX'20. DO

GIVTOL- 1. i)-7 NF.~2 Y(1)al.5DO YC2)-O.DO F(1)-O.DO F(2)"1 .DO ltLN,,. IDO C C C C C C C C

SUBROUTINESOL, WRITTEN S~ARATELY,, ALLOWS PAINTING THE SOLUTION AT INTERMEDIARYHWR SPACED POINTS. THE VALUES OF HWR, ICONTO AND ICONTNFOLLOW. THEY WILL ALLOWWRITING THE SOLUTION VALUES AT THE EQUIDISTANT POINTS X-XHIN+MeH~, M~ICONT0+I. . . . . ICONTN. T H E S EDATA ARE TRANSMITTED TO SOL THROUGH ~ / C O M S O L / RESTRICTION: (ICOMT0+I)*INR SHOULD BE GREA.~ER THAN 2*HIN. INlt-.ID0 ICO~rl'O-5 ICONTN-25 CALL F..XPFXT4(NEQ,XHIN, Xl~LqX,NHS,Y, F, GI%~k'OL,HIN)

C C C

50 C C C 51

THE SOLUTION AT XMAX I~ RETURNED IN VECTOR Y. IT IS WRITTEN BELOW.

WRTTE(*,50)XHAX, (¥(N) ,N:,I,~Q) FORMAT(/,2X,'X - ' , P 6 . 3 , 2 X , ~ Y - ' , 2 ( 2 X , F 1 0 . 6 ) ) FINALLY, THE STATISTICAL PARA~IErERS NFEV,NSTEP AND ~FAIL

ARE WRITTEN.

WRITE(*, 51) NFEV,NSTEP, NFAIL FOI~4AT(/,2X, 'NFEV-' ,14, ' NSTEP-', 14, ' NFAIL-' ,I3) STOP END

C SUBROUTINe: ILHS(NEll, X,Y, F)

C C C C C

TiilS FURNISHES THE It.H.S. 0P THE SPRING PENDULUM SYSTEM OF ETJATIONS. CONMON/CONT/NgEVIS ADDED TO KEEP TRACK FOR THE NUHBER OF ,~'NCTION EVALUATIONS. IMPLICIT DOUBLE PRECISION(A-H,0-Z)

DIMENS!0N F(1).Y(!) COMMON/CONT/NFEV

N~-~FEV+I D~,DSQ~LT(Y(1)*Y(1)+Y(2)*Y(2)) F (1)-20. DO-40.DO*Y(1)* (1.00-1. DO/D) F(2),,-40. DO*Y(2)* (1.D0-1 DO/D)

78

L.Cr lxaru et cl./Computer Physics Communications I00 (1997) 71-80 P~CTVgN F2CD

C C C C C C C C

THIS IS A SL~R~JTINE C A L L ~ !~ F~PFIT4,, hH~CH ALLOWS KNOWING THE SOLVT~ON AT I ~ I A g Y PCI~TS. IF 1-fl~ IS NOT NECESSARY, WRITE IT SI~LOLT AS A D~IMY SD~ROUTIh~. AS IT I~ WRITTEN NOW, IT ~ D ~ I S ~ S THE SOLLerIO~ AT EO~!IDIdT&NT POINTS WITH STFPSIZE HWR AS ~2R~ISHED I~ MAI~ PROGraM Ah~D TRANSMTTTED THROUGH COMMON/COF~0L/. IMPLICIT DOU~LE PRECISION(A-H,0-Z~ DI~E~SIO~ YCNP, I) ,F(~P, I) DI~SI~ Y'~(51) C~/CO~SOL/~IN, ~ , ICONTO, ICO~TN

C

I

X~ffiXHI~ ICO~'=C gh~gffiX~W~ ICO~T=ICO~T+ I IF(XW~. ~T.X) RETURN IF(X~.LT.X-2.DO*H) GO IF( ICO~T. L~. ICO~ITO)R ~ IF(IC~T. GT. I C O n ) R ~

TO 1

DX=X~fl~-X÷H CALL T ~ ( D X , ~ I ~ L ( ~ ) ,¥(~, i ) ,Y(~,2) ,Y(N ,3) ,Y(~,4) ,F(N, I) ,F(N,2), 2 50

CO~TI~E ~ITE(*,50)XWR~ (YWR(~),~=1 ,~EQ) FORMAT(2X,'X =',F~.2,2X,'SOL =',2(2X,F8.4)) ICC~O=ICO~T g~dR~ END

4. ~

n

The progrdm EXPFIT4 is written in, FORTRAN and can be run on any system providing a FORTRAN compiler. We have tested the program on a range of m~'chines (such ,as 386/486/Pentium-based l:'Cs running MS-DOS, a Sun Sparc running Unix, and a single T800 transputer hosted by a PC running MS-DOS) using several compilers (such as MS-FORTRAN, SPARCompiler 1=77, and 3L Parallel FORTRAN [6] ). The method was ~ested on a large sample of problems mostly originating from physics and particularly from classical mechanics. For each problem, the program was repeatedly executed with given tolerances ran~.ng betw~n 10 .-4 and l0 - l ° and we kep~ track of the total number of steps and the total number of function cvalu~ions nfev required to integrate t~he system of differential equations up to a given endpoim. Wc. ~so c;-~h:~cd the absolute error err at th~s endpoint and the obtained results were represented as a curve in the l o g n f e v vs l o g e r r plane. We furtber c~mpared our sixth-order method with a popular method of seventh order, the Runge-Kut~-Nystr~m based routine D O P R I N [ 3]. all u:st problems, we noticed that EXPFIT4 performs significantly better than DOPRIN if the total n v r ~ r of function evaluations required to obtain the same absolute error is takcn as the main efficiency criterion (see also [ 1 ] where details about six representative test problems are presented). It is. however, fair adding that in noise of these test cases the vector function f ( x , y ) is complicated enough ~o dominate the computation. We

L.Gr. Ixaru et al. /Computer Physics Communications I00 (1997) 71-80

then made a separate investigation of the relative costs for those test cases to find that (a) The cost per step for F..XPFIT is typically smaller than for DOPRIN, namely at most 50% if cn~ly the solution is advanced and at most 90% if the solution is advanced and the local error is c o h e r e d . (b) The computation of the method weights implies some costly operations, among them the compuu~km and diagonalization of the Jacobian. However, these extra operations are activated w~st|y when the stepsize i~ modified and this happens too seldomly to change the total computation time significantly. (c) For one and the same input accuracy, the number of steps required by F..XPFIT4 is signific~mly larger then with DOPRIN (sometimes up to a factor of four) so that DOPRIN is faster by a factor of up to tJn'ee.

Acknowledgements The authors are grateful to the referees for their valuable comment~ end s~ggestions. One of the autho~ (L. [xaru) wants to express his gratitude to the Belgian National Fund for Scientific Research (N.F.W.O.) for the financial support which allowed him to stay at the Department of Applied Mathematics and Computer Science of the University of Gent.

References I ! I L.Gr. ixaru, G. Vanden Berghe, H. De Meyer and M. Van. Dne~, Comput. Phys. Commun. |00 (1997) 56. 12] L.Gr. lxam, H. De Meyer, G. Vanden Berghe and M. Van D~le, Numer. Lin. AIg. Appl. 3 (1096) 81.

[3] E. Hatter, SE NOr~tt and G. Wanner,SolvingOrdinaryDifferentialEquations I, Non~iff Problems (Springer, Berlin, 1987). [4! L.Gr. lxam, NLmericalMethods for DifferentialEquations (Reidel. Dofdrecht, 1984). ;5] W.H. Press. B.P. Fiannery,S.A. Teukolskyand W.T. Vetterling,Numerical Recipes - the Art of ScientificComp~ing FORTRAN (Cambridge Univ [hess, Cambridge, 1986). [6] Parallel Fortran (version 2.1.3) User Guide, 3L Ltd., UK (1990).

L.Gr. I~ar. ei t:L /Computer Physics Commu._ic~ticr.s 100 {1997) 71-~,0

80 TE.g][" R U N O U T P U T X = 0.60

SOL =

1.4738

0.2190

X = 0.70

SOL =

1.4731

0.~47~

X = 0.80

SOL =

~.4808

0.0565

X = 0.90

SOL =

1.4956

-0.0415 -0.1340

X = 1.00

SOL =

1.5117

X = 1.10

SOL =

1.52!8

-0.2084

X = 1.20

SOL =

1.5199

-0.2540

X = 1.30

SOL =

1.5051

-0.2644

X = 1.40

SOL =

1.4826

-0.2388

X = 1.50

SOL =

1.4620

-0.1816

X = I.o60

SOL =

1.4532

-0.1012

X = 1.70

SOL =

1.4614

-0.0082

X = !.80

SOL =

1.4844

0.0858

X = 1.90

SOL =

1.5127

0.1686

X = 2,00

SOL =

1.5337

0.228.4

X = 2.10

SOL =

1.5375

0.2562

X = 2.20

SOL =

1.5216

0.2478

X = 2.30

SOL =

1.4925

0.2051

X = 2.40

SOL =

1.4627

0.1350

X = 2.50

SOL =

1.4456

0.0477

X =20.000

Y =

-0.257909

1.531259

~"EV= 879 NSTEP= 672 NFAIL=

3