NUMAL Section 5.2.1.1.1.2.C
BEGIN SECTION : 5.2.1.1.1.2.C (October, 1974)
AUTHOR: K.DEKKER.
INSTITUTE: MATHEMATICAL CENTRE.
RECEIVED: 1973/09/01.
BRIEF DESCRIPTION:
LINIGER1VS SOLVES INITIAL VALUE PROBLEMS , GIVEN AS AN AUTONOMOUS
SYSTEM OF FIRST ORDER DIFFERENTIAL EQUATIONS , BY MEANS OF AN
IMPLICIT,FIRST ORDER ACCURATE, EXPONENTIALLY FITTED ONESTEP METHOD.
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 LINIGER1VS READS:
"PROCEDURE" LINIGER1VS (X,XE,M,Y,SIGMA,DERIVATIVE,J,JACOBIAN,
ITMAX,HMIN,HMAX,AETA,RETA,INFO,OUTPUT);
"VALUE" M;
"INTEGER" M,ITMAX;
"REAL" X,XE,SIGMA,HMIN,HMAX,AETA,RETA;
"ARRAY" Y,J,INFO;
"PROCEDURE" DERIVATIVE,JACOBIAN,OUTPUT;
"CODE" 33132;
THE MEANING OF THE FORMAL PARAMETERS IS:
X: <VARIABLE>;
THE INDEPENDENT VARIABLE X;
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;
SIGMA: <ARITHMETIC EXPRESSION>;
THE MODULUS OF THE POINT AT WHICH EXPONENTIAL FITTING IS
DESIRED, FOR EXAMPLE THE LARGEST NEGATIVE EIGENVALUE OF THE
JACOBIAN OF THE SYSTEM OF DIFFERENTIAL EQUATIONS;
DERIVATIVE: <PROCEDURE IDENTIFIER>;
THE HEADING OF THIS PROCEDURE READS:
"PROCEDURE" DERIVATIVE(Y); "ARRAY" Y;
THIS PROCEDURE SHOULD DELIVER THE RIGHT HAND SIDE OF THE
I-TH DIFFERENTIAL EQUATION AT THE POINT (Y) AS Y[I];
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 (AN APPROXIMATION OF) THE JACOBIAN HAS TO
BE ASSIGNED TO THE ARRAY J;
ITMAX: <ARITHMETIC EXPRESSION>;
AN UPPERBOUND FOR THE NUMBER OF ITERATIONS IN NEWTON'S
PROCESS, USED TO SOLVE THE IMPLICIT EQUATIONS;
HMIN: <ARITHMETIC EXPRESSION>;
MINIMAL STEPSIZE BY WHICH THE INTEGRATION IS PERFORMED;
HMAX: <ARITHMETIC EXPRESSION>;
MAXIMAL STEPSIZE BY WHICH THE INTEGRATION IS PERFORMED;
AETA: <ARITHMETIC EXPRESSION>;
REQUIRED ABSOLUTE PRECISION IN THE INTEGRATION PROCESS;
RETA: <ARITHMETIC EXPRESSION>;
REQUIRED RELATIVE PRECISION IN THE INTEGRATION PROCESS;
IF BOTH AETA AND RETA HAVE NEGATIVE VALUES , INTEGRATION
WILL BE PERFORMED WITH A STEPSIZE EQUAL TO HMAX , WHICH MAY
BE VARIATED BY USER ; IN THIS CASE THE ABSOLUTE VALUES OF
AETA AND RETA WILL CONTROL THE NEWTON ITERATION;
INFO: <ARRAY IDENTIFIER>;
"ARRAY" INFO[1:9];
DURING INTEGRATION AND UPON EXIT THIS ARRAY CONTAINS THE
FOLLOWING INFORMATION:
INFO[1]: NUMBER OF INTEGRATION STEPS TAKEN;
INFO[2]: NUMBER OF DERIVATIVE EVALUATIONS USED;
INFO[3]: NUMBER OF JACOBIAN EVALUATIONS USED;
INFO[4]: NUMBER OF INTEGRATION STEPS EQUAL TO HMIN TAKEN ;
INFO[5]: NUMBER OF INTEGRATION STEPS EQUAL TO HMAX TAKEN ;
INFO[6]: MAXIMAL NUMBER OF ITERATIONS TAKEN IN THE NEWTON
PROCESS;
INFO[7]: LOCAL ERROR TOLERANCE;
INFO[8]: ESTIMATED LOCAL ERROR;
INFO[9]: MAXIMUM VALUE OF THE ESTIMATED LOCAL ERROR;
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, Y, INFO.
DATA AND RESULTS: SEE EXAMPLE OF USE, AND REF[2].
PROCEDURES USED:
INIVEC= CP31010,
MULVEC= CP31020,
MULROW= CP31021,
DUPVEC= CP31030,
MATVEC= CP34011,
ELMVEC= CP34020,
VECVEC= CP34010,
DEC = CP34300,
SOL = CP34051.
REQUIRED CENTRAL MEMORY:
EXECUTION FIELD LENGTH : CIRCA 20 + M * (5 + M).
RUNNING TIME: DEPENDS STRONGLY ON THE DIFFERENTIAL EQUATION TO SOLVE.
LANGUAGE: ALGOL 60.
METHOD AND PERFORMANCE:
LINIGER1VS: INTEGRATES THE SYSTEM OF DIFFERENTIAL EQUATIONS FROM X0
UNTIL XE, BY MEANS OF A FIRST ORDER FORMULA.
THE INTEGRATION METHOD IS BASED ON THE F(1) FORMULA DESCRIBED BY
LINIGER AND WILLOUGHBY (SEE REF[1]). ERROR ESTIMATES AND A STEPSIZE
STRATEGY FOR THIS METHOD ARE DESCRIBED IN [2] , AND A VARIABLE STEP
METHOD IS CONSTRUCTED FOR THE CONVENIENCE OF THE USER. HOWEVER, THE
STEPSIZE STRATEGY REQUIRES MANY EXTRA ARRAY OPERATIONS. THE USER
MAY AVOID THIS EXTRA WORK BY GIVING AETA AND RETA A NEGATIVE VALUE
AND PRESCRIBING A STEPSIZE (HMAX) HIMSELF.
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]. K.DEKKER.
ERROR ESTIMATES AND STEPSIZE STRATEGIES FOR THE LINIGER-
WILLOUGHBY FORMULAE.
(TO APPEAR IN 1974).
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.
THE SOLUTION AT X = 50 IS APPROXIMATELY:
Y[1] = .765 878 320 2487 AND Y[2] = .433 710 353 5768.
THE FOLLOWING PROGRAM SHOWS INTEGRATION OF THIS PROBLEM WITH
VARIABLE AND CONSTANT STEPSIZES:
"BEGIN" "COMMENT" TEST LINIGER1VS;
"INTEGER" ITMAX;
"REAL" X,SIGMA,RETA,TIME;
"REAL" "ARRAY" Y[1:2],J[1:2,1:2],INFO[1:9];
"PROCEDURE" F(A); "ARRAY" A;
"BEGIN" "REAL" A1,A2; A1:=A[1]; A2:=A[2];
A[1]:=(A1+.99)*(A2-1)+.99;
A[2]:=1000*((1+A1)*(1-A2)-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;
"END" JACOBIAN;
"PROCEDURE" OUT;
"IF" X=50 "THEN"
OUTPUT(61,"("6(3ZDB),2BD"-ZD,2(2B+.3DB3D),-3ZD.3D,/")",
INFO[1],INFO[2],INFO[3],INFO[4],INFO[5],INFO[6],INFO[9],Y[1],
Y[2],CLOCK-TIME);
"FOR" RETA:="-2,"-4,"-6 "DO"
"BEGIN" X:=Y[2]:=0; Y[1]:=1; TIME:=CLOCK;
LINIGER1VS(X,50.0,2,Y,SIGMA,F,J,JACOBIAN,10,.1,50.0,RETA,
RETA,INFO,OUT);
"END"; OUTPUT(61,"("//")");
"FOR" RETA:=-"-2,-"-4,-"-6 "DO"
"BEGIN" X:=Y[2]:=0; Y[1]:=1; TIME:=CLOCK;
LINIGER1VS(X,50.0,2,Y,SIGMA,F,J,JACOBIAN,10,.1,1.0,RETA,
RETA,INFO,OUT);
"END";
"END"
17 21 8 2 0 2 2" -2 +.772 017 +.435 672 0.525
13 25 23 2 0 2 2" -2 +.767 414 +.434 202 0.717
105 210 105 2 0 2 2" -2 +.766 027 +.433 758 4.741
50 52 1 0 50 2 0" 0 +.766 670 +.433 081 0.549
50 104 3 0 50 3 0" 0 +.766 183 +.433 811 1.158
50 152 12 0 50 4 0" 0 +.766 185 +.433 809 1.653
SOURCE TEXT(S):
"CODE" 33132;