NUMAL Section 5.2.1.1.1.3.B
BEGIN SECTION : 5.2.1.1.1.3.B (August, 1974)
AUTHORS: P.J. VAN DER HOUWEN AND K.DEKKER.
INSTITUTE: MATHEMATICAL CENTRE.
RECEIVED: 740416.
BRIEF DESCRIPTION:
EXPONENTIALLY FITTED TAYLOR SOLVES AN INITIAL VALUE PROBLEM, GIVEN
AS A SYSTEM OF FIRST ORDER DIFFERENTIAL EQUATIONS , BY MEANS OF A
ONE-STEP TAYLOR-METHOD . AUTOMATIC STEPSIZE CONTROL IS PROVIDED.
IN PARTICULAR THIS METHOD IS SUITABLE FOR THE INTEGRATION OF STIFF
DIFFERENTIAL EQUATIONS , PROVIDED THAT HIGHER ORDER DERIVATIVES CAN
BE EASILY OBTAINED.
KEYWORDS:
DIFFERENTIAL EQUATIONS,
INITIAL VALUE PROBLEMS,
EXPONENTIAL FITTING,
STIFF EQUATIONS,
THREE-CLUSTER METHOD,
ONE-STEP TAYLOR-METHOD.
CALLING SEQUENCE:
THE HEADING OF THE PROCEDURE READS:
"PROCEDURE" EXPONENTIALLY FITTED TAYLOR (T, TE, M0, M, U, SIGMA,
PHI, DIAMETER, DERIVATIVE, I, K, ALFA, NORM,
AETA, RETA, ETA, RHO, HMIN, HSTART, OUTPUT);
"INTEGER" M0,M,I,K,NORM;
"REAL" T,TE,SIGMA,PHI,DIAMETER,ALFA,AETA,RETA,ETA,RHO,HMIN,HSTART;
"ARRAY" U;
"PROCEDURE" DERIVATIVE,OUTPUT;
"CODE" 33050;
THE MEANING OF THE FORMAL PARAMETERS IS:
T: <VARIABLE>;
THE INDEPENDENT VARIABLE T;
MAY BE USED IN DERIVATIVE, SIGMA ETC.;
ENTRY: THE INITIAL VALUE T0;
EXIT : THE FINAL VALUE TE;
TE: <ARITHMETIC EXPRESSION>;
THE FINAL VALUE OF T (TE >= T);
M0: <ARITHMETIC EXPRESSION>;
INDEX OF THE FIRST EQUATION OF THE SYSTEM TO BE SOLVED;
M: <ARITHMETIC EXPRESSION>;
INDEX OF THE LAST EQUATION OF THE SYSTEM TO BE SOLVED;
U: <ARRAY IDENTIFIER>;
"ARRAY" U[M0:M];
THE DEPENDENT VARIABLE;
ENTRY: THE INITIAL VALUES OF THE SOLUTION OF THE SYSTEM OF
DIFFERENTIAL EQUATIONS AT T = T0;
EXIT : THE VALUES OF THE SOLUTION AT T = TE;
SIGMA: <ARITHMETIC EXPRESSION>;
THE MODULUS OF THE (COMPLEX) POINT AT WHICH EXPONENTIAL
FITTING IS DESIRED , FOR EXAMPLE AN APPROXIMATION OF THE
MODULUS OF THE CENTRE OF THE LEFT HAND CLUSTER;
PHI: <ARITHMETIC EXPRESSION>;
THE ARGUMENT OF THE (COMPLEX) POINT AT WHICH EXPONENTIAL
FITTING IS DESIRED;
PHI SHOULD HAVE A VALUE FROM THE RANGE [PI/2,PI];
DIAMETER: <ARITHMETIC EXPRESSION>;
THE DIAMETER OF THE LEFT HAND CLUSTER;
DERIVATIVE: <PROCEDURE IDENTIFIER>;
THE HEADING OF THIS PROCEDURE READS:
"PROCEDURE" DERIVATIVE(I,A); "INTEGER" I; "ARRAY" A;
I ASSUMES THE VALUES 1,2,3 AND A IS A ONE-DIMENSIONAL ARRAY
A[M0:M];
WHEN THIS PROCEDURE IS CALLED , ARRAY A CONTAINS THE
COMPONENTS OF THE (I-1)-ST DERIVATIVE OF U AT THE POINT T;
UPON COMPLETION OF DERIVATIVE, ARRAY A SHOULD CONTAIN THE
COMPONENTS OF THE I-TH DERIVATIVE OF U AT THE POINT T;
I: <VARIABLE>;
A JENSEN PARAMETER FOR PROCEDURE DERIVATIVE;
MAY BE USED IN M0 AND M;
K: <VARIABLE>;
INDICATES THE NUMBER OF INTEGRATION STEPS PERFORMED;
ENTRY: K = 0;
EXIT : THE NUMBER OF INTEGRATION STEPS PERFORMED;
ALFA: <ARITHMETIC EXPRESSION>;
MAXIMAL GROWTH FACTOR FOR THE INTEGRATION STEP LENGTH;
NORM: <ARITHMETIC EXPRESSION>;
IF NORM = 1 DISCREPANCY AND TOLERANCE ARE ESTIMATED IN THE
MAXIMUM NORM, OTHERWISE IN THE EUCLIDIAN NORM;
AETA: <ARITHMETIC EXPRESSION>;
DESIRED ABSOLUTE LOCAL ACCURACY ; AETA SHOULD BE POSITIVE;
RETA: <ARITHMETIC EXPRESSION>;
DESIRED RELATIVE LOCAL ACCURACY ; RETA SHOULD BE POSITIVE;
ETA: <VARIA2LE>;
COMPUTED TOLERANCE;
RHO: <VARIABLE>;
COMPUTED DISCREPANCY;
HMIN: <ARITHMETIC EXPRESSION>;
MINIMAL STEPSIZE BY WHICH THE INTEGRATION IS PERFORMED;
HOWEVER, A SMALLER STEP WILL BE TAKEN IF HMIN EXCEEDS THE
STEPSIZE HSTAB , PRESCRIBED BY THE STABILITY CONDITIONS
(SEE REF[2], FORMULA 6.12);
IF HSTAB TENDS TO ZERO, THE PROCEDURE TERMINATES;
HSTART: <VARIABLE>;
ENTRY: THE INTITIAL STEPSIZE ; HOWEVER, IF K = 0 ON ENTRY,
THE VALUE OF HSTART IS NOT TAKEN INTO CONSIDERATION;
EXIT: A SUGGESTION FOR THE STEPSIZE , IF THE INTEGRATION
SHOULD BE CONTINUED FOR T>TE;
HSTART MAY BE USED IN SUCCESSIVE CALLS OF THE PROCEDURE, IN
ORDER TO OBTAIN THE SOLUTION IN SEVERAL POINTS TE1,TE2,ETC;
OUTPUT: <PROCEDURE IDENTIFIER>;
THE HEADING OF THIS PROCEDURE READS:
"PROCEDURE" OUTPUT;
THROUGH THIS PROCEDURE THE VALUES AFTER EACH INTEGRATION
STEP OF FOR INSTANCE T, U, ETA AND RHO ARE ACCESSIBLE;
DATA AND RESULTS:
FOR FURTHER EXPLANATION OF THE PARAMETERS SIGMA,PHI,DIAMETER,AETA,
RETA,ETA,RHO,M0,M SEE REF[2];
FOR RESULTS: SEE EXAMPLE OF USE AND REF[2];
PROCEDURES USED:
INIVEC = CP 31010;
DUPVEC = CP 31030;
VECVEC = CP 34010;
ELMVEC = CP 34020;
ZEROIN = CP 34150.
REQUIRED CENTRAL MEMORY:
EXECUTION FIELD LENGTH: CIRCA 40 + 2 * (M - M0).
RUNNING TIME:
DEPENDS STRONGLY ON THE DIFFERENTIAL EQUATION TO BE SOLVED.
LANGUAGE: ALGOL 60.
METHOD AND PERFORMANCE: SEE REFERENCES.
REFERENCES:
[1]. P.J. VAN DER HOUWEN.
ONE-STEP METHODS FOR LINEAR INITIAL VALUE PROBLEMS II,
POLYNOMIAL METHODS.
TW REPORT 122, (1970) MATHEMATICAL CENTRE.
[2]. P.J. VAN DER HOUWEN, P.BEENTJES, K.DEKKER AND E.SLAGT.
ONE-STEP METHODS FOR LINEAR INITIAL VALUE PROBLEMS III,
NUMERICAL EXAMPLES.
TW REPORT 130, (1971) MATHEMATICAL CENTRE.
EXAMPLE OF USE:
THE SOLUTION AT T=EXP(1) AND T=EXP(2) OF THE DIFFERENTIAL EQUATION
DU/DT=-EXP(T)*(U-LN(T)) + 1/T WITH INITIAL CONDITION U(.01)=LN(.01)
AND ANALYTICAL SOLUTION U(T) = LN(T), MAY BE OBTAINED AS FOLLOWS:
"BEGIN" "INTEGER" I,K;
"REAL" T,TE,TE1,TE2,RETA,ETA,RHO,PI,HS,EXPT,LNT,TIME,U0,U1,U2;
"REAL" "ARRAY" U[0:0];
"PROCEDURE" DERIVATIVE(I,U); "INTEGER" I; "ARRAY" U;
"IF" I=1 "THEN" "BEGIN" EXPT:=EXP(T); LNT:=LN(T); U0:=U[0];
U1:=U[0]:=EXPT*(LNT-U0)+1/T
"END" "ELSE"
"IF" I=2 "THEN" U2:=U[0]:=EXPT*(LNT-U0-U1+1/T)-1/T/T
"ELSE" U[0]:=EXPT*(LNT-U0-2*U1-U2+2/T-1/T/T)+2/T/T/T;
"PROCEDURE" OUT;
"IF" T=TE "THEN" OUTPUT(61,"("6ZD,+3ZD.3DB3DB3D")",K,U[0]);
"PROCEDURE" OUT1;
OUTPUT(61,"("4BD"-D,3Z.3D,/")",RETA,CLOCK-TIME);
OUTPUT(61,"(""(" THIS LINE AND THE FOLLOWING TEXT IS ")"
"("PRINTED BY THIS PROGRAM")",//,
"(" THE RESULTS WITH EFT ARE -CONFER REF[2]- :")",/,
"(" K U(TE1) K U(TE2)")"
"(" RETA TIME")",/")");
PI:=4*ARCTAN(1); TE1:=EXP(1); TE2:=EXP(2);
"FOR" RETA:="-1,"-2,"-3,"-4 "DO"
"BEGIN" T:=.01; U[0]:=LN(T); K:=0; HS:=0; TIME:=CLOCK;
"FOR" TE:=TE1,TE2 "DO"
EXPONENTIALLY FITTED TAYLOR
(T,TE,0,0,U,EXP(T),PI,2*EXP(2*T/3),DERIVATIVE,I,K,1.5,2,
RETA/10,RETA,ETA,RHO,"-4,HS,OUT); OUT1
"END";
OUTPUT(61,"("//,"(" WITH RELAXED ACCURACY CONDITIONS FOR ")"
"("T>3:")",/,"(" K U(TE1) K U(TE2)")"
"(" RETA TIME")",/")");
"FOR" RETA:="-1,"-2,"-3,"-4 "DO"
"BEGIN" T:=.01; U[0]:=LN(T); K:=0; HS:=0; TIME:=CLOCK;
"FOR" TE:=TE1,TE2 "DO"
EXPONENTIALLY FITTED TAYLOR
(T,TE,0,0,U,EXP(T),PI,2*EXP(2*T/3),DERIVATIVE,I,K,1.5,2,
RETA/10*("IF" T<3 "THEN" 1 "ELSE" EXP(2*(T-3))),
RETA*("IF" T<3 "THEN" 1 "ELSE" EXP(2*(T-3))),
ETA,RHO,"-4,HS,OUT); OUT1
"END"
"END"
THIS LINE AND THE FOLLOWING TEXT IS PRINTED BY THIS PROGRAM
THE RESULTS WITH EFT ARE -CONFER REF[2]- :
K U(TE1) K U(TE2) RETA TIME
15 +1.003 845 001 42 +2.000 076 417 1"-1 .938
22 +1.001 211 286 52 +2.000 066 067 1"-2 1.121
36 +1.000 108 738 92 +2.000 020 495 1"-3 1.872
56 +1.000 045 271 171 +2.000 000 925 1"-4 3.493
WITH RELAXED ACCURACY CONDITIONS FOR T>3:
K U(TE1) K U(TE2) RETA TIME
15 +1.003 845 001 42 +2.000 076 417 1"-1 1.037
22 +1.001 211 286 50 +2.000 049 978 1"-2 1.154
36 +1.000 108 738 68 +2.000 023 330 1"-3 1.419
56 +1.000 045 271 98 +2.000 065 056 1"-4 2.008
SOURCE TEXT(S):
"CODE" 33050;