NUMAL Section 5.2.1.1.1.1.I
BEGIN SECTION : 5.2.1.1.1.1.I (February, 1979)
PROCEDURE : EFRK.
AUTHOR: K.DEKKER.
INSTITUTE: MATHEMATICAL CENTRE.
RECEIVED: 740710.
BRIEF DESCRIPTION:
EFRK SOLVES AN INITIAL VALUE PROBLEM FOR A SYSTEM OF FIRST
ORDER ORDINARY DIFFERENTIAL EQUATIONS DU / DT = H(T,U) .
EFRK IS A SPECIAL PURPOSE PROCEDURE FOR STIFF EQUATIONS WITH A
KNOWN, CLUSTERED EIGENVALUE SPECTRUM.
KEYWORDS:
INITIAL VALUE PROBLEM,
SYSTEM OF FIRST ORDER ORDINARY DIFFERENTIAL EQUATIONS,
STIFF EQUATION.
CALLING SEQUENCE:
THE HEADING OF THE PROCEDURE EFRK READS:
"PROCEDURE" EFRK ( T,TE, M0,M, U, SIGMA, PHI, DIAMETER, DERIVATIVE,
K, STEP, R, L, BETA, THIRDORDER, TOL, OUTPUT );
"VALUE" R,L;
"INTEGER" M0,M,K,R,L;
"REAL" T,TE,SIGMA,PHI,DIAMETER,STEP,TOL;
"ARRAY" U,BETA;
"BOOLEAN" THIRDORDER;
"PROCEDURE" DERIVATIVE,OUTPUT;
"CODE" 33070;
THE MEANING OF THE FORMAL PARAMETERS IS:
T: <VARIABLE>;
THE INDEPENDENT VARIABLE T;
ENTRY: THE INITIAL VALUE T0;
EXIT : THE FINAL VALUE TE;
TE: <ARITHMETIC ETPRESSION>;
THE FINAL VALUE OF T (TE>=T);
M0: <ARITHMETIC EXPRESSION>;
THE INDEX OF THE FIRST EQUATION;
M: <ARITHMETIC EXPRESSION>;
THE INDEX OF THE LAST EQUATION;
U: <ARRAY IDENTIFIER>;
"REAL" "ARRAY" U[M0:M];
THE DEPENDENT VARIABLE;
ENTRY: THE INITIAL VALUES OF THE SOLUTION OF THE SYSTEM OF
DIFFERENTIAL EQUATIONS AT T = T0;
EXIT : THE VALUES OF THE SOLUTION AT T = TE;
SIGMA: <ARITHMETIC EXPRESSION>;
THE MODULUS OF THE POINT AT WHICH EXPONENTIAL FITTING IS
DESIRED , FOR EXAMPLE AN APPROXIMATION OF THE CENTRE OF THE
LEFT HAND CLUSTER;
PHI: <ARITHMETIC EXPRESSION>;
THE ARGUMENT OF THE CENTRE OF THE LEFT HAND CLUSTER; IN THE
CASE OF TWO COMPLEX CONJUGATED CLUSTERS , THE ARGUMENT OF
THE CENTRE IN THE SECOND QUADRANT SHOULD BE TAKEN;
DIAMETER: <ARITHMETIC EXPRESSION>;
THE DIAMETER OF THE LEFT HAND CLUSTER OF EIGENVALUES OF THE
JACOBIAN MATRIX OF THE SYSTEM OF DIFFERENTIAL EQUATIONS;
IN CASE OF NON-LINEAR EQUATIONS DIAMETER SHOULD HAVE SUCH A
VALUE THAT THE VARIATION OF THE EIGENVALUES IN THIS CLUSTER
IN THE PERIOD ( T ,T+STEP ) IS LESS THAN HALF THE DIAMETER;
DERIVATIVE: <PROCEDURE IDENTIFIER>;
THE HEADING OF THIS PROCEDURE READS:
"PROCEDURE" DERIVATIVE(T,U); "REAL" T; "ARRAY" U;
THIS PROCEDURE SHOULD DELIVER THE VALUE OF H(T,U) IN THE
POINT (T,U) IN THE ARRAY U;
K: <VARIABLE>;
COUNTS THE NUMBER OF INTEGRATION STEPS TAKEN;
FOR EXAMPLE, K MAY BE USED IN THE EXPRESSION FOR TE;
ENTRY: AN (ARBIRARY) CHOSEN VALUE K0, E.G. K0=0;
EXIT : K0 + THE NUMBER OF INTEGRATION STEPS PERFORMED;
STEP: <ARITHMETIC EXPRESSION>;
THE STEPSIZE CHOSEN WILL BE AT MOST EQUAL TO STEP ;
THIS STEPSIZE MAY BE REDUCED BY STABILITY CONSTRAINTS,
IMPOSED BY A POSITIVE DIAMETER , OR BY CONSIDERATIONS OF
INTERNAL STABILITY (SEE REF[1], PAGE 11);
R: <ARITHMETIC EXPRESSION>;
R + L: THE NUMBER OF EVALUATIONS OF H(T, U) ON WHICH THE
RUNGE-KUTTA SCHEME IS BASED;
FOR R=1,2,>=3 FIRST, SECOND AND THIRD ORDER ACCURACY MAY BE
OBTAINED BY AN APPROPRIATE CHOICE OF THE ARRAY BETA;
L: <ARITHMETIC EXPRESSION>;
ENTRY:
IF PHI = 4*ARCTAN(1): THE ORDER OF THE EXPONENTIAL FITTING,
ELSE TWICE THE ORDER OF THE EXPONENTIAL FITTING;
NOTE THAT L SHOULD BE EVEN IN THE LATTER CASE;
BETA: <ARRAY IDENTIFIER>;
"REAL" "ARRAY" BETA[0:R+L];
ENTRY: THE ELEMENTS BETA[I] , I=0, ... ,R SHOULD HAVE THE
VALUE OF THE R+1 FIRST COEFFICIENTS OF THE STABILITY
POLYNOMIAL;
THIRDORDER: <BOOLEAN EXPRESSION>;
IF THIRD ORDER ACCURACY IS DESIRED , THIRDORDER SHOULD HAVE
THE VALUE "TRUE" , IN COMBINATION WITH APPROPRIATE CHOICES
OF R (R>=3) AND THE ARRAY BETA ( BETA[I]=1/I!, I=0,1,2,3 );
IN ALL OTHER CASES THIRDORDER MUST HAVE THE VALUE "FALSE";
TOL: <ARITHMETIC EXPRESSION>;
AN UPPERBOUND FOR THE ROUNDING ERRORS IN THE COMPUTATIONS
IN ONE RUNGE-KUTTA STEP ; IN SOME CASES ( E.G. LARGE VALUES
OF SIGMA AND R ) TOL WILL CAUSE A DECREASE OF THE STEPSIZE;
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 T, K, U, AND COMPUTE NEW VALUES FOR SIGMA, PHI, AND
DIAMETER.
PROCEDURES USED:
INIVEC= CP31010,
ELMVEC= CP34020,
DEC = CP34300,
SOL = CP34051.
REQUIRED CENTRAL MEMORY : CIRCA 30 + (M - M0) + L * (5 + L).
METHOD AND PERFORMANCE:
EFRK IS BASED ON EXPLICIT RUNGE-KUTTA METHODS OF ORDER 1, 2 AND 3,
WHICH MAKE USE OF EXPONENTIAL FITTING. AUTOMATIC ERROR CONTROL
IS NOT PROVIDED.
A DETAILED DESCRIPTION OF THE METHOD AND SOME NUMERICAL EXAMPLES
ARE GIVEN IN REF[1]. REF [3], PAGE 170 REPRESENTS A BRIEF SURVEY. A
COMPARATIVE TEST OVER A LARGE CLASS OF DIFFERENTIAL EQUATIONS IS
GIVEN IN REF [4].
FROM THESE RESULTS IT APPEARS THAT CALLS WITH THIRDORDER = "TRUE"
ARE LESS ADVISABLE.
REFERENCES:
[1]. K. DEKKER.
AN ALGOL 60 VERSION OF EXPONENTIALLY FITTED RUNGE-KUTTA
METHODS (DUTCH).
NR 25 (1972), MATHEMATICAL CENTRE.
[2]. T. J. DEKKER, P. W. HEMKER AND P. J. VAN DER HOUWEN.
COLLOQUIUM STIFF DIFFERENTIAL EQUATIONS 1 (DUTCH).
MC SYLLABUS 15.1, (1972) MATHEMATICAL CENTRE.
[3]. P. A. BEENTJES, K. DEKKER, H. C. HEMKER, S.P.N. VAN KAMPEN AND
G. M. WILLEMS.
COLLOQUIUM STIFF DIFFERENTIAL EQUATIONS 2 (DUTCH).
MC SYLLABUS 15.2, (1973) MATHEMATICAL CENTRE.
[4]. P. A. BEENTJES, K. DEKKER, H. C. HEMKER AND M. V. VELDHUIZEN.
COLLOQUIUM STIFF DIFFERENTIAL EQUATIONS 3 (DUTCH).
MC SYLLABUS 15.3, (1974) MATHEMATICAL CENTRE.
EXAMPLE OF USE:
CONSIDER THE SYSTEM OF DIFFERENTIAL EQUATIONS:
DY[1]/DX = -Y[1] + Y[1] * Y[2] + .99 * Y[2]
DY[2]/DX = -1000 * ( -Y[1] + Y[1] * Y[2] + Y[2] )
WITH THE INITIAL CONDITIONS AT X = 0:
Y[1] = 1 AND Y[2] = 0. (SEE REF[2], PAGE 11).
THE SOLUTION AT X = 50 IS APPROXIMATELY:
Y[1] = .765 878 320 487 AND Y[2] = .433 710 353 5768.
THE FOLLOWING PROGRAM SHOWS SOME DIFFERENT CALLS OF THE PROCEDURE
EFRK:
"BEGIN"
"INTEGER" K,R,L,PASSES;
"REAL" X,SIGMA,PHI,TIME,STEP,DIAMETER;
"REAL" "ARRAY" Y[1:2],BETA[0:6];
"PROCEDURE" DER(X,Y); "REAL" X; "ARRAY" Y;
"BEGIN" "REAL" Y1,Y2; Y1:=Y[1]; Y2:=Y[2];
Y[1]:=(Y1+.99)*(Y2-1)+.99;
Y[2]:=1000*((1+Y1)*(1-Y2)-1);
PASSES:=PASSES+1
"END";
"PROCEDURE" OUT;
"BEGIN" "REAL" S;
S:=(-1000*Y[1]-1001+Y[2])/2;
SIGMA:=ABS(S-SQRT(S*S+10*(Y[2]-1)));
DIAMETER:=2*STEP*ABS(1000*(1.99*Y[2]-2*Y[1]*(1-Y[2])));
"IF" X=50 "THEN"
OUTPUT(61,"("4BD,2BD,2(-5ZD),2(4B+.3DB3DB3D),-5ZD.3D,/")",
R,L,K,PASSES,Y[1],Y[2],CLOCK-TIME)
"END";
OUTPUT(61,"(""(" THIS LINE AND THE FOLLOWING TEXT IS ")"
"("PRINTED BY THIS PROGRAM")",//,
"(" THE RESULTS WITH EFRK ARE:")",/,
"(" R L K DER.EV. Y[1] Y[2]")"
"(" TIME")",/")");
PHI:=4*ARCTAN(1); BETA[0]:=BETA[1]:=1;
"FOR" R:=1,2,3 "DO" "FOR" L:=1,2,3 "DO"
"BEGIN" "FOR" K:=2 "STEP" 1 "UNTIL" R "DO"
BETA[K]:=BETA[K-1]/K;
"FOR" STEP:=1,.1 "DO"
"BEGIN" PASSES:=K:=0; X:=Y[2]:=0; Y[1]:=1; TIME:=CLOCK;
OUT;
EFRK(X,50.0,1,2,Y,SIGMA,PHI,DIAMETER,DER,K,STEP,R,L,BETA,
R>=3,"-4,OUT);
"END"; OUTPUT(61,"("/")");
"END";
"END"
THIS LINE AND THE FOLLOWING TEXT IS PRINTED BY THIS PROGRAM:
THE RESULTS WITH EFRK ARE:
R L K DER.EV. Y[1] Y[2] TIME
1 1 237 474 +.765 812 555 +.433 689 306 1.395
1 1 501 1002 +.765 847 870 +.433 700 619 3.381
1 2 52 156 +.765 570 874 +.433 615 119 0.465
1 2 501 1503 +.765 848 220 +.433 700 709 4.200
1 3 52 208 +.765 571 278 +.433 615 202 0.531
1 3 500 2000 +.765 848 512 +.433 700 827 4.879
2 1 3317 9951 +.765 878 320 +.433 710 353 21.808
2 1 1050 3150 +.765 878 321 +.433 710 330 7.153
2 2 174 696 +.765 878 335 +.433 710 335 1.385
2 2 501 2004 +.765 878 323 +.433 709 211 4.915
2 3 57 285 +.765 881 339 +.433 817 185 0.642
2 3 501 2505 +.765 878 323 +.433 709 725 5.756
3 1 7010 28040 +.765 878 320 +.433 710 354 55.298
3 1 3255 13020 +.765 878 320 +.433 710 374 25.772
3 2 949 4745 +.765 878 319 +.433 711 893 8.499
3 2 1384 6920 +.765 862 498 +.449 724 830 13.452
3 3 917 5502 +.765 878 018 +.434 105 184 9.143
3 3 1166 6996 +.765 861 696 +.433 705 641 15.512
SOURCE TEXT(S):
"CODE" 33070;