NUMAL Section 5.2.1.1.1.2.D
BEGIN SECTION : 5.2.1.1.1.2.D (August, 1974)
AUTHOR: K.DEKKER.
INSTITUTE: MATHEMATICAL CENTRE.
RECEIVED: 1973/07/16.
BRIEF DESCRIPTION:
LINIGER2 SOLVES INITIAL VALUE PROBLEMS , GIVEN AS AN AUTONOMOUS
SYSTEM OF FIRST ORDER DIFFERENTIAL EQUATIONS , BY MEANS OF AN
EXPONENTIALLY FITTED ONESTEP METHOD.
NO 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,
EXPONENTIAL FITTING,
IMPLICIT ONESTEP METHODS.
CALLING SEQUENCE :
THE HEADING OF THE PROCEDURE LINIGER2 READS:
"PROCEDURE" LINIGER2(X,XE,M,Y,SIGMA1,SIGMA2,F,EVALUATE,J,
JACOBIAN,K,ITMAX,STEP,AETA,RETA,OUTPUT);
"INTEGER" M,K,ITMAX;
"REAL" X,XE,SIGMA1,SIGMA2,STEP,AETA,RETA;
"ARRAY" Y,J;
"BOOLEAN" "PROCEDURE" EVALUATE;
"REAL" "PROCEDURE" F;
"PROCEDURE" JACOBIAN,OUTPUT;
"CODE" 33131;
THE MEANING OF THE FORMAL PARAMETERS IS:
X: <VARIABLE>;
THE INDEPENDENT VARIABLE X;
CAN BE USED IN F, 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>;
"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;
SIGMA1: <ARITHMETIC EXPRESSION>;
THE MODULUS OF THE POINT AT WHICH EXPONENTIAL FITTING IS
DESIRED; THIS POINT MAY BE COMPLEX OR REAL AND NEGATIVE;
SIGMA2: <ARITHMETIC EXPRESSION>;
SIGMA2 MAY DEFINE THREE DIFFERENT TYPES OF EXPONENTIAL
FITTING: FITTING IN TWO COMPLEX CONJUGATED POINTS , FITTING
IN TWO REAL NEGATIVE POINTS , OR FITTING IN ONE POINT
COMBINED WITH THIRD ORDER ACCURACY;
IF THIRD ORDER ACCURACY IS DESIRED , SIGMA2 SHOULD HAVE THE
VALUE 0;
IF FITTING IN A SECOND NEGATIVE POINT IS DESIRED , SIGMA2
SHOULD HAVE THE VALUE OF THE MODULUS OF THIS POINT;
IF FITTING IN TWO COMPLEX CONJUGATED POINTS IS DESIRED ,
THEN SIGMA2 SHOULD BE MINUS THE VALUE OF THE ARGUMENT OF
THE POINT IN THE SECOND QUADRANT (THUS A VALUE BETWEEN - PI
AND - PI/2);
F: <PROCEDURE IDENTIFIER>;
THE HEADING OF THIS PROCEDURE READS:
"REAL" "PROCEDURE" F(I); "INTEGER" I;
THIS PROCEDURE SHOULD DELIVER THE RIGHT HAND SIDE OF THE
I-TH DIFFERENTIAL EQUATION AS F;
EVALUATE: <PROCEDURE IDENTIFIER>;
THE HEADING OF THIS PROCEDURE READS:
"BOOLEAN" "PROCEDURE" EVALUATE(ITNUM); "INTEGER" ITNUM;
EVALUATE SHOULD HAVE THE VALUE "TRUE", IF IT IS DESIRED
THAT THE JACOBIAN OF THE SYSTEM IS UPDATED IN THE ITNUM-TH
ITERATION STEP OF THE NEWTON PROCESS;
J: <ARRAY IDENTIFIER>;
"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 HAS TO BE ASSIGNED TO THE
THE ARRAY J , OR AN APPROXIMATION OF THE JACOBIAN , IF ONLY
SECOND ORDER ACCURACY IS REQUIRED;
K: <VARIABLE>;
COUNTS THE NUMBER OF INTEGRATION STEPS TAKEN;
FOR EXAMPLE, CAN BE USED IN EVALUATE;
ITMAX: <ARITHMETIC EXPRESSION>;
AN UPPERBOUND FOR THE NUMBER OF ITERATIONS IN NEWTON'S
PROCESS, USED TO SOLVE THE IMPLICIT EQUATIONS;
STEP: <ARITHMETIC EXPRESSION>;
THE LENGTH OF THE INTEGRATION STEP, TO BE PRESCRIBED BY THE
THE USER;
AETA: <ARITHMETIC EXPRESSION>;
REQUIRED ABSOLUTE PRECISION IN THE NEWTON PROCESS, USED TO
SOLVE THE IMPLICIT EQUATIONS;
RETA: <ARITHMETIC EXPRESSION>;
REQUIRED RELATIVE PRECISION IN THE NEWTON PROCESS, USED TO
SOLVE THE IMPLICIT EQUATIONS;
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 THE PARAMETERS , FOR
EXAMPLE X, K, Y.
DATA AND RESULTS: SEE EXAMPLE OF USE, AND REF[3].
PROCEDURES USED:
VECVEC= CP34010,
MATVEC= CP34011,
MATMAT= CP34013,
DEC = CP34300,
SOL = CP34051.
REQUIRED CENTRAL MEMORY:
EXECUTION FIELD LENGTH: CIRCA 20 + M * (4+M).
RUNNING TIME: DEPENDS STRONGLY ON THE DIFFERENTIAL EQUATION TO SOLVE.
LANGUAGE: ALGOL 60.
METHOD AND PERFORMANCE:
LINIGER2: INTEGRATES THE SYSTEM OF DIFFERENTIAL EQUATIONS FROM X0
UNTIL XE, BY MEANS OF A SECOND ORDER FORMULA (IF SIGMA2=0
AND EVALUATE="TRUE" EVEN THIRD ORDER).
SEE ALSO REF[1] AND REF[2].
REFERENCES:
[1]. W.LINIGER AND R.A.WILLOUGHBY.
EFFICIENT INTEGRATION METHODS FOR STIFF SYSTEMS OF ORDINARY
DIFFERENTIAL EQUATIONS.
SIAM J. NUM. ANAL. 7 (1970) PAGE 47.
[2]. T.J.DEKKER, P.W.HEMKER AND P.W.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.
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 2487 AND Y[2] = .433 710 353 5768.
THE FOLLOWING PROGRAM SHOWS SOME DIFFERENT CALLS OF THE PROCEDURE
LINIGER2, AND ILLUSTRATES THE ACCURACY WHICH MAY BE OBTAINED BY IT:
"BEGIN"
"INTEGER" K,ITMAX,PASSES,PASJAC;
"REAL" X,SIGMA,STEP,TIME;
"REAL" "ARRAY" Y[1:2],J[1:2,1:2];
"REAL" "PROCEDURE" F(I); "INTEGER" I;
"IF" I=1 "THEN" F:=(Y[1]+.99)*(Y[2]-1)+.99 "ELSE"
"BEGIN" PASSES:=PASSES+1; F:=1000*((1+Y[1])*(1-Y[2])-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;
"BOOLEAN" "PROCEDURE" EVALUATE1(I); "INTEGER" I;
EVALUATE1:= I=1;
"BOOLEAN" "PROCEDURE" EVALUATE2(I); "INTEGER" I;
EVALUATE2:= "TRUE";
"PROCEDURE" OUT;
"IF" X=50 "THEN"
OUTPUT(61,"("3(-4ZDB),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 LINIGER2 -SECOND ORDER- ARE:")",/,
"(" K DER.EV. JAC.EV. Y[1] Y[2]")"
"(" TIME")",/")");
"FOR" STEP:=10,1 "DO"
"FOR" ITMAX:=1,3 "DO"
"BEGIN" PASSES:=PASJAC:=0; X:=Y[2]:=0; Y[1]:=1; TIME:=CLOCK;
LINIGER2(X,50.0,2,Y,SIGMA,0.0,F,EVALUATE1,J,JACOBIAN,K,ITMAX,
STEP,"-4,"-4,OUT);
"END";
OUTPUT(61,"("//,
"(" THE RESULTS WITH LINIGER2 -THIRD ORDER- ARE:")",/,
"(" K DER.EV. JAC.EV. Y[1] Y[2]")"
"(" TIME")",/")");
"FOR" STEP:=10,1 "DO"
"FOR" ITMAX:=1,3 "DO"
"BEGIN" PASSES:=PASJAC:=0; X:=Y[2]:=0; Y[1]:=1; TIME:=CLOCK;
LINIGER2(X,50.0,2,Y,SIGMA,0.0,F,EVALUATE2,J,JACOBIAN,K,ITMAX,
STEP,"-4,"-4,OUT);
"END";
"END"
THIS LINE AND THE FOLLOWING TEXT IS PRINTED BY THIS PROGRAM:
THE RESULTS WITH LINIGER2 -SECOND ORDER- ARE:
K DER.EV. JAC.EV. Y[1] Y[2] TIME
5 5 5 +.766 392 210 +.434 218 863 0.092
5 15 5 +.765 755 853 +.433 671 223 0.175
50 50 50 +.765 884 310 +.433 715 687 0.949
50 101 50 +.765 877 388 +.433 710 059 1.494
THE RESULTS WITH LINIGER2 -THIRD ORDER- ARE:
K DER.EV. JAC.EV. Y[1] Y[2] TIME
5 5 5 +.766 392 210 +.434 218 863 0.092
5 15 15 +.765 882 250 +.433 711 614 0.300
50 50 50 +.765 884 310 +.433 715 687 0.949
50 101 101 +.765 878 873 +.433 710 531 2.080
SOURCE TEXT(S):
"CODE" 33131;