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;