NUMAL Section 5.2.1.1.1.2.E
BEGIN SECTION : 5.2.1.1.1.2.E (October, 1974)
AUTHOR: J.G. VERWER.
INSTITUTE: MATHEMATICAL CENTRE.
RECEIVED: 740809.
BRIEF DESCRIPTION:
GMS SOLVES AN INITIAL VALUE PROBLEM, GIVEN AS AN AUTONOMOUS
SYSTEM OF FIRST ORDER DIFFERENTIAL EQUATIONS DY/DX = F(Y), BY
MEANS OF A THIRD ORDER GENERALIZED LINEAR MULTISTEP METHOD;
IN PARTICULAR THIS PROCEDURE IS SUITABLE FOR THE INTEGRATION
OF STIFF EQUATIONS.
KEYWORDS:
DIFFERENTIAL EQUATIONS,
INITIAL VALUE PROBLEM,
AUTONOMOUS SYSTEM,
STIFF EQUATIONS,
GENERALIZED LINEAR MULTISTEP METHOD.
CALLING SEQUENCE:
THE HEADING OF THE PROCEDURE READS:
"PROCEDURE" GMS(X, XE, R, Y, H, HMIN, HMAX, DELTA, DERIVATIVE,
JACOBIAN, AETA, RETA, N, JEV, LU, NSJEV,
LINEAR, OUT);
"VALUE" R;
"REAL" X, XE, H, HMIN, HMAX, DELTA, AETA, RETA;
"INTEGER" R, N, JEV, NSJEV, LU;
"BOOLEAN" LINEAR;
"ARRAY" Y;
"PROCEDURE" DERIVATIVE, JACOBIAN, OUT;
"CODE" 33191;
GMS INTEGRATES THE SYSTEM OF DIFFERENTIAL EQUATIONS DY/DX = F(Y)
FROM X = X0 TO X = XE;
THE MEANING OF THE FORMAL PARAMETERS IS:
X: <VARIABLE>;
THE INDEPENDENT VARIABLE X;
ENTRY: THE INITIAL VALUE OF X;
EXIT : THE END VALUE OF X;
XE: <ARITHMETIC EXPRESSION>;
ENTRY: THE END VALUE OF X;
R: <ARITHMETIC EXPRESSION>;
ENTRY: THE NUMBER OF DIFFERENTIAL EQUATIONS;
Y: <ARRAY IDENTIFIER>;
"ARRAY" Y[1:R];
THE DEPENDENT VARIABLE;
ENTRY: THE INITIAL VALUE OF Y;
EXIT : THE SOLUTION Y AT THE POINT X AFTER EACH
INTEGRATION STEP;
H: <ARITHMETIC EXPRESSION>;
ENTRY: THE STEPLENGTH WHEN THE INTEGRATION HAS TO BE
PERFORMED WITHOUT THE STEPSIZE MECHANISM, OTHER-
WISE THE INITIAL STEPLENGTH (SEE THE PARAMETERS
HMIN AND HMAX);
HMIN, HMAX: <ARITHMETIC EXPRESSION>;
ENTRY: MINIMAL AND MAXIMAL STEPLENGTH BY WHICH THE INTE-
GRATION IS ALLOWED TO BE PERFORMED;
BY PUTTING HMIN = HMAX THE STEPSIZE MECHANISM IS
ELIMINATED; IN THIS CASE THE GIVEN VALUES FOR HMIN AND
HMAX ARE IRRELEVANT, WHILE THE INTEGRATION IS PERFORMED
WITH THE STEPLENGTH GIVEN BY H;
DELTA: <ARITHMETIC EXPRESSION>;
ENTRY: THE REAL PART OF THE POINT AT WHICH EXPONENTIAL
FITTING IS DESIRED;
(SEE METHOD AND PERFORMANCE);
ALTERNATIVES:
DELTA = (AN ESTIMATE OF) THE REAL PART OF THE LARGEST
EIGENVALUE IN MODULUS OF THE JACOBIAN MATRIX OF THE
SYSTEM;
DELTA <= -10**15, IN ORDER TO OBTAIN ASYMPTOTIC STABILITY;
DELTA = 0, IN ORDER TO OBTAIN A HIGHER ORDER OF ACCURACY
IN CASE OF LINEAR EQUATIONS;
DERIVATIVE: <PROCEDURE IDENTIFIER>;
"PROCEDURE" DERIVATIVE(Y); "ARRAY" Y;
<REPLACEMENT OF THE I-TH COMPONENT OF THE SOLUTION Y BY
THE I-TH COMPONENT OF THE DERIVATIVE F(Y), I = 1,..., R>;
JACOBIAN: <PROCEDURE IDENTIFIER>;
"PROCEDURE" JACOBIAN(J, Y); "ARRAY" J, Y;
WHEN IN GMS 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);
AETA, RETA: <ARITHMETIC EXPRESSION>;
ENTRY: MEASURE OF THE ABSOLUTE AND RELATIVE LOCAL
ACCURACY REQUIRED;
THESE VALUES ARE IRRELEVANT WHEN THE INTEGRATION IS PER-
FORMED WITHOUT THE STEPSIZE MECHANISM;
N: <VARIABLE>;
EXIT : THE NUMBER OF INTEGRATION STEPS;
JEV: <VARIABLE>;
EXIT: THE NUMBER OF JACOBIAN EVALUATIONS;
LU: <VARIABLE>;
EXIT: THE NUMBER OF LU-DECOMPOSITIONS;
NSJEV: <VARIABLE>;
ENTRY: NUMBER OF INTEGRATION STEPS PER
JACOBIAN EVALUATION;
THE VALUE OF NSJEV IS RELEVANT ONLY WHEN THE INTEGRATION
IS PERFORMED WITHOUT THE STEPSIZE MECHANISM AND THE
SYSTEM TO BE SOLVED IS NON-LINEAR;
LINEAR: <BOOLEAN EXPRESSION>;
ENTRY: TRUE WHEN THE SYSTEM TO BE INTEGRATED IS LINEAR,
OTHERWISE FALSE;
IF LINEAR IS TRUE THE STEPSIZE MECHANISM IS AUTOMATICALLY
ELIMINATED;
OUT: <PROCEDURE IDENTIFIER>;
"PROCEDURE" OUT;
<BY MEANS OF OUT ONE MAY PRINT THE VALUES OF THE RELEVANT
PARAMETERS OCCURRING IN THE PARAMETERLIST; OUT IS CALLED
AFTER EACH INTEGRATION STEP>;
DATA AND RESULTS: SEE REF[2].
PROCEDURES USED:
VECVEC = CP34010,
MATVEC = CP34011,
MATMAT = CP34013,
ELMROW = CP34024,
ELMVEC = CP34020,
DUPVEC = CP31030,
GSSELM = CP34231,
SOLELM = CP34061,
COLCST = CP31131,
MULVEC = CP31020.
REQUIRED CENTRAL MEMORY:
EXECUTION FIELD LENGTH: 8 * R + 3 * R * R;
RUNNING TIME:
DEPENDS STRONGLY ON THE DIFFERENTIAL EQUATION TO BE SOLVED.
LANGUAGE: ALGOL 60.
METHOD AND PERFORMANCE:
THE PROCEDURE GMS DESCRIBES AN IMPLEMENTATION OF A THIRD ORDER
THREE-STEP GENERALIZED LINEAR MULTISTEP METHOD WITH QUASI-ZERO
PARASITIC ROOTS AND QUASI-ADAPTIVE STABILITY FUNCTION. IN PARTI-
CULAR THE ALGORITHM IS DEVELOPED FOR THE INTEGRATION OF STIFF
SYSTEMS OF ORDINARY DIFFERENTIAL EQUATIONS. THE PROCEDURE SUPPLIES
THE ADDITIONAL STARTING VALUES AND PERFORMS A STEPSIZE CONTROL
WHICH IS BASED ON THE NON-LINEARITY OF THE DIFFERENTIAL EQUATION.
BY THIS CONTROL THE JACOBIAN MATRIX IS INCIDENTALLY EVALUATED. IT
IS POSSIBLE TO ELIMINATE THE STEPSIZE CONTROL. THEN, ONE HAS
TO GIVE THE NUMBER OF INTEGRATION STEPS PER JACOBIAN EVALUATION.
FOR LINEAR EQUATIONS THE STEPSIZE CONTROL IS AUTOMATICALLY ELIMIN-
ATED, WHILE THE PROCEDURE PERFORMS ONE EVALUATION OF THE JACOBIAN.
MOREOVER, IN THIS CASE THE THREE-STEP SCHEME IS REDUCED TO A ONE-
STEP SCHEME. THE PROCEDURE USES ONE FUNCTION EVALUATION PER INTE-
GRATION STEP AND IT DOES NOT REJECT INTEGRATION STEPS. EACH CHANGE
IN THE STEPLENGTH OR EACH REEVALUATION OF THE JACOBIAN COSTS ONE
LU-DECOMPOSITION. IT IS POSSIBLE TO FIT EXPONENTIALLY, THIS FIT-
TING IS EQUIVALENT TO FITTING IN THE SENSE OF LINIGER AND
WILLOUGHBY, ONLY WHEN THE JACOBIAN MATRIX IS EVALUATED AT EACH IN-
TEGRATION STEP. WHEN THE SYSTEM TO BE INTEGRATED IS NON-LINEAR
AND THE JACOBIAN MATRIX IS NOT EVALUATED AT EACH INTEGRATION STEP,
IT IS RECOMMENDED TO FIT AT INFINITY (DELTA <= -10**15).
DETAILS ARE GIVEN IN REFERENCE 2.
REFERENCES:
[1] HOUWEN, P. J. VAN DER AND VERWER, J. G.,
GENERALIZED LINEAR MULTISTEP METHODS 1, DEVELOPMENT OF ALGO-
RITHMS WITH ZERO-PARASITIC ROOTS,
REPORT NW 10/74, MATHEMATISCH CENTRUM, AMSTERDAM 1974.
[2] VERWER, J. G.,
GENERALIZED LINEAR MULTISTEP METHODS 2, NUMERICAL
APPLICATIONS, REPORT NW 12/74, MATHEMATISCH CENTRUM,
AMSTERDAM, 1974.
EXAMPLE OF USE:
WE CONSIDER THE DIFFERENTIAL EQUATION
DY1/DX = -1000 * Y1 * (Y1 + Y2 -1.999987),
DY2/DX = -2500 * Y2 * (Y1 + Y2 - 2),
ON THE INTERVAL [0,50], WITH INITIAL VALUE Y1(0) = Y2(0) = 1.
THE REFERENCE SOLUTION AT X = 50 IS GIVEN BY:
Y1(50) = .5976546988,
Y2(50) = 1.4023434075.
"BEGIN"
"PROCEDURE" DER(Y); "ARRAY" Y;
"BEGIN" "REAL" Y1, Y2;
Y1:= Y[1]; Y2:= Y[2];
Y[1]:= -1000 * Y1 * (Y1 + Y2 - 1.999987);
Y[2]:= -2500 * Y2 * (Y1 + Y2 - 2)
"END" DER;
"PROCEDURE" JAC(J, Y); "ARRAY" J, Y;
"BEGIN" "REAL" Y1, Y2; Y1:= Y[1]; Y2:= Y[2];
J[1,1]:= 1999.987 - 1000 * (2 * Y1 + Y2);
J[1,2]:= -1000 * Y1;
J[2,1]:= -2500 * Y2;
J[2,2]:= 2500 * (2 - Y1 - 2 * Y2)
"END" JAC;
"PROCEDURE" OUTP;
"IF" X = 50 "THEN"
"BEGIN" "REAL" YE1, YE2;
YE1:= .5976546988; YE2:= 1.4023434075;
OUTPUT(61, "("
"("X = ")", 2D2B,
"("N = ")", 3ZD2B,
"("JEV = ")", 3ZD2B,
"("LU = ")", 3ZD, 2/,
"("Y1 = ")", +.13D"+2D2B,
"("REL. ERR. = ")", .2D"+2D, /,
"("Y2 = ")", +.13D"+2D2B,
"("REL. ERR. = ")", .2D"+2D")",
X, N, JEV, LU, Y[1], ABS((Y[1] - YE1) / YE1),
Y[2], ABS((Y[2] - YE2) / YE2))
"END" OUTP;
"INTEGER" N, JEV, LU; "REAL" X;
"ARRAY" Y[1:2]; X:= 0.0; Y[1]:= Y[2]:= 1.0;
GMS(X, 50.0, 2, Y, .01, .001, .5, -"15,
DER, JAC, "-5, "-5, N, JEV,
LU, 0, "FALSE", OUTP)
"END"
THIS PROGRAM DELIVERS:
X = 50 N = 109 JEV = 3 LU = 12
Y1 = +.5976547958004"+00 REL. ERR. = .16"-06
Y2 = +.1402343310813"+01 REL. ERR. = .69"-07
SOURCE TEXT:
"CODE" 33191;