NUMAL Section 5.2.1.1.1.1.G
BEGIN SECTION : 5.2.1.1.1.1.G (February, 1979)
PROCEDURE : DIFFSYS.
AUTHORS : R.BULIRSCH AND J.STOER.
CONTRIBUTOR: K.DEKKER.
INSTITUTE: MATHEMATICAL CENTRE.
RECEIVED: 731231.
BRIEF DESCRIPTION:
DIFFSYS SOLVES AN INITIAL VALUE PROBLEM ,FOR A SYSTEM OF FIRST
ORDER ORDINARY DIFFERENTIAL EQUATIONS DY / DX = F(X,Y).
THE METHOD IS RECOMMENDED IF HIGH ACCURACY IS DESIRED.
DIFFSYS IS NOT SUITED FOR STIFF EQUATIONS.
KEYWORDS:
INITIAL VALUE PROBLEMS,
SYSTEM OF FIRST ORDER ORDINARY DIFFERENTIAL EQUATIONS.
CALLING SEQUENCE:
THE HEADING OF THE PROCEDURE DIFFSYS READS:
"PROCEDURE" DIFFSYS(X,XE,N,Y,DERIVATIVE,AETA,RETA,S,H0,OUTPUT);
"VALUE" N;
"INTEGER" N;
"REAL" X,XE,AETA,RETA,H0;
"ARRAY" Y,S;
"PROCEDURE" DERIVATIVE,OUTPUT;
"CODE" 33180;
THE MEANING OF THE FORMAL PARAMETERS IS:
X: <VARIABLE>;
THE INDEPENDENT VARIABLE ;
ENTRY: THE INITIAL VALUE X0;
EXIT : THE FINAL VALUE XE;
XE: <ARITHMETIC EXPRESSION>;
THE FINAL VALUE OF X (XE>=X);
N: <ARITHMETIC EXPRESSION>;
THE NUMBER OF EQUATIONS;
Y: <ARRAY IDENTIFIER>;
"REAL" "ARRAY" Y[1:N];
THE DEPENDENT VARIABLE;
ENTRY: THE INITIAL VALUES OF THE SYSTEM OF DIFFERENTIAL
EQUATIONS AT X=X0;
EXIT : THE FINAL VALUES OF THE SOLUTION AT X=XE;
DERIVATIVE: <PROCEDURE IDENTIFIER>;
THE HEADING OF THIS PROCEDURE READS:
"PROCEDURE" DERIVATIVE(X,Y,DY); "REAL" X; "ARRAY" Y,DY;
THIS PROCEDURE SHOULD DELIVER THE RIGHT HAND SIDE OF
THE I-TH DIFFERENTIAL EQUATION AT THE POINT (X,Y) AS DY[I],
I=1,...,N;
AETA: <ARITHMETIC EXPRESSION>;
REQUIRED ABSOLUTE PRECISION IN THE INTEGRATION PROCESS;
RETA: <ARITHMETIC EXPRESSION>;
REQUIRED RELATIVE PRECISION IN THE INTEGRATION PROCESS;
S: <ARRAY IDENTIFIER>;
"REAL" "ARRAY" S[1:N];
THE ARRAY S IS USED TO CONTROL THE ACCURACY OF THE COMPUTED
VALUES OF Y;
ENTRY: IT IS ADVISABLE TO SET S[I]=0, I=1,...,N;
EXIT : THE MAXIMUM VALUE OF ABS(Y[I]), ENCOUNTERED DURING
INTEGRATION, IF THIS VALUE EXCEEDS THE VALUE OF S[I]
ON ENTRY;
H0: <VARIABLE>;
THE INITIAL STEP TO BE TAKEN;
OUTPUT: <PROCEDURE IDENTIFIER>;
THE HEADING OF THIS PROCEDURE READS;
"PROCEDURE" OUTPUT;
THIS PROCEDURE IS CALLED AT THE END OF EACH INTEGRATION
STEP ; THE USER CAN ASK FOR OUTPUT OF SOME PARAMETERS , FOR
EXAMPLE X, Y, S.
PROCEDURES USED: NONE.
REQUIRED CENTRAL MEMORY : CIRCA 28* N MEMORY PLACES.
METHOD AND PERFORMANCE:
THE PROCEDURE DIFFSYS IS A SLIGHT MODIFICATION OF THE ALGORITHM
PUBLISHED BY BULIRSCH AND STOER (SEE REF[1]) . BY THIS MODIFICATION
INTEGRATION FROM X0 UNTIL XE CAN BE PERFORMED BY ONE CALL OF
DIFFSYS. A NUMBER OF INTEGRATION STEPS ARE TAKEN, STARTING WITH THE
INITIAL STEP H0. IN EACH INTEGRATION STEP A NUMBER OF SOLUTIONS ARE
COMPUTED BY MEANS OF THE MODIFIED MIDPOINT RULE . EXTRAPOLATION IS
USED TO IMPROVE THESE SOLUTIONS,UNTIL THE REQUIRED ACCURACY IS MET.
AN INTEGRATION STEP IS REJECTED, IF THE ACCURACY REQUIREMENTS ARE
NOT FULFILLED AFTER NINE EXTRAPOLATION STEPS . IN THESE CASES THE
INTEGRATION STEP IS REJECTED , AND INTEGRATION IS TRIED AGAIN WITH
THE INTEGRATION STEP HALVED.
THE ALGORITHM IS FOR EACH STEP A VARIABLE ORDER METHOD (THE HIGHEST
ORDER IS 14 ), AND USES A VARIABLE NUMBER OF FUNCTION EVALUATIONS,
DEPENDING ON THE ORDER (MINIMUM IS 3, MAXIMUM IS 217).
THE ALGORITHM IS LESS SENSITIVE TO TOO SMALL VALUES OF THE INITIAL
STEPSIZE THAN THE ORIGINAL ALGORITHM ( SEE REF [2] ); HOWEVER BAD
GUESSES REQUIRE STILL SOME MORE COMPUTATIONS.
REFERENCES:
[1]. R.BULIRSCH AND J.STOER.
NUMERICAL TREATMENT OF ORDINARY DIFFERENTIAL EQUATIONS BY
EXTRAPOLATION METHODS.
NUMERISCHE MATHEMATIK, VOLUME 8, PAGE 1-13, 1965.
[2]. PHYLLIS FOX.
A COMPARATIVE STUDY OF COMPUTER PROGRAMS FOR INTEGRATING
DIFFERENTIAL EQUATIONS.
COMMUNICATIONS OF THE A.C.M., VOLUME 15, PAGE 941-948, 1972.
[3]. T.E.HULL, W.H.ENRIGHT, B.M.FELLEN AND A.E.SEDGWICK.
COMPARING NUMERICAL METHODS FOR ORDINARY DIFFERENTIAL
EQUATIONS.
SIAM JOURNAL ON NUMERICAL ANALYSIS,VOLUME 9,PAGE 603-635,1972.
EXAMPLE OF USE:
THE FOLLOWING PROGRAM ILLUSTRATES THE COSTS AND THE ACCURACIES
WHICH ARE OBTAINED WHEN SOLVING A SYSTEM OF DIFFERENTIAL EQUATIONS
ARISING FROM THE RESTRICTED PROBLEM OF THREE BODIES ( SEE REF[1] ).
THE SOLUTION IS A CLOSED ORBIT WITH PERIOD T=6.192169331396.
"BEGIN"
"INTEGER" PASSES,K;
"REAL" X,XE,TIME,TOL,H0;
"REAL" "ARRAY" Y,S[1:4];
"PROCEDURE" DER(X,Y,DY); "REAL" X; "ARRAY" Y,DY;
"BEGIN" "REAL" MU,MU1,Y1,Y2,Y3,Y4,S1,S2;
MU:=1/82.45; MU1:=1-MU;
PASSES:=PASSES+1;
Y1:=Y[1]; Y2:=DY[1]:=Y[2]; Y3:=Y[3]; Y4:=DY[3]:=Y[4];
S1:=(Y1+MU)**2+Y3**2; S2:=(Y1-MU1)**2+Y3**2;
S1:=S1*SQRT(S1); S2:=S2*SQRT(S2);
DY[2]:=Y1+2*Y4-MU1*(Y1+MU)/S1-MU*(Y1-MU1)/S2;
DY[4]:=Y3-2*Y2-MU1*Y3/S1-MU*Y3/S2
"END";
"PROCEDURE" OUT;
"BEGIN" K:=K+1;
"IF" X>=XE "THEN"
OUTPUT(61,"("2(-5ZD),2(4B+Z.3DB3DB3DB3D),-5ZD.3D,/")",K,
PASSES,Y[1],Y[3],CLOCK-TIME)
"END";
OUTPUT(61,"(""(" THIS LINE AND THE FOLLOWING TEXT IS ")"
"("PRINTED BY THIS PROGRAM")",//,
"(" THE RESULTS WITH DIFFSYS - H0=.2 - ARE: ")",/,
"(" K DER.EV. Y[1] Y[3] ")",
"(" TIME")",/")");
"FOR" TOL:="-4,"-6,"-8,"-10,"-12 "DO"
"BEGIN" PASSES:=K:=0; X:=0; XE:=6.192169331396;
Y[1]:=1.2; Y[2]:=Y[3]:=0; Y[4]:=-1.04935750983;
S[1]:=S[2]:=S[3]:=S[4]:=0; H0:=.2; TIME:=CLOCK;
DIFFSYS(X,XE,4,Y,DER,TOL,TOL,S,H0,OUT);
"END"
"END"
THIS LINE AND THE FOLLOWING TEXT IS PRINTED BY THIS PROGRAM:
THE RESULTS WITH DIFFSYS - H0=.2 - ARE:
K DER.EV. Y[1] Y[3] TIME
30 2591 +1.320 357 347 741 -.032 645 454 836 5.686
33 3414 +1.200 078 037 878 -.000 053 906 067 7.455
37 4213 +1.200 003 282 801 -.000 002 363 741 9.267
44 4618 +1.199 999 999 711 -.000 000 000 095 10.242
56 6299 +1.200 000 000 003 -.000 000 000 090 13.827
SOURCE TEXT:
"CODE" 33180;