NUMAL Section 5.2.1.1.1.2.A
BEGIN SECTION : 5.2.1.1.1.2.A (August, 1974)
AUTHOR: S.P.N. VAN KAMPEN.
INSTITUTE: MATHEMATICAL CENTRE.
RECEIVED: 730529.
BRIEF DESCRIPTION:
EFSIRK SOLVES AN INITIAL VALUE PROBLEM, GIVEN AS AN AUTONOMOUS
SYSTEM OF FIRST ORDER DIFFERENTIAL EQUATIONS DY/DX = F(Y), BY
MEANS OF AN EXPONENTIALLY FITTED, SEMI-IMPLICIT RUNGE-KUTTA
METHOD; IN PARTICULAR THIS PROCEDURE IS SUITABLE FOR THE
INTEGRATION OF STIFF EQUATIONS.
KEYWORDS:
DIFFERENTIAL EQUATIONS,
INITIAL VALUE PROBLEM,
AUTONOMOUS SYSTEM,
STIFF EQUATIONS,
SEMI-IMPLICIT RUNGE-KUTTA METHOD,
EXPONENTIAL FITTING.
CALLING SEQUENCE:
HEADING:
"PROCEDURE" EFSIRK(X, XE, M, Y, DELTA, DERIVATIVE, JACOBIAN, J,
N, AETA, RETA, HMIN, HMAX, LINEAR, OUTPUT);
"VALUE" M; "INTEGER" M, N;
"REAL" X, XE, DELTA, AETA, RETA, HMIN, HMAX;
"PROCEDURE" DERIVATIVE, JACOBIAN, OUTPUT;
"ARRAY" Y, J;
"BOOLEAN" LINEAR;
"CODE" 33160;
THE MEANING OF THE FORMAL PARAMETERS IS:
X: <VARIABLE>;
THE INDEPENDENT VARIABLE X;
ENTRY: THE INITIAL VALUE X0;
EXIT : THE END VALUE XE;
XE: <ARITHMETIC EXPRESSION>;
THE END VALUE OF X;
M: <ARITHMETIC EXPRESSION>;
THE NUMBER OF DIFFERENTIAL EQUATIONS;
Y: <ARRAY IDENTIFIER>;
"ARRAY" Y[1 : M];
THE DEPENDENT VARIABLE;
DURING THE INTEGRATION PROCESS THE COMPUTED SOLUTION
AT THE POINT X IS ASSIGNED TO THE ARRAY Y;
ENTRY: THE INITIAL VALUES OF THE SOLUTION OF THE SYSTEM;
DELTA: <ARITHMETIC EXPRESSION>;
DELTA DENOTES THE REAL PART OF THE POINT AT WHICH
EXPONENTIAL FITTING IS DESIRED;
ALTERNATIVES:
DELTA = (AN ESTIMATE OF) THE REAL PART OF THE, IN ABSOLUTE
VALUE, LARGEST EIGENVALUE OF THE JACOBIAN MATRIX OF THE
SYSTEM;
DELTA < -10**14, IN ORDER TO OBTAIN ASYMPTOTIC
STABILITY;
DELTA = 0, IN ORDER TO OBTAIN A HIGHER ORDER OF ACCURACY IN
CASE OF LINEAR OR ALMOST LINEAR EQUATIONS;
DERIVATIVE: <PROCEDURE IDENTIFIER>;
"PROCEDURE" DERIVATIVE(A); "ARRAY" A;
WHEN IN EFSIRK DERIVATIVE IS CALLED, A[I] CONTAINS THE
VALUES OF Y[I];
UPON COMPLETION OF A CALL OF DERIVATIVE, THE ARRAY A
SHOULD CONTAIN THE VALUES OF F(Y);
NOTE THAT THE VARIABLE X SHOULD NOT BE USED IN DERIVATIVE,
BECAUSE THE DIFFERENTIAL EQUATION IS SUPPOSED TO BE
AUTONOMOUS;
JACOBIAN: <PROCEDURE IDENTIFIER>;
"PROCEDURE" JACOBIAN(J, Y); "ARRAY" J, Y;
WHEN IN EFSIRK JACOBIAN IS CALLED THE ARRAY Y CONTAINS
THE VALUES OF THE DEPENDENT VARIABLE;
UPON COMPLETION OF A CALL OF JACOBIAN THE ARRAY J SHOULD
CONTAIN THE VALUES OF THE JACOBIAN MATRIX OF F(Y);
J: <ARRAY IDENTIFIER>;
J[1 : M, 1 : M];
J IS AN AUXILLIARY ARRAY WHICH IS USED IN THE PROCEDURE
JACOBIAN;
N: <VARIABLE>;
AN INTEGER WHICH COUNTS THE INTEGRATION STEPS;
AETA, RETA:
<ARITHMETIC EXPRESSION>;
REQUIRED ABSOLUTE AND RELATIVE LOCAL ACCURACY;
HMIN, HMAX:
<ARITHMETIC EXPRESSION>;
MINIMAL AND MAXIMAL STEPSIZE BY WHICH THE INTEGRATION IS
PERFORMED;
LINEAR: <BOOLEAN EXPRESSION>;
IF LINEAR = "TRUE" THE PROCEDURE JACOBIAN WILL ONLY BE
CALLED IF N = 1; THE INTEGRATION WILL THEN BE PERFORMED
WITH A STEPSIZE HMAX; THE CORRESPONDING REDUCTION
OF COMPUTING TIME CAN BE EXPLOITED IN CASE OF LINEAR OR
ALMOST LINEAR EQUATIONS;
OUTPUT: <PROCEDURE IDENTIFIER>;
"PROCEDURE" OUTPUT;
IN OUTPUT ONE MAY PRINT THE VALUES OF E.G. X,
Y[I], J[K, L] AND N.
DATA AND RESULTS: SEE REF[2] AND [3].
PROCEDURES USED:
VECVEC = CP34010,
MATVEC = CP34011,
MATMAT = CP34013,
GSSELM = CP34231,
SOLELM = CP34061.
REQUIRED CENTRAL MEMORY:
EXECUTION FIELD LENGTH: CIRCA M * M + 5 * M.
RUNNING TIME:
DEPENDS STRONGLY ON THE DIFFERENTIAL EQUATION TO BE SOLVED
LANGUAGE: ALGOL 60.
METHOD AND PERFORMANCE:
THE PROCEDURE EFSIRK IS AN EXPONENTIALLY FITTED, A-STABLE,
SEMI-IMPLICIT RUNGE-KUTTA METHOD OF THIRD ORDER (SEE REF[1]
AND [2]). THE ALGORITHM USES FOR EACH STEP TWO FUNCTION EVALUATIONS
AND IF LINEAR = "FALSE" ONE EVALUATION OF THE JACOBIAN MATRIX.
THE STEPSIZE IS NOT DETERMINED BY THE ACCURACY OF THE NUMERICAL
SOLUTION, BUT BY THE AMOUNT BY WHICH THE GIVEN DIFFERENTIAL
EQUATION DIFFERS FROM A LINEAR EQUATION (SEE REF[2]).
THE PROCEDURE DOES NOT REJECT INTEGRATION STEPS.
REFERENCES:
[1].P.J. VAN DER HOUWEN.
ONE-STEP METHODS WITH ADAPTIVE STABILITY FUCNTIONS
FOR THE INTEGRATION OF DIFFERENTIAL EQUATIONS.
LECTURE NOTES OF THE CONFERENCE ON
NUMERISCHE, INSBESONDERE APPROXIMATIONSTHEORETISCHE
BEHANDLUNG VON FUNKTIONALGLEICHUNGEN.
OBERWOLFACH, DECEMBER, 3 - 12, 1972.
[2].SYLLABUS COLLOQUIUM STIFF DIFFERENTIAL EQUATIONS 2 (DUTCH).
MATH.CENTR. SYLLABUS 15.2/73.
[3].SYLLABUS COLLOQUIUM STIFF DIFFERENTIAL EQUATIONS 3 (DUTCH).
MATH.CENTR. SYLLABUS 15.3/73.
TO APPEAR IN 1973.
EXAMPLE OF USE:
WE CONSIDER THE DIFFERENTIAL EQUATION
DY / DX = -EXP(X) * (Y - LN(X)) + 1 / X,
ON THE INTERVAL [0.01, 8], WITH INITIAL VALUE Y(0.01) = LN(0.01)
AND ANALYTICAL SOLUTION Y(X) = LN(X);
FOR THE FIT POINT WE USE THE EIGENVALUE OF THE JACOBIAN MATRIX,
I.E. DELTA = -EXP(X);
"BEGIN"
"PROCEDURE" DER(Y); "ARRAY" Y;
"BEGIN" "REAL" Y2; Y2:= Y[2];
DELTA:= -EXP(Y2); LNX:= LN(Y2);
Y[1]:= (Y[1] - LNX) * DELTA + 1 / Y2;
Y[2]:= 1
"END" DER;
"PROCEDURE" JAC(J, Y); "ARRAY" J, Y;
"BEGIN" "REAL" Y2; Y2:= Y[2];
J[1, 1]:= DELTA;
J[1, 2]:= (Y[1] - LNX - 1 / Y2) * DELTA - 1 / (Y2 * Y2);
J[2, 1]:= J[2, 2]:= 0
"END" JAC;
"PROCEDURE" OUTP;
"IF" X = XE "THEN"
"BEGIN" "REAL" Y1; Y1:= Y[1]; LNX:= LN(X);
OUTPUT(61, "(""("N = ")", 2ZD,
"(" X = ")", +D.D,
"(" Y(X) = ")", +D.5D,
"(" DELTA = ")", +3ZD.2D, /,
"("ABS. ERR. = ")", .2D"+2D,
"(" REL. ERR. = ")", .2D"+2D, //")",
N, X, Y1, DELTA,
ABS(Y1 - LNX), ABS((Y1 - LNX) / LNX));
"IF" X = 0.4 "THEN" XE:= 8
"END" OUTP;
"INTEGER" N;
"REAL" X, XE, DELTA, LNX;
"ARRAY" Y[1 : 2], J[1 : 2, 1 : 2];
XE:= 0.4; X:= 0.01; Y[1]:= LN(0.01); Y[2]:= X;
EFSIRK(X, XE, 2, Y, DELTA, DER, JAC, J,
N, "-2, "-2, 0.005, 1.5, "FALSE", OUTP)
"END"
THIS PROGRAM DELIVERS:
N = 10 X = +0.4 Y(X) = -0.91099 DELTA = -1.44
ABS. ERR. = .53"-02 REL. ERR. = .58"-02
N = 98 X = +8.0 Y(X) = +2.07911 DELTA = -2980.02
ABS. ERR. = .33"-03 REL. ERR. = .16"-03
SOURCE TEXT(S):
"CODE" 33160;