NUMAL Section 5.2.1.1.1.1.F
BEGIN SECTION : 5.2.1.1.1.1.F (February, 1979)
PROCEDURE : MULTISTEP.
AUTHOR: P.W.HEMKER.
INSTITUTE: MATHEMATICAL CENTRE.
RECEIVED: 730515.
BRIEF DESCRIPTION:
MULTISTEP SOLVES AN INITIAL VALUE PROBLEM, FOR A SYSTEM OF
FIRST ORDER ORDINARY DIFFERENTIAL EQUATIONS DY = F(Y).
IN PARTICULAR THIS PROCEDURE IS SUITABLE FOR THE INTEGRATION OF
STIFF DIFFERENTIAL EQUATIONS. IT CAN ALSO BE USED FOR
NON-STIFF PROBLEMS.
KEYWORDS:
INITIAL VALUE PROBLEM,
SYSTEM OF FIRST ORDER ORDINARY DIFFERENTIAL EQUATIONS.
STIFF EQUATIONS,
CALLING SEQUENCE:
THE HEADING OF THE PROCEDURE READS:
"BOOLEAN" "PROCEDURE" MULTISTEP(X,XEND,Y,HMIN,HMAX,YMAX,EPS,
FIRST,SAVE,DERIV,AVAILABLE,JACOBIAN,STIFF,N,OUT);
"VALUE" HMIN,HMAX,EPS,XEND,N,STIFF;
"BOOLEAN" FIRST,AVAILABLE,STIFF;
"INTEGER" N;
"REAL" X,XEND,HMIN,HMAX,EPS;
"ARRAY" Y,YMAX,SAVE,JACOBIAN;
"PROCEDURE" DERIV,OUT;
"CODE" 33080;
MULTISTEP DELIVERS THE FOLLOWING BOOLEAN VALUE:
IF DIFFICULTIES ARE ENCOUNTERED DURING THE INTEGRATION
(I.E. SAVE[-1] ^= 0 "OR" SAVE[-2] ^= 0 ) MULTISTEP IS SET
TO "FALSE", OTHERWISE MULTISTEP IS SET "TRUE".
THE MEANING OF THE FORMAL PARAMETERS IS:
X: <VARIABLE>;
THE INDEPENDENT VARIABLE X.
CAN BE USED IN DERIV, AVAILABLE ETC.;
ENTRY: THE INITIAL VALUE X0;
EXIT : THE FINAL VALUE 'XEND';
XEND: <ARITHMETIC EXPRESSION>;
THE FINAL VALUE OF X ( XEND >= X );
Y: <ARRAY IDENTIFIER>;
"ARRAY" Y[1:6*N];
THE DEPENDENT VARIABLE;
ENTRY Y[1:N] : THE INITIAL VALUES OF THE SOLUTION OF THE
SYSTEM OF DIFFERENTIAL EQUATIONS AT X = X0;
EXIT Y[1:N] : THE FINAL VALUES OF THE SOLUTION AT X = XEND;
HMIN,HMAX: <ARITHMETIC EXPRESSION>;
ENTRY : MINIMAL RESP. MAXIMAL STEPLENGTH ALLOWED;
YMAX: <ARRAY IDENTIFIER>;
"ARRAY" YMAX[1:N];
ENTRY: THE ABSOLUTE LOCAL ERROR BOUND DIVIDED BY EPS
EXIT : YMAX[1] GIVES THE MAXIMAL VALUE OF THE ENTRY VALUE
OF YMAX[I] AND THE VALUES OF ABS(Y[I]) DURING
INTEGRATION;
EPS: <ARITHMETIC EXPRESSION>;
THE RELATIVE LOCAL ERROR BOUND;
FIRST: <IDENTIFIER>;
IF FIRST = "TRUE" THEN THE PROCEDURE STARTS THE INTEGRATION
WITH A FIRST ORDER ADAMS METHOD AND A STEPLENGTH EQUAL
TO HMIN. UPON COMPLETION OF A CALL FIRST:= "FALSE";
IF FIRST = "FALSE" THEN THE PROCEDURE CONTINUES
INTEGRATION;
SAVE: <ARRAY IDENTIFIER>;
"ARRAY" SAVE[-38:6*N];
IN THIS ARRAY THE PROCEDURE STORES INFORMATION WHICH CAN BE
USED IN A CONTINUING CALL WITH FIRST="FALSE";
BESIDES THE FOLLOWING MESSAGES ARE DELIVERED:
SAVE[ 0]=0 : AN ADAMS METHOD HAS BEEN USED;
1 : THE PROCEDURE SWITCHED TO GEARS METHOD;
SAVE[-1]=0 : NO ERROR MESSAGE;
1 : WITH THE HMIN SPECIFIED THE PROCEDURE CANNOT
HANDLE THE NONLINEARITY (DECREASE HMIN! );
SAVE[-2] NUMBER OF TIMES THAT THE REQUESTED LOCAL ERROR
BOUND WAS EXCEEDED;
SAVE[-3] IF SAVE[-2] IS NONZERO THEN SAVE[-3] GIVES AN
ESTIMATE OF THE MAXIMAL LOCAL ERROR BOUND,
OTHERWISE SAVE[-3]=0;
DERIV: <PROCEDURE IDENTIFIER>;
THE HEADING OF THIS PROCEDURE READS:
"PROCEDURE" DERIV(DF); "ARRAY" DF;
THIS PROCEDURE SHOULD DELIVER DY[I]/DX IN DF[I];
AVAILABLE: <BOOLEAN EXPRESSION>;
IF AN ANALYTICAL EXPRESSION OF THE JACOBIAN MATRIX IS NOT
AVAILABLE THIS EXPRESSION IS SET TO "FALSE";
OTHERWISE THIS EXPRESSION IS SET TO "TRUE" AND THE
EVALUATION OF THIS BOOLEAN EXPRESSION MUST EFFECT THE
FOLLOWING SIDE-EFFECT:
THE ENTRIES OF THE JACOBIAN MATRIX D(DY[I]/DX)/DY[J] ARE
DELIVERED IN THE ARRAY ELEMENTS JACOBIAN[I,J];
JACOBIAN: <ARRAY IDENTIFIER>;
"ARRAY" JACOBIAN[1:N,1:N];
AT EACH EVALUATION OF THE BOOLEAN EXPRESSION AVAILABLE WITH
THE RESULT AVAILABLE:="TRUE", THE JACOBIAN MATRIX HAS TO BE
ASSIGNED TO THIS ARRAY (SEE THE EXAMPLE OF USE);
STIFF: <BOOLEAN EXPRESSION>;
IF STIFF = "TRUE" THE PROCEDURE SKIPS AN ATTEMPT TO SOLVE
THE PROBLEM WITH ADAMS-BASHFORTH- OR ADAMS-MOULTON
METHODS, DIRECTLY USING GEARS METHOD;
N: <ARITHMETIC EXPRESSION>;
THE NUMBER OF EQUATIONS;
OUT: <PROCEDURE IDENTIFIER>;
THE HEADING OF THIS PROCEDURE READS:
"PROCEDURE" OUT(H,K); "VALUE" H,K; "REAL" H; "INTEGER" K;
AT THE END OF EACH ACCEPTED STEP OF THE INTEGRATION PROCESS
THIS PROCEDURE IS CALLED. THE LAST STEPLENGTH USED (H) AND
THE ORDER OF THE METHOD (K) ARE DELIVERED.
AT EACH CALL OF THE PROCEDURE OUT, THE CURRENT VALUES OF
THE INDEPENDENT VARIABLE (X) AND OF THE SOLUTION (Y[I](X) )
ARE AVAILABLE FOR USE. MOREOVER, IN THE NEIGHBOURHOOD OF
THE CURRENT VALUE OF X, ANY VALUE OF Y[I](X SPECIFIED) CAN
BE COMPUTED BY MEANS OF THE FOLLOWING INTERPOLATION FORMULA
Y[I](X SPECIFIED) =
SUM(J,0,K, Y[I+J*N] * ((X SPECIFIED - X)/H) ** J ).
PROCEDURES USED:
MATVEC = CP34011 ,
DEC = CP34300 ,
SOL = CP34051 .
REQUIRED CENTRAL MEMORY: CIRCA N * ( 2 * N + 5 ) MEMORY PLACES.
METHOD AND PERFORMANCE :
MULTISTEP IS BASED ON TWO LINEAR MULTISTEP METHODS. FOR STIFF
PROBLEMS IT USES THE BACKWARD DIFFERENTIATION METHODS, FOR
FOR NON-STIFF PROBLEMS THE ADAMS-BASHFORTH-MOULTON METHODS.
MULTISTEP IS PROVIDED WITH ORDER, STEPSIZE AND ERROR CONTROL.
REFERENCES:
[1].P.W.HEMKER.
AN ALGOL 60 PROCEDURE FOR THE SOLUTION OF STIFF DIFFERENTIAL
EQUATIONS.
MATH. CENTRE, AMSTERDAM. REPORT MR 128/71;
EXAMPLE OF USE:
THE SOLUTION AT X=1 AND AT X=10 OF THE DIFFERENTIAL EQUATIONS:
DY[1]/DX = 0.04 * (1-Y[1]-Y[2]) - Y[1] * ("4*Y[2] + 3"7*Y[1])
DY[2]/DX = 3"7 Y[1]**2
WITH THE INITIAL CONDITIONS AT X = 0 :
Y[1] = 0 AND Y[2] = 0
MAY BE OBTAINED BY THE FOLLOWING PROGRAM:
"BEGIN"
"BOOLEAN" FIRST;
"INTEGER" I,J,CF,CJ,CA;
"REAL" X,XEND,HMIN,EPS,R;
"ARRAY" Y[1:12],YMAX[1:2],D[-40:12],JAC[1:2,1:2];
"PROCEDURE" DER (F);"ARRAY" F;
"BEGIN" "REAL" R; CF:=CF+1;
F[2]:= R:= 3"7*Y[1]*Y[1];
F[1]:= 0.04*(1-Y[1]-Y[2]) - "4*Y[1]*Y[2] - R;
"END" F;
"BOOLEAN" "PROCEDURE" AVAIL;
"BEGIN" "REAL" R; CJ:= CJ+1;
AVAIL:= "TRUE";
JAC[2,1]:= R:= 6"7*Y[1];
JAC[1,1]:= -0.04 - "4*Y[2] - R;
JAC[1,2]:= -0.04 - "4*Y[1];
JAC[2,2]:= 0
"END" JAC AVAIL;
"PROCEDURE" OUT(H,K);
"VALUE" H, K; "REAL" H; "INTEGER" K; CA:= CA+1;
LABEL:
OUTPUT(61,"("/,"("HMIN,EPS?")",/")");
INREAL(70,HMIN); INREAL(70,EPS);
"IF" HMIN<0 "THEN" "GOTO" ESCAPE;
FIRST:= "TRUE"; CA:=CF:=CJ:=0;
X:=0; Y[1]:= Y[2]:= 0;
YMAX[1]:= 0.0001; YMAX[2]:= 1;
"FOR" XEND:= 1, 10 "DO"
"BEGIN"
MULTISTEP(X,XEND,Y,HMIN,5,YMAX,EPS,FIRST,D,DER,AVAIL,
JAC,"TRUE",2,OUT);
OUTPUT(61,"("3(5ZD,2B),2(+.13D"+2D,2B),/")",
CA,CF,CJ,Y[1],Y[2]);
"END";
"GOTO" LABEL;
ESCAPE:
"END"
IT DELIVERS WITH HMIN = "-10 AND EPS = "-9:
240 648 2 +.3074626[602000]"-04 +.3350951[493111]"-01
315 902 3 +.16233909[62091]"-04 +.15861383[92015]"+00
(NON-SIGNIFICANT DIGITS ARE PLACED BETWEEN [ ] ).
SOURCE TEXT(S):
"CODE" 33080;