NUMAL Section 5.2.1.1.1.1.I

BEGIN SECTION : 5.2.1.1.1.1.I (February, 1979)

PROCEDURE  : EFRK.

AUTHOR: K.DEKKER.

INSTITUTE: MATHEMATICAL CENTRE.

RECEIVED: 740710.

BRIEF DESCRIPTION:

    EFRK  SOLVES  AN INITIAL VALUE PROBLEM  FOR A SYSTEM OF FIRST
    ORDER ORDINARY DIFFERENTIAL EQUATIONS   DU / DT = H(T,U) .
    EFRK IS A SPECIAL PURPOSE PROCEDURE FOR STIFF EQUATIONS WITH A
    KNOWN, CLUSTERED EIGENVALUE SPECTRUM.

KEYWORDS:

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

CALLING SEQUENCE:

    THE HEADING OF THE PROCEDURE EFRK READS:
    "PROCEDURE" EFRK ( T,TE, M0,M, U, SIGMA, PHI, DIAMETER, DERIVATIVE,
                       K, STEP, R, L, BETA,  THIRDORDER, TOL, OUTPUT );
    "VALUE" R,L;
    "INTEGER" M0,M,K,R,L;
    "REAL" T,TE,SIGMA,PHI,DIAMETER,STEP,TOL;
    "ARRAY" U,BETA;
    "BOOLEAN" THIRDORDER;
    "PROCEDURE" DERIVATIVE,OUTPUT;
    "CODE" 33070;

    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 ETPRESSION>;
            THE FINAL VALUE OF T (TE>=T);
    M0:     <ARITHMETIC EXPRESSION>;
            THE INDEX OF THE FIRST EQUATION;
    M:      <ARITHMETIC EXPRESSION>;
            THE INDEX OF THE LAST EQUATION;

    U:      <ARRAY IDENTIFIER>;
            "REAL" "ARRAY" U[M0:M];
            THE DEPENDENT VARIABLE;
            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;
    SIGMA:  <ARITHMETIC EXPRESSION>;
            THE  MODULUS  OF  THE POINT AT WHICH EXPONENTIAL FITTING IS
            DESIRED , FOR EXAMPLE AN APPROXIMATION OF THE CENTRE OF THE
            LEFT HAND CLUSTER;
    PHI:    <ARITHMETIC EXPRESSION>;
            THE ARGUMENT OF THE CENTRE OF THE LEFT HAND CLUSTER; IN THE
            CASE OF  TWO COMPLEX CONJUGATED CLUSTERS , THE  ARGUMENT OF
            THE CENTRE IN THE SECOND QUADRANT SHOULD BE TAKEN;
    DIAMETER: <ARITHMETIC EXPRESSION>;
            THE DIAMETER OF THE LEFT HAND CLUSTER OF EIGENVALUES OF THE
            JACOBIAN MATRIX  OF THE SYSTEM  OF DIFFERENTIAL  EQUATIONS;
            IN CASE OF NON-LINEAR EQUATIONS DIAMETER SHOULD HAVE SUCH A
            VALUE THAT THE VARIATION OF THE EIGENVALUES IN THIS CLUSTER
            IN THE PERIOD ( T ,T+STEP ) IS LESS THAN HALF THE DIAMETER;
    DERIVATIVE: <PROCEDURE IDENTIFIER>;
            THE HEADING OF THIS PROCEDURE READS:
            "PROCEDURE" DERIVATIVE(T,U); "REAL" T; "ARRAY" U;
            THIS PROCEDURE SHOULD DELIVER THE VALUE OF H(T,U) IN THE
            POINT (T,U) IN THE ARRAY U;
    K:      <VARIABLE>;
            COUNTS THE NUMBER OF INTEGRATION STEPS TAKEN;
            FOR EXAMPLE, K MAY BE USED IN THE EXPRESSION FOR TE;
            ENTRY: AN (ARBIRARY) CHOSEN VALUE K0, E.G. K0=0;
            EXIT : K0 + THE NUMBER OF INTEGRATION STEPS PERFORMED;
    STEP:   <ARITHMETIC EXPRESSION>;
            THE  STEPSIZE  CHOSEN  WILL  BE  AT  MOST  EQUAL TO  STEP ;
            THIS  STEPSIZE  MAY BE  REDUCED  BY  STABILITY CONSTRAINTS,
            IMPOSED BY  A  POSITIVE DIAMETER , OR  BY CONSIDERATIONS OF
            INTERNAL STABILITY (SEE REF[1], PAGE 11);
    R:      <ARITHMETIC EXPRESSION>;
            R + L: THE NUMBER OF  EVALUATIONS OF  H(T, U)  ON WHICH THE
            RUNGE-KUTTA SCHEME IS BASED;
            FOR R=1,2,>=3 FIRST, SECOND AND THIRD ORDER ACCURACY MAY BE
            OBTAINED  BY  AN  APPROPRIATE  CHOICE  OF THE  ARRAY  BETA;
    L:      <ARITHMETIC EXPRESSION>;
            ENTRY:
            IF PHI = 4*ARCTAN(1): THE ORDER OF THE EXPONENTIAL FITTING,
            ELSE TWICE THE ORDER OF THE EXPONENTIAL FITTING;
            NOTE THAT L SHOULD BE EVEN IN THE LATTER CASE;
    BETA:   <ARRAY IDENTIFIER>;
            "REAL" "ARRAY" BETA[0:R+L];
            ENTRY: THE ELEMENTS  BETA[I] , I=0, ... ,R  SHOULD HAVE THE
                   VALUE OF THE R+1 FIRST COEFFICIENTS OF THE STABILITY
                   POLYNOMIAL;

    THIRDORDER: <BOOLEAN EXPRESSION>;
            IF THIRD ORDER ACCURACY IS DESIRED , THIRDORDER SHOULD HAVE
            THE VALUE "TRUE" , IN COMBINATION WITH  APPROPRIATE CHOICES
            OF R (R>=3) AND THE ARRAY BETA ( BETA[I]=1/I!, I=0,1,2,3 );
            IN ALL OTHER CASES THIRDORDER MUST  HAVE THE VALUE "FALSE";
    TOL:    <ARITHMETIC EXPRESSION>;
            AN UPPERBOUND FOR  THE ROUNDING ERRORS  IN THE COMPUTATIONS
            IN ONE RUNGE-KUTTA STEP ; IN SOME CASES ( E.G. LARGE VALUES
            OF SIGMA AND R ) TOL WILL CAUSE A DECREASE OF THE STEPSIZE;
    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 SOME PARAMETERS , FOR
            EXAMPLE T, K, U, AND COMPUTE NEW VALUES FOR SIGMA, PHI, AND
            DIAMETER.

PROCEDURES USED:

    INIVEC= CP31010,
    ELMVEC= CP34020,
    DEC   = CP34300,
    SOL   = CP34051.

REQUIRED CENTRAL MEMORY : CIRCA 30 + (M - M0) + L * (5 + L).

METHOD AND PERFORMANCE:
    EFRK IS BASED ON EXPLICIT RUNGE-KUTTA METHODS OF ORDER 1, 2 AND 3,
    WHICH MAKE USE OF EXPONENTIAL FITTING. AUTOMATIC ERROR CONTROL
    IS NOT PROVIDED.
    A DETAILED DESCRIPTION  OF THE METHOD  AND  SOME NUMERICAL EXAMPLES
    ARE GIVEN IN REF[1]. REF [3], PAGE 170 REPRESENTS A BRIEF SURVEY. A
    COMPARATIVE TEST  OVER A LARGE CLASS  OF  DIFFERENTIAL EQUATIONS IS
    GIVEN IN REF [4].
    FROM THESE RESULTS IT APPEARS  THAT CALLS WITH  THIRDORDER = "TRUE"
    ARE LESS ADVISABLE.

REFERENCES:

    [1]. K. DEKKER.
         AN  ALGOL 60  VERSION  OF  EXPONENTIALLY  FITTED   RUNGE-KUTTA
         METHODS (DUTCH).
         NR 25 (1972), MATHEMATICAL CENTRE.

    [2]. T. J. DEKKER, P. W. HEMKER AND P. J. 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.

    [4]. P. A. BEENTJES, K. DEKKER, H. C. HEMKER AND M. V. VELDHUIZEN.
         COLLOQUIUM STIFF DIFFERENTIAL EQUATIONS 3 (DUTCH).
         MC SYLLABUS 15.3, (1974) 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 487  AND  Y[2] = .433 710 353  5768.
    THE FOLLOWING PROGRAM  SHOWS SOME DIFFERENT CALLS  OF THE PROCEDURE
    EFRK:

    "BEGIN"

        "INTEGER" K,R,L,PASSES;
        "REAL" X,SIGMA,PHI,TIME,STEP,DIAMETER;
        "REAL" "ARRAY" Y[1:2],BETA[0:6];

        "PROCEDURE" DER(X,Y); "REAL" X; "ARRAY" Y;
        "BEGIN" "REAL" Y1,Y2; Y1:=Y[1]; Y2:=Y[2];
            Y[1]:=(Y1+.99)*(Y2-1)+.99;
            Y[2]:=1000*((1+Y1)*(1-Y2)-1);
            PASSES:=PASSES+1
        "END";

        "PROCEDURE" OUT;
        "BEGIN" "REAL" S;
            S:=(-1000*Y[1]-1001+Y[2])/2;
            SIGMA:=ABS(S-SQRT(S*S+10*(Y[2]-1)));
            DIAMETER:=2*STEP*ABS(1000*(1.99*Y[2]-2*Y[1]*(1-Y[2])));
            "IF" X=50 "THEN"
            OUTPUT(61,"("4BD,2BD,2(-5ZD),2(4B+.3DB3DB3D),-5ZD.3D,/")",
            R,L,K,PASSES,Y[1],Y[2],CLOCK-TIME)
        "END";

        OUTPUT(61,"(""("    THIS LINE AND THE FOLLOWING TEXT IS ")"
        "("PRINTED BY THIS PROGRAM")",//,
        "("    THE RESULTS WITH EFRK ARE:")",/,
        "("    R  L     K   DER.EV.        Y[1]             Y[2]")"
        "("         TIME")",/")");
        PHI:=4*ARCTAN(1); BETA[0]:=BETA[1]:=1;
        "FOR" R:=1,2,3 "DO" "FOR" L:=1,2,3 "DO"
        "BEGIN" "FOR" K:=2 "STEP" 1 "UNTIL" R "DO"
            BETA[K]:=BETA[K-1]/K;
            "FOR" STEP:=1,.1 "DO"
            "BEGIN" PASSES:=K:=0; X:=Y[2]:=0; Y[1]:=1; TIME:=CLOCK;
                OUT;
           EFRK(X,50.0,1,2,Y,SIGMA,PHI,DIAMETER,DER,K,STEP,R,L,BETA,
                R>=3,"-4,OUT);
            "END"; OUTPUT(61,"("/")");
        "END";
    "END"

    THIS LINE AND THE FOLLOWING TEXT IS PRINTED BY THIS PROGRAM:

    THE RESULTS WITH EFRK ARE:
    R  L     K   DER.EV.         Y[1]             Y[2]         TIME
    1  1    237    474     +.765 812 555    +.433 689 306      1.395
    1  1    501   1002     +.765 847 870    +.433 700 619      3.381

    1  2     52    156     +.765 570 874    +.433 615 119      0.465
    1  2    501   1503     +.765 848 220    +.433 700 709      4.200

    1  3     52    208     +.765 571 278    +.433 615 202      0.531
    1  3    500   2000     +.765 848 512    +.433 700 827      4.879

    2  1   3317   9951     +.765 878 320    +.433 710 353     21.808
    2  1   1050   3150     +.765 878 321    +.433 710 330      7.153

    2  2    174    696     +.765 878 335    +.433 710 335      1.385
    2  2    501   2004     +.765 878 323    +.433 709 211      4.915

    2  3     57    285     +.765 881 339    +.433 817 185      0.642
    2  3    501   2505     +.765 878 323    +.433 709 725      5.756

    3  1   7010  28040     +.765 878 320    +.433 710 354     55.298
    3  1   3255  13020     +.765 878 320    +.433 710 374     25.772

    3  2    949   4745     +.765 878 319    +.433 711 893      8.499
    3  2   1384   6920     +.765 862 498    +.449 724 830     13.452

    3  3    917   5502     +.765 878 018    +.434 105 184      9.143
    3  3   1166   6996     +.765 861 696    +.433 705 641     15.512

SOURCE TEXT(S):
"CODE" 33070;