NUMAL Section 5.2.1.1.1.2.A

BEGIN SECTION : 5.2.1.1.1.2.A (August, 1974)

AUTHOR: S.P.N. VAN KAMPEN.

INSTITUTE: MATHEMATICAL CENTRE.

RECEIVED: 730529.

BRIEF DESCRIPTION:

    EFSIRK SOLVES AN INITIAL VALUE PROBLEM, GIVEN AS AN AUTONOMOUS
    SYSTEM OF FIRST ORDER DIFFERENTIAL EQUATIONS DY/DX = F(Y), BY
    MEANS OF AN EXPONENTIALLY FITTED, SEMI-IMPLICIT RUNGE-KUTTA
    METHOD; IN PARTICULAR THIS PROCEDURE IS SUITABLE FOR THE
    INTEGRATION OF STIFF EQUATIONS.

KEYWORDS:

    DIFFERENTIAL EQUATIONS,
    INITIAL VALUE PROBLEM,
    AUTONOMOUS SYSTEM,
    STIFF EQUATIONS,
    SEMI-IMPLICIT RUNGE-KUTTA METHOD,
    EXPONENTIAL FITTING.

CALLING SEQUENCE:

    HEADING:
    "PROCEDURE" EFSIRK(X, XE, M, Y, DELTA, DERIVATIVE, JACOBIAN, J,
                       N, AETA, RETA, HMIN, HMAX, LINEAR, OUTPUT);
    "VALUE" M; "INTEGER" M, N;
    "REAL" X, XE, DELTA, AETA, RETA, HMIN, HMAX;
    "PROCEDURE" DERIVATIVE, JACOBIAN, OUTPUT;
    "ARRAY" Y, J;
    "BOOLEAN" LINEAR;
    "CODE" 33160;

    THE MEANING OF THE FORMAL PARAMETERS IS:
    X:      <VARIABLE>;
            THE INDEPENDENT VARIABLE X;
            ENTRY: THE INITIAL VALUE X0;
            EXIT : THE END VALUE XE;
    XE:     <ARITHMETIC EXPRESSION>;
            THE END VALUE OF X;
    M:      <ARITHMETIC EXPRESSION>;
            THE NUMBER OF DIFFERENTIAL EQUATIONS;
    Y:      <ARRAY IDENTIFIER>;
            "ARRAY" Y[1 : M];
            THE DEPENDENT VARIABLE;
            DURING THE INTEGRATION PROCESS THE COMPUTED SOLUTION
            AT THE POINT X IS ASSIGNED TO THE ARRAY Y;
            ENTRY: THE INITIAL VALUES OF THE SOLUTION OF THE SYSTEM;

    DELTA:  <ARITHMETIC EXPRESSION>;
            DELTA DENOTES THE REAL PART OF THE POINT AT WHICH
            EXPONENTIAL FITTING IS DESIRED;
            ALTERNATIVES:
            DELTA = (AN ESTIMATE OF) THE REAL PART OF THE, IN  ABSOLUTE
            VALUE, LARGEST EIGENVALUE OF THE JACOBIAN MATRIX OF THE
            SYSTEM;
            DELTA < -10**14, IN ORDER TO OBTAIN ASYMPTOTIC
            STABILITY;
            DELTA = 0, IN ORDER TO OBTAIN A HIGHER ORDER OF ACCURACY IN
            CASE OF LINEAR OR ALMOST LINEAR EQUATIONS;
    DERIVATIVE: <PROCEDURE IDENTIFIER>;
            "PROCEDURE" DERIVATIVE(A); "ARRAY" A;
            WHEN IN EFSIRK DERIVATIVE IS CALLED, A[I] CONTAINS THE
            VALUES OF Y[I];
            UPON COMPLETION OF A CALL OF DERIVATIVE, THE ARRAY A
            SHOULD CONTAIN THE VALUES OF F(Y);
            NOTE THAT THE VARIABLE X SHOULD NOT BE USED IN DERIVATIVE,
            BECAUSE THE DIFFERENTIAL EQUATION IS SUPPOSED TO BE
            AUTONOMOUS;
    JACOBIAN: <PROCEDURE IDENTIFIER>;
            "PROCEDURE" JACOBIAN(J, Y); "ARRAY" J, Y;
            WHEN IN EFSIRK JACOBIAN IS CALLED THE ARRAY Y CONTAINS
            THE VALUES OF THE DEPENDENT VARIABLE;
            UPON COMPLETION OF A CALL OF JACOBIAN THE ARRAY J SHOULD
            CONTAIN THE VALUES OF THE JACOBIAN MATRIX OF F(Y);
    J:      <ARRAY IDENTIFIER>;
            J[1 : M, 1 : M];
            J IS AN AUXILLIARY ARRAY WHICH IS USED IN THE PROCEDURE
            JACOBIAN;
    N:      <VARIABLE>;
            AN INTEGER WHICH COUNTS THE INTEGRATION STEPS;
    AETA, RETA:
            <ARITHMETIC EXPRESSION>;
            REQUIRED ABSOLUTE AND RELATIVE LOCAL ACCURACY;
    HMIN, HMAX:
            <ARITHMETIC EXPRESSION>;
            MINIMAL  AND MAXIMAL  STEPSIZE BY WHICH  THE INTEGRATION IS
            PERFORMED;
    LINEAR: <BOOLEAN EXPRESSION>;
            IF LINEAR = "TRUE" THE PROCEDURE JACOBIAN WILL ONLY BE
            CALLED IF N = 1; THE INTEGRATION WILL THEN BE PERFORMED
            WITH A STEPSIZE HMAX; THE CORRESPONDING REDUCTION
            OF COMPUTING TIME CAN BE EXPLOITED IN CASE OF LINEAR OR
            ALMOST LINEAR EQUATIONS;
    OUTPUT: <PROCEDURE IDENTIFIER>;
            "PROCEDURE" OUTPUT;
            IN OUTPUT ONE MAY PRINT THE VALUES OF E.G. X,
            Y[I], J[K, L] AND N.

DATA AND RESULTS: SEE REF[2] AND [3].

PROCEDURES USED:

    VECVEC = CP34010,
    MATVEC = CP34011,
    MATMAT = CP34013,
    GSSELM = CP34231,
    SOLELM = CP34061.

REQUIRED CENTRAL MEMORY:

    EXECUTION FIELD LENGTH: CIRCA M * M + 5 * M.

RUNNING TIME:

    DEPENDS  STRONGLY  ON  THE  DIFFERENTIAL  EQUATION  TO  BE   SOLVED

LANGUAGE: ALGOL 60.

METHOD AND PERFORMANCE:

    THE PROCEDURE EFSIRK IS AN EXPONENTIALLY FITTED, A-STABLE,
    SEMI-IMPLICIT RUNGE-KUTTA METHOD OF THIRD ORDER (SEE REF[1]
    AND [2]). THE ALGORITHM USES FOR EACH STEP TWO FUNCTION EVALUATIONS
    AND IF LINEAR = "FALSE" ONE EVALUATION OF THE JACOBIAN MATRIX.
    THE STEPSIZE IS NOT DETERMINED BY THE ACCURACY OF THE NUMERICAL
    SOLUTION, BUT BY THE AMOUNT BY WHICH THE GIVEN DIFFERENTIAL
    EQUATION DIFFERS FROM A LINEAR EQUATION (SEE REF[2]).
    THE PROCEDURE DOES NOT REJECT INTEGRATION STEPS.

REFERENCES:

    [1].P.J. VAN DER HOUWEN.
        ONE-STEP METHODS WITH ADAPTIVE STABILITY FUCNTIONS
        FOR THE INTEGRATION OF DIFFERENTIAL EQUATIONS.
        LECTURE NOTES OF THE CONFERENCE ON
        NUMERISCHE, INSBESONDERE APPROXIMATIONSTHEORETISCHE
        BEHANDLUNG VON FUNKTIONALGLEICHUNGEN.
        OBERWOLFACH, DECEMBER, 3 - 12, 1972.

    [2].SYLLABUS COLLOQUIUM STIFF DIFFERENTIAL EQUATIONS 2 (DUTCH).
        MATH.CENTR. SYLLABUS 15.2/73.

    [3].SYLLABUS COLLOQUIUM STIFF DIFFERENTIAL EQUATIONS 3 (DUTCH).
        MATH.CENTR. SYLLABUS 15.3/73.
        TO APPEAR IN 1973.

EXAMPLE OF USE:

    WE CONSIDER THE DIFFERENTIAL EQUATION
    DY / DX = -EXP(X) * (Y - LN(X)) + 1 / X,
    ON THE INTERVAL  [0.01, 8], WITH  INITIAL VALUE  Y(0.01) = LN(0.01)
    AND ANALYTICAL SOLUTION Y(X) = LN(X);
    FOR THE FIT POINT WE USE THE  EIGENVALUE  OF THE JACOBIAN MATRIX,
    I.E. DELTA = -EXP(X);

    "BEGIN"
        "PROCEDURE" DER(Y); "ARRAY" Y;
        "BEGIN" "REAL" Y2; Y2:= Y[2];
            DELTA:= -EXP(Y2); LNX:= LN(Y2);
            Y[1]:= (Y[1] - LNX) * DELTA + 1 / Y2;
            Y[2]:= 1
        "END" DER;
        "PROCEDURE" JAC(J, Y); "ARRAY" J, Y;
        "BEGIN" "REAL" Y2; Y2:= Y[2];
            J[1, 1]:= DELTA;
            J[1, 2]:= (Y[1] - LNX - 1 / Y2) * DELTA - 1 / (Y2 * Y2);
            J[2, 1]:= J[2, 2]:= 0
        "END" JAC;
        "PROCEDURE" OUTP;
        "IF" X = XE "THEN"
        "BEGIN" "REAL" Y1; Y1:= Y[1]; LNX:= LN(X);
            OUTPUT(61, "(""("N = ")", 2ZD,
                       "("   X = ")", +D.D,
                       "("  Y(X) = ")", +D.5D,
                       "("  DELTA = ")", +3ZD.2D, /,
                       "("ABS. ERR. = ")", .2D"+2D,
                       "("  REL. ERR. = ")", .2D"+2D, //")",
                       N, X, Y1, DELTA,
                       ABS(Y1 - LNX), ABS((Y1 - LNX) / LNX));
            "IF" X = 0.4 "THEN" XE:= 8
        "END" OUTP;
        "INTEGER" N;
        "REAL" X, XE, DELTA, LNX;
        "ARRAY" Y[1 : 2], J[1 : 2, 1 : 2];

        XE:= 0.4; X:= 0.01; Y[1]:= LN(0.01); Y[2]:= X;
        EFSIRK(X, XE, 2, Y, DELTA, DER, JAC, J,
               N, "-2, "-2, 0.005, 1.5, "FALSE", OUTP)
    "END"

    THIS PROGRAM DELIVERS:

        N =  10   X = +0.4  Y(X) = -0.91099  DELTA =    -1.44
        ABS. ERR. = .53"-02  REL. ERR. = .58"-02

        N = 98    X = +8.0  Y(X) = +2.07911  DELTA = -2980.02
        ABS. ERR. = .33"-03  REL. ERR. = .16"-03

SOURCE TEXT(S):
"CODE" 33160;