NUMAL Section 5.2.1.1.1.1.H

BEGIN SECTION : 5.2.1.1.1.1.H (February, 1979)

PROCEDURE : ARK.

AUTHOR:  P.A. BEENTJES.

INSTITUTE:  MATHEMATICAL CENTRE.

RECEIVED: 740510.

BRIEF DESCRIPTION:

    ARK SOLVES AN INITIAL VALUE PROBLEM, FOR A SYSTEM  OF FIRST ORDER
    ORDINARY DIFFERENTIAL EQUATIONS . ARK IS RECOMMENDED FOR THE
    INTEGRATION OF SEMI-DISCRETE PARABOLIC AND HYPERBOLIC
    INITIAL-BOUNDARY PROBLEMS.

KEYWORDS:
    INITIAL VALUE PROBLEM,
    SEMI-DISCRETE PARABOLIC AND HYPERBOLIC PROBLEM.

CALLING SEQUENCE:

    THE HEADING OF THE PROCEDURE READS:
    "PROCEDURE" ARK (T, TE, M0, M, U, DERIVATIVE, DATA, OUT);
    "INTEGER" M0, M; "REAL" T, TE; "ARRAY" U, DATA;
    "PROCEDURE" DERIVATIVE, OUT;
    "CODE" 33061;

    ARK :   INTEGRATES THE SYSTEM OF ORDINARY DIFFERENTIAL  EQUATIONS
            DU / DT = H(T, U),  U = U0   AT   T = T0.

    THE MEANING OF THE FORMAL PARAMETERS IS:
    T:      <VARIABLE>;
            THE  INDEPENDENT  VARIABLE  T;
            ENTRY: THE INITIAL VALUE T0;
            EXIT : THE FINAL VALUE TE;
    TE:     <ARITHMETIC EXPRESSION>;
            THE FINAL VALUE OF T (TE >= T);
    M0,M:   <ARITHMETIC EXPRESSION>;
            INDICES  OF  THE  FIRST  AND  LAST EQUATION OF THE SYSTEM;
    U:      <ARRAY IDENTIFIER>;
            "ARRAY" U[M0 : M];
            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;

    DERIVATIVE: <PROCEDURE IDENTIFIER>;
            THE HEADING OF THIS PROCEDURE READS:
            "PROCEDURE" DERIVATIVE(T, V); "REAL" T; "ARRAY" V;
            THIS PROCEDURE  PERFORMS  AN  EVALUATION OF THE RIGHT HAND
            SIDE  OF THE SYSTEM WITH DEPENDENT VARIABLES V[M0 : M] AND
            INDEPENDENT VARIABLE T; UPON COMPLETION OF DERIVATIVE, THE
            RIGHT HAND SIDE   SHOULD  BE   OVERWRITTEN  ON  V[M0 : M];
    DATA:   <ARRAY IDENTIFIER>;
            "ARRAY" DATA[1 : 10 + DATA[1]];
            IN ARRAY DATA ONE SHOULD GIVE:
            DATA[1]: THE   NUMBER   OF   EVALUATIONS  OF   H(T, U) PER
                     INTEGRATION STEP(DATA[1] >= DATA[2]);
            DATA[2]: THE ORDER OF ACCURACY OF THE METHOD (DATA[2]<=3);
            DATA[3]: STABILITY BOUND(SEE REFERENCE [3]);
            DATA[4]: THE SPECTRAL RADIUS  OF THE JACOBIAN MATRIX  WITH
                     RESPECT TO  THOSE EIGENVALUES, WHICH ARE  LOCATED
                     IN THE NON-POSITIVE HALF PLANE;
            DATA[5]: THE MINIMAL STEPSIZE;
            DATA[6]: THE ABSOLUTE TOLERANCE;
            DATA[7]: THE RELATIVE TOLERANCE;
                     IF  BOTH  DATA[6]  AND  DATA[7] ARE NEGATIVE, THE
                     INTEGRATION  IS  PERFORMED  WITH  A CONSTANT STEP
                     DATA[5];
            DATA[8]: DATA[8]  SHOULD  BE  0  IF  ARK  IS  CALLED   FOR
                     A FIRST TIME;  FOR  CONTINUED  INTEGRATION
                     DATA[8] SHOULD NOT BE CHANGED;
            DATA[11], ..., DATA[10 + DATA[1]]: POLYNOMIAL COEFFICIENTS
                     (SEE REFERENCE [3]);
            AFTER  EACH  STEP THE FOLLOWING BY-PRODUCTS ARE DELIVERED:
            DATA[8]: THE   NUMBER  OF   INTEGRATION  STEPS  PERFORMED;
            DATA[9]: AN  ESTIMATION  OF THE  LOCAL  ERROR  LAST  MADE;
            DATA[10]: INFORMATIVE MESSAGES:
                     DATA[10] = 0: NO DIFFICULTIES;
                     DATA[10] = 1: MINIMAL  STEPLENGTH   EXCEEDS   THE
                             STEPLENGTH   PRESCRIBED   BY    STABILITY
                             THEORY, I.E. DATA[5] > DATA[3] / DATA[4];
                             (TERMINATION OF ARK);
                             DECREASE  MINIMAL  STEPLENGTH;
            IF NECESSARY, DATA[I],I = 4(1)7, CAN BE UPDATED(AFTER EACH
            STEP) BY MEANS OF PROCEDURE OUT;
    OUT:    <PROCEDURE IDENTIFIER>;
            THE HEADING OF THIS PROCEDURE READS:
            "PROCEDURE" OUT;
            AFTER  EACH  INTEGRATION STEP PERFORMED INFORMATION CAN BE
            OBTAINED  OR UPDATED BY THIS PROCEDURE, E.G. THE VALUES OF
            T, U[M0 : M] AND DATA[I], I = 4(1)10.

DATA AND RESULTS:

    FOR THE INDICES  M0  AND  M  THE  FOLLOWING  REMARKS  CAN BE MADE:
    WHEN THE METHOD OF LINES  IS  APPLIED  TO  HYPERBOLIC DIFFERENTIAL
    EQUATIONS THE NUMBER OF RELEVANT ORDINARY  DIFFERENTIAL  EQUATIONS
    DECREASES  DURING  THE  INTEGRATION  PROCESS;  IN   PROCEDURE  ARK
    THIS   MAY   BE   REALIZED   BY INTEGERS  M0  AND  M,  WHICH   ARE
    DEFINED AS FUNCTIONS OF THE NUMBER OF RIGHT HAND SIDE EVALUATIONS.
    A SELECTION OF POSSIBLE ENTRIES FOR ARRAY DATA (DEPENDENT ON THE
    KIND OF INITIAL VALUE PROBLEM) IS GIVEN IN REFERENCE [4],SECTION 8.

PROCEDURES USED:

    INIVEC = CP31010,
    MULVEC = CP31020,
    DUPVEC = CP31030,
    VECVEC = CP34010,
    ELMVEC = CP34020,
    DECSOL = CP34301.

REQUIRED CENTRAL MEMORY : CIRCA 75 + 2 * (M - M0) MEMORY PLACES.

METHOD AND PERFORMANCE:

    ARK IS  AN  IMPLEMENTATION  OF  LOW  ORDER  STABILIZED RUNGE KUTTA
    METHODS (SEE REFERENCE [1]);
    AUTOMATIC STEPSIZE CONTROL IS PROVIDED BUT STEP-REJECTION HAS BEEN
    EXCLUDED IN ORDER TO SAVE STORAGE;
    BECAUSE OF ITS LIMITED STORAGE REQUIREMENTS AND ADAPTIVE STABILITY
    FACILITIES THE METHOD IS WELL SUITED FOR THE SOLUTION  OF  INITIAL
    BOUNDARY  VALUE  PROBLEMS   FOR  PARTIAL  DIFFERENTIAL  EQUATIONS;
    NUMERICAL   RESULTS,   OBTAINED   WITH   A    SLIGHTLY   DIFFERENT
    IMPLEMENTATION CAN BE FOUND IN REFERENCE [2].

REFERENCES:

    [1]. P.J. VAN DER HOUWEN.
            STABILIZED RUNGE KUTTA METHOD WITH LIMITED
            STORAGE REQUIREMENTS.
            MATH. CENTR. REPORT  TW 124/71;

    [2]. P.A. BEENTJES.
            AN ALGOL 60 VERSION OF STABILIZED RUNGE KUTTA
            METHODS (DUTCH).
            MATH. CENTR. REPORT NR 23/72;

    [3]. P.J. VAN DER HOUWEN, J. KOK.
            NUMERICAL SOLUTION OF A MINIMAX PROBLEM.
            MATH. CENTR. REPORT TW 123/71;

    [4]. P.J. VAN DER HOUWEN ET AL.
            ONE STEP METHODS FOR LINEAR INITIAL VALUE PROBLEMS, I.I.I.
            NUMERICAL EXAMPLES,
            MATH. CENTR. REPORT TW 130/71.

EXAMPLE OF USE:

    THE VALUES OF

    1. Y(1) AND Y(2) OF THE INITIAL VALUE PROBLEM
       DY / DX = Y - 2 * X / Y,    Y(0) = 1

     AND

    2. U(.6, 0) OF THE CAUCHY PROBLEM (SEE REFERENCE [2]):
       DU / DT = .5 * DU / DX,  U(0, X) = EXP(-X * X)

    MAY BE OBTAINED BY THE FOLLOWING PROGRAM:

"BEGIN" "INTEGER" M0, M, I; "REAL" T, TE, DAT;
    "ARRAY" Y[1 : 1], U[-150 : 150], DATA[1 : 14];

    "PROCEDURE" DER1(T, V); "REAL" T; "ARRAY" V;
    V[1]:= V[1] - 2 * T / V[1];

    "PROCEDURE" DER2(T, V); "REAL" T; "ARRAY" V;
    "BEGIN" "INTEGER" J; "REAL" V1, V2, V3;
       V2:= V[M0]; M0:= M0 + 1; M:= M - 1; V3:= V[M0];
       "FOR" J:= M0 "STEP" 1 "UNTIL" M "DO"
       "BEGIN" V1:= V2; V2:= V3; V3:= V[J + 1];
           V[J]:= 250 * (V3 - V1) / 3
       "END"
    "END" DER2;

    "PROCEDURE" OUT1;
    "IF" T = TE "THEN"
    "BEGIN" "IF" T = 1 "THEN" OUTPUT(61, "("/, "(" PROBLEM 1")", //,
        "(" X  NUMBER OF INTEGRATION STEPS Y(COMPUTED)  Y(EXACT)")",
        //")");
        OUTPUT(61, "("ZD, 13ZD,12B,2(-3ZD.7D),"("...")", /")",
        T, DATA[8], Y[1], SQRT(2 * T + 1));
        TE:= 2
    "END" OUT1;

    "PROCEDURE" OUT2;
    "IF" T = .6 "THEN"
    OUTPUT(61, "("//, "(" PROBLEM 2")", //,
       "(" NUMBER OF DERIVATIVE CALLS")",
       "("  U(.6, 0)COMPUTED   U(.6, 0)EXACT")", //, 13ZD,
       2(-10Z.7D), "("...")"")", DATA[1] * DATA[8], U[0], EXP(-.09));

    I:= 1;
    "FOR" DAT:= 3, 3, 1, 1, "-3, "-6, "-6, 0, 0, 0, 1, .5, 1 / 6 "DO"
    "BEGIN" DATA[I]:= DAT; I:= I + 1 "END";
    T:= 0; Y[1]:= 1; TE:= 1;
    ARK(T, TE, 1, 1, Y, DER1, DATA, OUT1);
    I:= 1;
    "FOR" DAT:= 4, 3, SQRT(8), 500 / 3, DATA[3] / DATA[4], -1, -1,
       0, 0, 0, 1, .5, 1 / 6, 1 / 24 "DO"
    "BEGIN" DATA[I]:= DAT; I:= I + 1 "END";
    M0:= -150; M:= 150; T:= 0; U[0]:= 1;
    "FOR" I:= 1 "STEP" 1 "UNTIL" M "DO"
    U[I]:= U[-I]:= EXP(-(.003 * I) ** 2);
    ARK(T, .6, M0, M, U, DER2, DATA, OUT2)
"END"

    THIS PROGRAM DELIVERS:

 PROBLEM 1

 X  NUMBER OF INTEGRATION STEPS Y(COMPUTED)  Y(EXACT)

 1            38                1.7320535    1.7320508...
 2            56                2.2360928    2.2360680...

 PROBLEM 2

 NUMBER OF DERIVATIVE CALLS  U(.6, 0)COMPUTED   U(.6, 0)EXACT

           144               .9139326           .9139312...

SOURCE TEXT(S):
"CODE" 33061;