NUMAL Section 5.2.1.1.1.2.B
BEGIN SECTION : 5.2.1.1.1.2.B (August, 1974)
AUTHOR: K.DEKKER.
INSTITUTE: MATHEMATICAL CENTRE.
RECEIVED: 1973/07/31.
BRIEF DESCRIPTION:
EFERK SOLVES INITIAL VALUE PROBLEMS , GIVEN AS AN AUTONOMOUS SYSTEM
OF FIRST ORDER DIFFERENTIAL EQUATIONS, BY MEANS OF AN EXPONENTIALLY
FITTED, EXPLICIT RUNGE KUTTA METHOD OF THIRD ORDER, WHICH INVOLVES
THE USE OF THE JACOBIAN MATRIX. AUTOMATIC STEP CONTROL IS PROVIDED.
IN PARTICULAR THIS METHOD IS SUITABLE FOR THE INTEGRATION OF STIFF
DIFFERENTIAL EQUATIONS.
KEYWORDS:
DIFFERENTIAL EQUATIONS,
INITIAL VALUE PROBLEMS,
STIFF EQUATIONS,
EXPONENTIAL FITTING,
EXPLICIT RUNGE KUTTA METHODS.
CALLING SEQUENCE:
THE HEADING OF THE PROCEDURE EFERK READS:
"PROCEDURE" EFERK(X,XE,M,Y,SIGMA,PHI,DERIVATIVE,J,JACOBIAN,
K,L,AUT,AETA,RETA,HMIN,HMAX,LINEAR,OUTPUT);
"VALUE" L;
"INTEGER" M,K,L;
"REAL" X,XE,SIGMA,PHI,AETA,RETA,HMIN,HMAX;
"ARRAY" Y,J;
"BOOLEAN" AUT,LINEAR;
"PROCEDURE" DERIVATIVE,JACOBIAN,OUTPUT;
"CODE" 33120;
THE MEANING OF THE FORMAL PARAMETERS IS:
X: <VARIABLE>;
THE INDEPENDENT VARIABLE X;
CAN BE USED IN DERIVATIVE, JACOBIAN, OUTPUT, ETC.;
ENTRY: THE INITIAL VALUE X0;
EXIT : THE FINAL VALUE XE;
XE: <ARITHMETIC EXPRESSION>;
THE FINAL VALUE OF X (XE>=X);
M: <ARITHMETIC EXPRESSION>;
THE NUMBER OF EQUATIONS;
Y: <ARRAY IDENTIFIER>;
"REAL" "ARRAY" Y[1:M];
THE DEPENDENT VARIABLE;
ENTRY: THE INITIAL VALUES OF THE SYSTEM OF DIFFERENTIAL
EQUATIONS: Y[I] AT X=X0;
EXIT : THE FINAL VALUES OF THE SOLUTION: Y[I] AT X=XE;
SIGMA: <ARITHMETIC EXPRESSION>;
THE MODULUS OF THE POINT AT WHICH EXPONENTIAL FITTING IS
DESIRED, FOR EXAMPLE THE LARGEST NEGATIVE EIGENVALUE OF THE
JACOBIAN MATRIX OF THE SYSTEM OF DIFFERENTIAL EQUATIONS;
PHI: <ARITHMETIC EXPRESSION>;
THE ARGUMENT OF THE COMPLEX POINT AT WHICH EXPONENTIAL
FITTING IS DESIRED;
DERIVATIVE: <PROCEDURE IDENTIFIER>;
THE HEADING OF THIS PROCEDURE READS:
"PROCEDURE" DERIVATIVE(Y); "ARRAY" Y;
THIS PROCEDURE SHOULD DELIVER THE RIGHT HAND SIDE OF THE
I-TH DIFFERENTIAL EQUATION AT THE POINT (Y) AS Y[I];
J: <ARRAY IDENTIFIER>;
"REAL" "ARRAY" J[1:M,1:M];
THE JACOBIAN MATRIX OF THE SYSTEM;
THE ARRAY J SHOULD BE UPDATED IN THE PROCEDURE JACOBIAN;
JACOBIAN: <PROCEDURE IDENTIFIER>;
THE HEADING OF THIS PROCEDURE READS:
"PROCEDURE" JACOBIAN(J,Y); "ARRAY" J,Y;
IN THIS PROCEDURE THE JACOBIAN AT THE POINT (Y) HAS TO BE
ASSIGNED TO THE ARRAY J;
K: <VARIABLE>;
COUNTS THE NUMBER OF INTEGRATION STEPS TAKEN;
FOR EXAMPLE, MAY BE USED IN THE EXPRESSION FOR XE;
L: <ARITHMETIC EXPRESSION>;
ENTRY:
IF PHI = 4*ARCTAN(1): THE ORDER OF THE EXPONENTIAL FITTING,
ELSE TWICE THE ORDER OF THE EXPONENTIAL FITTING;
AUT: <BOOLEAN EXPRESSION>;
IF THE SYSTEM HAS BEEN WRITTEN IN AUTONOMOUS FORM BY ADDING
THE EQUATION DY[M]/DX = 1 TO THE SYSTEM , THEN AUT MAY HAVE
THE VALUE "FALSE", ELSE AUT SHOULD HAVE THE VALUE "TRUE";
AETA: <ARITHMETIC EXPRESSION>;
REQUIRED ABSOLUTE PRECISION IN THE INTEGRATION PROCESS;
AETA HAS TO BE POSITIVE;
RETA: <ARITHMETIC EXPRESSION>;
REQUIRED RELATIVE PRECISION IN THE INTEGRATION PROCESS;
RETA HAS TO BE POSITIVE;
HMIN: <ARITHMETIC EXPRESSION>;
THE STEPLENGTH CHOSEN WILL BE AT LEAST EQUAL TO HMIN;
HMAX: <ARITHMETIC EXPRESSION>;
THE STEPLENGTH CHOSEN WILL BE AT MOST EQUAL TO HMAX;
LINEAR: <ARITHMETIC EXPRESSION>;
THE PROCEDURE JACOBIAN IS CALLED ONLY IF LINEAR="FALSE" OR
K=0; SO IF THE SYSTEM IS LINEAR , LINEAR MAY HAVE THE VALUE
"TRUE";
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, K, Y.
DATA AND RESULTS: SEE EXAMPLE OF USE, AND REF[4].
PROCEDURES USED:
VECVEC= CP34010,
MATVEC= CP34011,
DEC = CP34300,
SOL = CP34051.
REQUIRED CENTRAL MEMORY:
EXECUTION FIELD LENGTH: CIRCA 30 + 4 * M + L * (5+L).
RUNNING TIME: DEPENDS STRONGLY ON THE DIFFERENTIAL EQUATION TO SOLVE.
LANGUAGE: ALGOL 60.
METHOD AND PERFORMANCE:
THE PROCEDURE EFERK IS AN EXPONENTIALLY FITTED, SEMI-EXPLICIT RUNGE
KUTTA METHOD OF THIRD ORDER ( SEE REF [1] AND [3] ) . THE ALGORITHM
USES FOR EACH STEP TWO FUNCTION EVALUATIONS AND IF LINEAR = "FALSE"
ONE EVALUATION OF THE JACOBIAN MATRIX.THE STEPSIZE IS DETERMINED BY
AN ESTIMATION OF THE LOCAL TRUNCATION ERROR BASED ON THE RESIDUAL
FUNCTION (SEE REF[3]). INTEGRATION STEPS ARE NOT REJECTED.
REFERENCES:
[1]. P.J.VAN DER HOUWEN.
ONE-STEP METHODS WITH ADAPTIVE STABILITY FUNCTIONS FOR THE
INTEGRATION OF DIFFERENTIAL EQUATIONS.
LECTURES NOTES OF THE CONFERENCE ON NUMERISCHE, INSBESONDERE
APPROXIMATIONSTHEORETISCHE BEHANDLUNG VON FUNKTIONAL-
GLEICHUNGEN.
OBERWOLFACH, DECEMBER, 3 -12, 1972.
[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]. (TO APPEAR).
COLLOQUIUM STIFF DIFFERENTIAL EQUATIONS 3 (DUTCH).
MC SYLLABUS 15.3, (1973) 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
EFERK,AND ILLUSTRATES THE ACCURACIES WHICH MAY BE OBTAINED BY THEM:
"BEGIN"
"INTEGER" K,PASSES,PASJAC;
"REAL" X,SIGMA,PHI,TIME,TOL;
"REAL" "ARRAY" Y[1:2],J[1:2,1:2];
"PROCEDURE" DER(Y); "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" JACOBIAN(J,Y); "ARRAY" J,Y;
"BEGIN" J[1,1]:=Y[2]-1; J[1,2]:=.99+Y[1];
J[2,1]:=1000*(1-Y[2]); J[2,2]:=-1000*(1+Y[1]);
SIGMA:=ABS(J[2,2]+J[1,1]-SQRT((J[2,2]-J[1,1])**2+
4*J[2,1]*J[1,2]))/2;
PASJAC:=PASJAC+1
"END" JACOBIAN;
"PROCEDURE" OUT;
"IF" X=50 "THEN"
OUTPUT(61,"("3(-5ZD),2(4B+.3DB3DB3D),-5ZD.3D,/")",K,PASSES,
PASJAC,Y[1],Y[2],CLOCK-TIME);
OUTPUT(61,"(""(" THIS LINE AND THE FOLLOWING TEXT IS ")"
"("PRINTED BY THIS PROGRAM")",//,
"(" THE RESULTS WITH EFERK -FIRST ORDER FIT- ARE:")",/,
"(" K DER.EV. JAC.EV. Y[1] Y[2]")"
"(" TIME")",/")");
PHI:=4*ARCTAN(1);
"FOR" TOL:=1,"-1,"-2,"-3 "DO"
"BEGIN" PASSES:=PASJAC:=0; X:=Y[2]:=0; Y[1]:=1; TIME:=CLOCK;
EFERK(X,50.0,2,Y,SIGMA,PHI,DER,J,JACOBIAN,K,1,"TRUE",TOL,
TOL,"-6,50.0,"FALSE",OUT);
"END";
"END"
THIS LINE AND THE FOLLOWING TEXT IS PRINTED BY THIS PROGRAM:
THE RESULTS WITH EFERK -FIRST ORDER FIT- ARE:
K DER.EV. JAC.EV. Y[1] Y[2] TIME
93 186 93 +.765 883 211 +.428 752 781 1.170
105 210 105 +.765 878 445 +.433 569 561 1.296
147 294 147 +.765 878 317 +.433 708 489 1.834
266 532 266 +.765 878 320 +.433 710 229 3.297
SOURCE TEXT(S):
"CODE" 33120;