NUMAL Section 5.2.1.1.1.1.F

BEGIN SECTION : 5.2.1.1.1.1.F (February, 1979)

PROCEDURE : MULTISTEP.

AUTHOR: P.W.HEMKER.

INSTITUTE: MATHEMATICAL CENTRE.

RECEIVED: 730515.

BRIEF DESCRIPTION:

    MULTISTEP  SOLVES  AN INITIAL VALUE  PROBLEM, FOR A SYSTEM OF
    FIRST ORDER ORDINARY DIFFERENTIAL EQUATIONS  DY = F(Y).
    IN PARTICULAR THIS PROCEDURE IS SUITABLE FOR THE INTEGRATION OF
    STIFF DIFFERENTIAL EQUATIONS. IT CAN ALSO BE USED FOR
    NON-STIFF PROBLEMS.

KEYWORDS:

    INITIAL VALUE PROBLEM,
    SYSTEM OF FIRST ORDER ORDINARY DIFFERENTIAL EQUATIONS.
    STIFF EQUATIONS,

CALLING SEQUENCE:

    THE HEADING OF THE PROCEDURE READS:
    "BOOLEAN" "PROCEDURE" MULTISTEP(X,XEND,Y,HMIN,HMAX,YMAX,EPS,
                      FIRST,SAVE,DERIV,AVAILABLE,JACOBIAN,STIFF,N,OUT);
    "VALUE" HMIN,HMAX,EPS,XEND,N,STIFF;
    "BOOLEAN" FIRST,AVAILABLE,STIFF;
    "INTEGER" N;
    "REAL" X,XEND,HMIN,HMAX,EPS;
    "ARRAY" Y,YMAX,SAVE,JACOBIAN;
    "PROCEDURE" DERIV,OUT;
    "CODE" 33080;

    MULTISTEP  DELIVERS THE FOLLOWING BOOLEAN VALUE:
            IF DIFFICULTIES  ARE  ENCOUNTERED  DURING  THE  INTEGRATION
            (I.E. SAVE[-1] ^= 0 "OR" SAVE[-2] ^= 0 )  MULTISTEP IS  SET
            TO "FALSE", OTHERWISE MULTISTEP IS SET "TRUE".

    THE MEANING OF THE FORMAL PARAMETERS IS:
    X:      <VARIABLE>;
            THE INDEPENDENT VARIABLE X.
            CAN BE USED IN DERIV, AVAILABLE ETC.;
            ENTRY: THE INITIAL VALUE X0;
            EXIT : THE FINAL VALUE 'XEND';
    XEND:   <ARITHMETIC EXPRESSION>;
            THE FINAL VALUE OF X ( XEND >= X );
    Y:      <ARRAY IDENTIFIER>;
            "ARRAY" Y[1:6*N];
            THE DEPENDENT VARIABLE;
            ENTRY  Y[1:N] : THE INITIAL VALUES OF THE SOLUTION OF THE
                   SYSTEM OF DIFFERENTIAL EQUATIONS AT X = X0;
            EXIT Y[1:N] : THE FINAL VALUES OF THE SOLUTION AT X = XEND;
    HMIN,HMAX: <ARITHMETIC EXPRESSION>;
            ENTRY : MINIMAL  RESP. MAXIMAL STEPLENGTH  ALLOWED;
    YMAX:   <ARRAY IDENTIFIER>;
            "ARRAY" YMAX[1:N];
            ENTRY: THE  ABSOLUTE  LOCAL  ERROR  BOUND  DIVIDED  BY  EPS
            EXIT : YMAX[1]  GIVES THE  MAXIMAL VALUE OF THE ENTRY VALUE
                   OF  YMAX[I]  AND  THE  VALUES  OF  ABS(Y[I])  DURING
                   INTEGRATION;
    EPS:    <ARITHMETIC EXPRESSION>;
            THE RELATIVE LOCAL ERROR BOUND;
    FIRST:  <IDENTIFIER>;
            IF FIRST = "TRUE" THEN THE PROCEDURE STARTS THE INTEGRATION
            WITH  A  FIRST  ORDER  ADAMS METHOD AND A STEPLENGTH  EQUAL
            TO HMIN. UPON COMPLETION OF A CALL FIRST:= "FALSE";
            IF   FIRST  =  "FALSE"   THEN   THE   PROCEDURE   CONTINUES
            INTEGRATION;

    SAVE:   <ARRAY IDENTIFIER>;
            "ARRAY" SAVE[-38:6*N];
            IN THIS ARRAY THE PROCEDURE STORES INFORMATION WHICH CAN BE
            USED IN A CONTINUING CALL WITH FIRST="FALSE";
            BESIDES THE FOLLOWING MESSAGES ARE DELIVERED:
            SAVE[ 0]=0 : AN ADAMS METHOD HAS BEEN USED;
                     1 : THE PROCEDURE SWITCHED TO GEARS METHOD;
            SAVE[-1]=0 : NO ERROR MESSAGE;
                     1 : WITH  THE HMIN SPECIFIED  THE PROCEDURE CANNOT
                         HANDLE  THE  NONLINEARITY  (DECREASE HMIN!  );
            SAVE[-2]  NUMBER OF TIMES  THAT THE  REQUESTED  LOCAL ERROR
                      BOUND WAS EXCEEDED;
            SAVE[-3]  IF SAVE[-2] IS NONZERO THEN SAVE[-3] GIVES AN
                      ESTIMATE OF THE MAXIMAL LOCAL ERROR BOUND,
                      OTHERWISE SAVE[-3]=0;
    DERIV:  <PROCEDURE IDENTIFIER>;
            THE HEADING OF THIS PROCEDURE READS:
            "PROCEDURE" DERIV(DF); "ARRAY" DF;
            THIS PROCEDURE SHOULD DELIVER DY[I]/DX IN DF[I];
    AVAILABLE: <BOOLEAN EXPRESSION>;
            IF AN  ANALYTICAL EXPRESSION  OF THE JACOBIAN MATRIX IS NOT
            AVAILABLE THIS EXPRESSION IS SET TO "FALSE";
            OTHERWISE THIS EXPRESSION IS SET TO "TRUE" AND THE
            EVALUATION  OF  THIS  BOOLEAN  EXPRESSION  MUST  EFFECT THE
            FOLLOWING SIDE-EFFECT:
            THE ENTRIES OF THE JACOBIAN MATRIX   D(DY[I]/DX)/DY[J]  ARE
            DELIVERED IN THE ARRAY ELEMENTS  JACOBIAN[I,J];
    JACOBIAN: <ARRAY IDENTIFIER>;
            "ARRAY" JACOBIAN[1:N,1:N];
            AT EACH EVALUATION OF THE BOOLEAN EXPRESSION AVAILABLE WITH
            THE RESULT AVAILABLE:="TRUE", THE JACOBIAN MATRIX HAS TO BE
            ASSIGNED TO THIS ARRAY (SEE THE  EXAMPLE OF USE);
    STIFF:  <BOOLEAN EXPRESSION>;
            IF STIFF = "TRUE"  THE PROCEDURE  SKIPS AN ATTEMPT TO SOLVE
            THE PROBLEM WITH ADAMS-BASHFORTH- OR ADAMS-MOULTON
            METHODS, DIRECTLY USING GEARS METHOD;
    N:      <ARITHMETIC EXPRESSION>;
            THE NUMBER OF EQUATIONS;
    OUT:    <PROCEDURE IDENTIFIER>;
            THE HEADING OF THIS PROCEDURE READS:
            "PROCEDURE" OUT(H,K);  "VALUE" H,K;  "REAL" H; "INTEGER" K;
            AT THE END OF EACH ACCEPTED STEP OF THE INTEGRATION PROCESS
            THIS PROCEDURE IS CALLED. THE LAST STEPLENGTH USED (H)  AND
            THE ORDER OF THE METHOD (K) ARE DELIVERED.
            AT EACH CALL OF THE  PROCEDURE OUT, THE  CURRENT VALUES  OF
            THE INDEPENDENT VARIABLE (X) AND OF THE SOLUTION (Y[I](X) )
            ARE AVAILABLE  FOR USE. MOREOVER, IN THE  NEIGHBOURHOOD  OF
            THE CURRENT VALUE OF X, ANY VALUE OF  Y[I](X SPECIFIED) CAN
            BE COMPUTED BY MEANS OF THE FOLLOWING INTERPOLATION FORMULA
            Y[I](X SPECIFIED) =
                 SUM(J,0,K, Y[I+J*N] * ((X SPECIFIED - X)/H) ** J ).

PROCEDURES USED:

    MATVEC = CP34011 ,
    DEC    = CP34300 ,
    SOL    = CP34051 .

REQUIRED CENTRAL MEMORY: CIRCA  N * ( 2 * N + 5 )  MEMORY PLACES.

METHOD AND PERFORMANCE :
    MULTISTEP IS BASED ON TWO LINEAR MULTISTEP METHODS. FOR STIFF
    PROBLEMS IT USES THE BACKWARD DIFFERENTIATION METHODS, FOR
    FOR NON-STIFF PROBLEMS THE ADAMS-BASHFORTH-MOULTON METHODS.
    MULTISTEP IS PROVIDED WITH ORDER, STEPSIZE AND ERROR CONTROL.

REFERENCES:

    [1].P.W.HEMKER.
        AN ALGOL 60 PROCEDURE  FOR THE SOLUTION  OF STIFF  DIFFERENTIAL
        EQUATIONS.
        MATH. CENTRE, AMSTERDAM. REPORT MR 128/71;

EXAMPLE OF USE:

    THE SOLUTION AT X=1 AND AT X=10 OF THE DIFFERENTIAL EQUATIONS:
    DY[1]/DX = 0.04 * (1-Y[1]-Y[2]) - Y[1] * ("4*Y[2] + 3"7*Y[1])
    DY[2]/DX = 3"7 Y[1]**2
    WITH THE INITIAL CONDITIONS AT  X = 0 :
    Y[1] = 0  AND Y[2] = 0
    MAY BE OBTAINED BY THE FOLLOWING PROGRAM:

 "BEGIN"

        "BOOLEAN" FIRST;
        "INTEGER" I,J,CF,CJ,CA;
        "REAL" X,XEND,HMIN,EPS,R;
        "ARRAY" Y[1:12],YMAX[1:2],D[-40:12],JAC[1:2,1:2];

        "PROCEDURE" DER (F);"ARRAY" F;
        "BEGIN" "REAL" R; CF:=CF+1;
             F[2]:= R:= 3"7*Y[1]*Y[1];
             F[1]:= 0.04*(1-Y[1]-Y[2]) - "4*Y[1]*Y[2] - R;
        "END" F;

        "BOOLEAN" "PROCEDURE" AVAIL;
        "BEGIN" "REAL" R; CJ:= CJ+1;
            AVAIL:= "TRUE";
            JAC[2,1]:= R:= 6"7*Y[1];
            JAC[1,1]:= -0.04 - "4*Y[2] - R;
            JAC[1,2]:= -0.04 - "4*Y[1];
            JAC[2,2]:= 0
        "END" JAC AVAIL;

        "PROCEDURE" OUT(H,K);
        "VALUE" H, K; "REAL" H; "INTEGER" K; CA:= CA+1;

        LABEL:
        OUTPUT(61,"("/,"("HMIN,EPS?")",/")");
        INREAL(70,HMIN); INREAL(70,EPS);
        "IF" HMIN<0 "THEN" "GOTO" ESCAPE;

        FIRST:= "TRUE"; CA:=CF:=CJ:=0;
        X:=0; Y[1]:= Y[2]:= 0;
        YMAX[1]:= 0.0001; YMAX[2]:= 1;

        "FOR" XEND:= 1, 10 "DO"
        "BEGIN"
            MULTISTEP(X,XEND,Y,HMIN,5,YMAX,EPS,FIRST,D,DER,AVAIL,
                                               JAC,"TRUE",2,OUT);
            OUTPUT(61,"("3(5ZD,2B),2(+.13D"+2D,2B),/")",
                                             CA,CF,CJ,Y[1],Y[2]);
        "END";

        "GOTO" LABEL;
        ESCAPE:
    "END"

    IT DELIVERS WITH  HMIN = "-10  AND  EPS = "-9:
    240   648   2   +.3074626[602000]"-04   +.3350951[493111]"-01
    315   902   3   +.16233909[62091]"-04   +.15861383[92015]"+00
    (NON-SIGNIFICANT DIGITS ARE PLACED BETWEEN [ ] ).

SOURCE TEXT(S):
"CODE" 33080;