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;