NUMAL Section 5.2.1.1.1.2.F
BEGIN SECTION : 5.2.1.1.1.2.F (October, 1975)
AUTHOR: B.LINDBERG.
CONTRIBUTOR: K.DEKKER.
INSTITUTE: MATHEMATICAL CENTRE.
RECEIVED: 741101.
BRIEF DESCRIPTION:
IMPEX SOLVES AN INITIAL VALUE PROBLEM,GIVEN AS AN AUTONOMOUS SYSTEM
OF FIRST ORDER DIFFERENTIAL EQUATIONS , BY MEANS OF THE IMPLICIT
MIDPOINT RULE WITH SMOOTHING AND EXTRAPOLATION.
AUTOMATIC STEPSIZE CONTROL IS PROVIDED.
IN PARTICULAR THIS METHOD IS SUITABLE FOR THE INTEGRATION OF STIFF
DIFFERENTIAL EQUATIONS.
KEYWORDS:
DIFFERENTIAL EQUATIONS,
INITIAL VALUE PROBLEMS,
STIFF EQUATIONS,
IMPLICIT MIDPOINT RULE,
SMOOTHING,
EXTRAPOLATION.
CALLING SEQUENCE:
THE HEADING OF THE PROCEDURE IMPEX READS:
"PROCEDURE" IMPEX (N, T0, TEND, Y0, DERIV, AVAILABLE, H0, HMAX,
PRESCH, EPS, WEIGHTS, UPDATE, FAIL, CONTROL);
"VALUE" N;
"INTEGER" N;
"REAL" T0,TEND,H0,HMAX,EPS;
"BOOLEAN" PRESCH,FAIL;
"ARRAY" Y0,WEIGHTS;
"BOOLEAN" "PROCEDURE" AVAILABLE;
"PROCEDURE" DERIV,UPDATE,CONTROL;
"CODE" 33135;
IMPEX: INTEGRATES THE SYSTEM OF DIFFERENTIAL EQUATIONS , GIVEN AS
DY/DT=F(T,Y), FROM T0 UNTIL TEND.
THE MEANING OF THE FORMAL PARAMETERS IS:
N: <ARITHMETIC EXPRESSION>;
THE NUMBER OF EQUATIONS;
T0: <ARITHMETIC EXPRESSION>;
THE INITIAL VALUE OF THE INDEPENDENT VARIABLE;
TEND: <ARITHMETIC EXPRESSION>;
THE FINAL VALUE OF THE INDEPENDENT VARIABLE;
Y0: <ARRAY IDENTIFIER>;
"REAL" "ARRAY" Y0[1:N];
ENTRY: THE INITIAL VALUES OF THE SYSTEM OF DIFFERENTIAL
EQUATIONS: Y0[I] AT T=T0;
DERIV: <PROCEDURE" IDENTIFIER>;
THE HEADING OF THIS PROCEDURE READS:
"PROCEDURE" DERIV(T,Y,F,N);
"INTEGER" N; "REAL" T; "ARRAY" Y,F;
THIS PROCEDURE SHOULD DELIVER THE VALUE OF F(T,Y) IN THE
ARRAY F[1:N];
AVAILABLE: <PROCEDURE IDENTIFIER>;
THE HEADING OF THIS PROCEDURE READS:
"BOOLEAN" "PROCEDURE" AVAILABLE(T,Y,A,N);
"INTEGER" N; "REAL" T; "ARRAY" Y,A;
IF AN ANALYTIC EXPRESSION OF THE JACOBIAN MATRIX AT THE
POINT (T,Y) IS NOT AVAILABLE, THIS PROCEDURE SHOULD DELIVER
THE VALUE "FALSE";
OTHERWISE THE PROCEDURE SHOULD DELIVER THE VALUE "TRUE",
AND THE JACOBIAN MATRIX SHOULD BE ASSIGNED TO THE ARRAY
A[1:N,1:N];
H0: <ARITHMETIC EXPRESSION>;
THE INITIAL STEPSIZE;
HMAX: <ARITHMETIC EXPRESSION>;
MAXIMAL STEPSIZE BY WHICH THE INTEGRATION IS PERFORMED;
PRESCH: <BOOLEAN EXPRESSION>;
INDICATOR FOR CHOICE OF STEPSIZE;
THE STEPSIZE IS AUTOMATICALLY CONTROLLED IF PRESCH="FALSE";
OTHERWISE THE STEPSIZE HAS TO BE PRESCRIBED , EITHER ONLY
INITIALLY OR ALSO BY THE PROCEDURE CONTROL;
EPS: <ARITHMETIC EXPRESSION>;
BOUND FOR THE ESTIMATE OF THE LOCAL ERROR;
WEIGHTS: <ARRAY IDENTIFIER>;
"REAL" "ARRAY" WEIGHTS[1:N];
WEIGHTS FOR THE COMPUTATION OF THE WEIGHTED EUCLIDEAN NORM
OF THE ERRORS;
ENTRY: INITIAL WEIGHTS;
NOTE THAT THE CHOICE WEIGHTS[I] = 1 IMPLIES AN ESTIMATION
OF THE ABSOLUTE ERROR , WHEREAS WEIGTHS[I] = Y[I] DEFINES A
RELATIVE ERROR;
UPDATE: <PROCEDURE IDENTIFIER>;
THE HEADING OF THIS PROCEDURE READS:
"PROCEDURE" UPDATE(WEIGHTS,Y2,N);
"INTEGER" N; "ARRAY" WEIGHTS,Y2;
THIS PROCEDURE MAY CHANGE THE ARRAY WEIGHTS , ACCORDING TO
THE VALUE OF AN APPROXIMATION FOR Y(T) , GIVEN IN THE ARRAY
Y2[1:N];
FAIL: <BOOLEAN EXPRESSION>;
EXIT :
IF THE PROCEDURE FAILS TO SOLVE THE SYSTEM OF EQUATIONS,
FAIL WILL HAVE THE VALUE "TRUE" UPON EXIT;
THIS MAY OCCUR UPON DIVERGENCE OF THE ITERATION PROCESS,
USED IN THE MIDPOINT RULE , WHILE INTEGRATION IS PERFORMED
WITH A USER DEFINED PRESCRIBED STEPSIZE;
CONTROL: <PROCEDURE IDENTIFIER>;
THE HEADING OF THIS PROCEDURE READS:
"PROCEDURE" CONTROL(TPRINT,T,H,HNEW,Y,ERROR,N);
"INTEGER" N; "REAL" TPRINT,T,H,HNEW; "ARRAY" Y,ERROR;
CONTROL IS CALLED ON ENTRY OF IMPEX, AND FURTHER AS SOON AS
THE INEQUALITY TPRINT <= T HOLDS;
DURING A CALL OF CONTROL PRINTING OF RESULTS AND
CHANGE OF STEPSIZE (IF PRESCH = "TRUE") IS THEN POSSIBLE;
THE MEANING OF THE FORMAL PARAMETERS IS:
TPRINT: <VARIABLE>;
ENTRY: THE VALUE OF THE INDEPENDENT VARIABLE AT
WHICH A CALL OF CONTROL WAS DESIRED;
EXIT: A NEW VALUE (TPRINT>T) AT WHICH A CALL OF
CONTROL IS DESIRED;
T: <VARIABLE>;
THE ACTUAL VALUE OF THE INDEPENDENT VARIABLE, UP TO
WHICH INTEGRATION HAS BEEN PERFORMED;
H: <VARIABLE>;
HALVE THE ACTUAL STEPSIZE;
HNEW: <VARIABLE>;
THE NEW STEPSIZE;
IF PRESCH="TRUE", THEN THE USER MAY PRESCRIBE A NEW
STEPSIZE BY CHANGING HNEW;
Y: <ARRAY IDENTIFIER>;
"REAL" "ARRAY" Y[1:5,1:N];
THE VALUE OF THE DEPENDENT VARIABLE AND ITS FIRST
FOUR DIVIDED DIFFERENCES AT THE POINT T ARE GIVEN
IN THIS ARRAY;
ERROR: <ARRAY IDENTIFIER>;
"REAL" "ARRAY" ERROR[1:3];
THE ELEMENTS OF THIS ARRAY CONTAIN THE FOLLOWING
ERRORS:
ERROR[1]: THE LOCAL ERROR;
ERROR[2]: THE GLOBAL ERROR OF SECOND ORDER IN H;
ERROR[3]: THE GLOBAL ERROR OF FOURTH ORDER IN H;
N: <VARIABLE>;
THE NUMBER OF EQUATIONS;
EXAMPLE OF USE: SEE EXAMPLE OF USE OF THE PROCEDURE IMPEX;
DATA AND RESULTS:
FOR DATA, SEE REF[1].
THE RESULTS OF THE INTEGRATION ARE ATTAINABLE THROUGH THE PROCEDURE
CONTROL , WHICH IS CALLED AT SPECIFIED, USER DEFINED, VALUES OF THE
INDEPENDENT VARIABLE . IN PARTICULAR , THE VALUES OF THE DEPENDENT
VARIABLE AT THE ENDPOINT OF INTEGRATION ARE OBTAINED BY A CALL OF
CONTROL WITH TPRINT=TEND.
PROCEDURES USED:
INIVEC = CP31010,
INIMAT = CP31011,
MULVEC = CP31020,
MULROW = CP31021,
DUPVEC = CP31030,
DUPROWVEC= CP31032,
DUPMAT = CP31035,
VECVEC = CP34010,
MATVEC = CP34011,
MATMAT = CP34013,
ELMVEC = CP34020,
ELMROW = CP34024,
DEC = CP34300,
SOL = CP34051.
REQUIRED CENTRAL MEMORY:
EXECUTION FIELD LENGTH: CIRCA N * ( 23 + 2 * N ) (DECIMAL).
RUNNING TIME: DEPENDS STRONGLY ON THE DIFFERENTIAL EQUATION TO SOLVE.
LANGUAGE: ALGOL 60.
METHOD AND PERFORMANCE:
THE INTEGRATION METHOD (REF[1]) IS BASED ON THE COMPUTATION OF TWO
INDEPENDENT SOLUTIONS Y(T,H) AND Y(T,H/2) BY THE IMPLICIT MIDPOINT
RULE. PASSIVE SMOOTHING AND PASSIVE EXTRAPOLATION IS PERFORMED TO
OBTAIN STABILITY AND HIGH ACCURACY . THE ALGORITHM USES FOR EACH
STEP AT LEAST THREE FUNCTION EVALUATIONS, AND ON CHANGE OF STEPSIZE
OR AT SLOW CONVERGENCE IN THE ITERATION PROCESS AN APPROXIMATION OF
THE JACOBIAN MATRIX (COMPUTED BY DIVIDED DIFFERENCES OR EXPLICITLY
SPECIFIED BY THE USER) . IF THE COMPUTED LOCAL ERROR EXCEEDS THE
TOLERANCE , THE LAST STEP IS REJECTED. MOREOVER , TWO GLOBAL ERRORS
ARE COMPUTED.
REFERENCES:
[1]. B.LINDBERG.
IMPEX 2 , A PROCEDURE FOR THE SOLUTION OF SYSTEMS OF STIFF
DIFFERENTIAL EQUATIONS.
ROYAL INSTITUTE OF TECHNOLOGY,STOCKHOLM. TRITA-NA-7303 (1973).
[2]. (TO APPEAR)
COLLOQUIUM STIFF DIFFERENTIAL EQUATIONS 3 (DUTCH).
M.C. SYLLABUS 15.3 (1974), MATHEMATICAL CENTRE.
EXAMPLE OF USE:
CONSIDER THE AUTONOMOUS SYSTEM OF DIFFERENTIAL EQUATIONS:
DY[1]/DX = .2 * ( Y[2] - Y[1] ),
DY[2]/DX = 10 * Y[1] - ( 60 - Y[3]/8 ) * Y[2] + Y[3]/8,
DY[3]/DX = 1,
WITH INITIAL CONDITIONS AT X=0 : Y[1]=Y[2]=Y[3]=0 (SEE REF[2]).
THE SOLUTION AT SEVERAL POINTS IN THE INTERVAL [0, 400] MAY BE
OBTAINED BY THE FOLLOWING PROGRAM:
(THE SOLUTION AT X=400 IS: Y[1]=22.24222011, Y[2]=27.11071335)
"BEGIN" "INTEGER" N,NFE,NJE,POINT;
"REAL" T,TEND,EPS,HMAX,L,H2,TIME;
"ARRAY" Y,SW[1:3],PRINT[1:5];
"BOOLEAN" FAIL;
"PROCEDURE" LIPEST(L,Y,EPS,T,F,N);
"REAL" T,L,EPS; "ARRAY" Y; "INTEGER" N; "PROCEDURE" F;
"BEGIN" "REAL" N1,N2; "INTEGER" I,IT; "ARRAY" F1,F2,Z,X[1:N];
"REAL" "PROCEDURE" NORM(Y); "ARRAY" Y;
NORM:=SQRT(VECVEC(1,N,0,Y,Y));
DUPVEC(1,N,0,Z,Y);
"FOR" I:=1 "STEP" 1 "UNTIL" N "DO"
X[I]:="IF" Y[I]=0 "THEN" EPS "ELSE" (1+EPS)*Y[I];
N1:=NORM(X)*EPS; F(T,X,F1,N);
"FOR" IT:=1 "STEP" 1 "UNTIL" 5 "DO"
"BEGIN" F(T,Z,F2,N);
ELMVEC(1,N,0,F2,F1,-1);
N2:=N1/NORM(F2);
DUPVEC(1,N,0,Z,X); ELMVEC(1,N,0,Z,F2,N2)
"END";
F(T,Z,F2,N);
ELMVEC(1,N,0,F2,F1,-1);
L:=NORM(F2)/N1
"END" LIPEST;
"PROCEDURE" F(T,Y,F1,N); "VALUE" T; "REAL" T; "ARRAY" Y,F1;
"INTEGER" N;
"BEGIN" NFE:=NFE+1;
F1[1]:=0.2*(Y[2]-Y[1]);
F1[2]:=10*Y[1]-(60-.125*Y[3])*Y[2]+.125*Y[3];
F1[3]:=1
"END";
"BOOLEAN" "PROCEDURE" AVAILABLE(T,Y,A,N);
"INTEGER" N; "REAL" T; "ARRAY" Y,A;
"BEGIN" NJE:=NJE+1; AVAILABLE:="TRUE";
A[1,1]:=-.2; A[1,2]:=.2; A[1,3]:=A[3,1]:=A[3,2]:=A[3,3]:=0;
A[2,1]:=10; A[2,2]:=.125*Y[3]-60; A[2,3]:=.125*(1+Y[2])
"END"
"PROCEDURE" UPDATE(SW,R1,N); "INTEGER" N; "ARRAY" SW,R1;
"BEGIN" "REAL" S1,S2; "INTEGER" I;
"FOR" I:=1 "STEP" 1 "UNTIL" N "DO"
"BEGIN" S1:=1/SW[I]; S2:=ABS(R1[I]);
"IF" S1<S2 "THEN" SW[I]:=1/S2
"END"
"END";
"PROCEDURE" CONTROL(TP,T,H,HNEW,Y,ERR,N);
"REAL" TP,T,H,HNEW; "ARRAY" Y,ERR; "INTEGER" N;
"BEGIN" "INTEGER" I;
"ARRAY" C[3:5],X[1:N];
"REAL" S,S2,S3,S4,C1;
NEXT: S:=(T-TP)/H;
S2:=S*S; S3:=S2*S; S4:=S3*S;
C[3]:=(S2-S)/2;
C[4]:=-S3/6+S2/2-S/3;
C[5]:=S4/24-S3/4+11*S2/24-S/4;
"FOR" I:=1 "STEP" 1 "UNTIL" N "DO"
X[I]:=Y[1,I]-S*Y[2,I]+C[3]*Y[3,I]+C[4]*Y[4,I]+C[5]*Y[5,I];
OUTPUT(61,"("3ZD.2D2B,D.D"-D2B,2(+D.8D"-D2B),2(4ZD),3ZD.2D,
/")",TP,ERR[3],X[1],X[2],NFE,NJE,CLOCK-TIME);
"IF" TP<TEND "THEN"
"BEGIN" POINT:=POINT+1; TP:=PRINT[POINT];
"IF" TP<=T "THEN" "GOTO" NEXT
"END"
"END" CONTROL;
N:=3; NJE:=NFE:=0; T:=0; TEND:=400; EPS:="-5; HMAX:=400;
Y[1]:=Y[2]:=Y[3]:=0; SW[1]:=SW[2]:=SW[3]:=1;
PRINT[1]:=.1; PRINT[2]:=1; PRINT[3]:=10; PRINT[4]:=100;
PRINT[5]:=400; POINT:= 0;
LIPEST(L,Y,"-5,T,F,N);
H2:=(EPS*320)**(1/5)/(4*L);
OUTPUT(61,"(""("EPS=")",D.2D"-D,/,"("INTERVAL OF INTEGRATION=(")",
3ZD,"(",")",3ZD,"(")")",/,"("MAXIMALLY ALLOWED STEPSIZE=")",
D.2D"-D,//")",EPS,T,TEND,HMAX);
OUTPUT(61,"(""("LIPSCHCONST=")",BD.3D"+D,/,"("STARTING STEPSIZE")"
"("=")",BD.2D"+D,/,"("FUNCTIONAL EVAL=")",4ZD,//")",L,H2,NFE);
TIME:=CLOCK;
OUTPUT(61,"(""(" X ERROR Y[1] Y[2]")"
"(" NFE NJE TIME")",/")");
IMPEX(N,T,TEND,Y,F,AVAILABLE,H2,HMAX,"FALSE",EPS,SW,UPDATE,FAIL,
CONTROL);
OUTPUT(61,"("/"("NO OF FUNCTIONAL EVALUATIONS= ")",3ZD,/,
"("NO OF JACOBEAN EVALUATIONS= ")",3ZD,/")",NFE,NJE)
"END"
THIS PROGRAM DELIVERS:
EPS=1.00"-5
INTERVAL OF INTEGRATION=( 0, 400)
MAXIMALLY ALLOWED STEPSIZE=4.00" 2
LIPSCHCONST= 6.003"+1
STARTING STEPSIZE= 1.32"-3
FUNCTIONAL EVAL= 7
X ERROR Y[1] Y[2] NFE NJE TIME
0.00 0.0" 0 +0.00000000" 0 +0.00000000" 0 7 0 0.01
0.10 6.3"-7 +1.49614151"-6 +1.74013792"-4 46 4 0.72
1.00 1.5"-6 +1.91041887"-4 +2.08361269"-3 85 8 1.48
10.00 8.7"-7 +1.30147663"-2 +2.34487800"-2 119 9 1.99
100.00 1.3"-5 +3.06302487"-1 +3.27552180"-1 225 13 3.47
400.00 1.4"-5 +2.22406546" 1 +2.71090507" 1 556 30 7.51
NO OF FUNCTIONAL EVALUATIONS= 556
NO OF JACOBEAN EVALUATIONS= 30
SOURCE TEXT(S):
"CODE" 33135;