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;