NUMAL Section 5.2.1.1.1.2.B

BEGIN SECTION : 5.2.1.1.1.2.B (August, 1974)

AUTHOR: K.DEKKER.

INSTITUTE: MATHEMATICAL CENTRE.

RECEIVED: 1973/07/31.

BRIEF DESCRIPTION:

    EFERK SOLVES INITIAL VALUE PROBLEMS , GIVEN AS AN AUTONOMOUS SYSTEM
    OF FIRST ORDER DIFFERENTIAL EQUATIONS, BY MEANS OF AN EXPONENTIALLY
    FITTED, EXPLICIT RUNGE KUTTA METHOD OF THIRD ORDER, WHICH INVOLVES
    THE USE OF THE JACOBIAN MATRIX. AUTOMATIC STEP 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,
    EXPLICIT RUNGE KUTTA METHODS.

CALLING SEQUENCE:

    THE HEADING OF THE PROCEDURE EFERK READS:
    "PROCEDURE" EFERK(X,XE,M,Y,SIGMA,PHI,DERIVATIVE,J,JACOBIAN,
                  K,L,AUT,AETA,RETA,HMIN,HMAX,LINEAR,OUTPUT);
    "VALUE" L;
    "INTEGER" M,K,L;
    "REAL" X,XE,SIGMA,PHI,AETA,RETA,HMIN,HMAX;
    "ARRAY" Y,J;
    "BOOLEAN" AUT,LINEAR;
    "PROCEDURE" DERIVATIVE,JACOBIAN,OUTPUT;
    "CODE" 33120;

    THE MEANING OF THE FORMAL PARAMETERS IS:
    X:      <VARIABLE>;
            THE INDEPENDENT VARIABLE X;
            CAN BE USED IN DERIVATIVE, JACOBIAN, OUTPUT, ETC.;
            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>;
            "REAL" "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 MATRIX OF THE SYSTEM OF DIFFERENTIAL EQUATIONS;
    PHI:    <ARITHMETIC EXPRESSION>;
            THE  ARGUMENT  OF THE  COMPLEX POINT  AT WHICH  EXPONENTIAL
            FITTING IS DESIRED;
    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>;
            "REAL" "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  THE JACOBIAN AT THE POINT (Y)  HAS TO BE
            ASSIGNED TO THE ARRAY J;
    K:      <VARIABLE>;
            COUNTS THE NUMBER OF INTEGRATION STEPS TAKEN;
            FOR EXAMPLE, MAY BE USED IN THE EXPRESSION FOR XE;
    L:      <ARITHMETIC EXPRESSION>;
            ENTRY:
            IF PHI = 4*ARCTAN(1): THE ORDER OF THE EXPONENTIAL FITTING,
            ELSE TWICE THE ORDER OF THE EXPONENTIAL FITTING;
    AUT:    <BOOLEAN EXPRESSION>;
            IF THE SYSTEM HAS BEEN WRITTEN IN AUTONOMOUS FORM BY ADDING
            THE EQUATION DY[M]/DX = 1 TO THE SYSTEM , THEN AUT MAY HAVE
            THE VALUE "FALSE", ELSE AUT SHOULD HAVE THE VALUE "TRUE";
    AETA:   <ARITHMETIC EXPRESSION>;
            REQUIRED  ABSOLUTE  PRECISION  IN THE  INTEGRATION PROCESS;
            AETA HAS TO BE POSITIVE;
    RETA:   <ARITHMETIC EXPRESSION>;
            REQUIRED  RELATIVE  PRECISION  IN THE  INTEGRATION PROCESS;
            RETA HAS TO BE POSITIVE;
    HMIN:   <ARITHMETIC EXPRESSION>;
            THE STEPLENGTH CHOSEN WILL BE AT LEAST EQUAL TO HMIN;
    HMAX:   <ARITHMETIC EXPRESSION>;
            THE STEPLENGTH CHOSEN WILL BE AT MOST EQUAL TO HMAX;
    LINEAR: <ARITHMETIC EXPRESSION>;
            THE PROCEDURE JACOBIAN IS CALLED  ONLY IF LINEAR="FALSE" OR
            K=0; SO IF THE SYSTEM IS LINEAR , LINEAR MAY HAVE THE VALUE
            "TRUE";
    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 X, K, Y.

DATA AND RESULTS: SEE EXAMPLE OF USE, AND REF[4].

PROCEDURES USED:

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

REQUIRED CENTRAL MEMORY:
    EXECUTION FIELD LENGTH: CIRCA 30 + 4 * M + L * (5+L).

RUNNING TIME: DEPENDS  STRONGLY  ON THE DIFFERENTIAL EQUATION TO SOLVE.

LANGUAGE: ALGOL 60.

METHOD AND PERFORMANCE:

   THE PROCEDURE EFERK IS AN EXPONENTIALLY FITTED, SEMI-EXPLICIT RUNGE
   KUTTA METHOD OF THIRD ORDER ( SEE REF [1] AND [3] ) . THE ALGORITHM
   USES FOR EACH STEP TWO FUNCTION EVALUATIONS AND IF LINEAR = "FALSE"
   ONE EVALUATION OF THE JACOBIAN MATRIX.THE STEPSIZE IS DETERMINED BY
   AN ESTIMATION  OF THE LOCAL TRUNCATION ERROR  BASED ON THE RESIDUAL
   FUNCTION (SEE REF[3]). INTEGRATION STEPS ARE NOT REJECTED.

REFERENCES:

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

   [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]. (TO APPEAR).
         COLLOQUIUM STIFF DIFFERENTIAL EQUATIONS 3 (DUTCH).
         MC SYLLABUS 15.3, (1973) 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
    EFERK,AND ILLUSTRATES THE ACCURACIES WHICH MAY BE OBTAINED BY THEM:

    "BEGIN"

        "INTEGER" K,PASSES,PASJAC;
        "REAL" X,SIGMA,PHI,TIME,TOL;
        "REAL" "ARRAY" Y[1:2],J[1:2,1:2];

        "PROCEDURE" DER(Y); "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" 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;
            PASJAC:=PASJAC+1
        "END" JACOBIAN;

        "PROCEDURE" OUT;
        "IF" X=50 "THEN"
       OUTPUT(61,"("3(-5ZD),2(4B+.3DB3DB3D),-5ZD.3D,/")",K,PASSES,
        PASJAC,Y[1],Y[2],CLOCK-TIME);

        OUTPUT(61,"(""("    THIS LINE AND THE FOLLOWING TEXT IS ")"
        "("PRINTED BY THIS PROGRAM")",//,
        "("    THE RESULTS WITH EFERK -FIRST ORDER FIT- ARE:")",/,
        "("     K   DER.EV. JAC.EV.       Y[1]             Y[2]")"
        "("         TIME")",/")");
        PHI:=4*ARCTAN(1);
        "FOR" TOL:=1,"-1,"-2,"-3 "DO"
        "BEGIN" PASSES:=PASJAC:=0; X:=Y[2]:=0; Y[1]:=1; TIME:=CLOCK;
          EFERK(X,50.0,2,Y,SIGMA,PHI,DER,J,JACOBIAN,K,1,"TRUE",TOL,
                                         TOL,"-6,50.0,"FALSE",OUT);
        "END";
    "END"

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

    THE RESULTS WITH EFERK -FIRST ORDER FIT- ARE:
     K   DER.EV. JAC.EV.       Y[1]             Y[2]         TIME
    93    186     93    +.765 883 211    +.428 752 781      1.170
   105    210    105    +.765 878 445    +.433 569 561      1.296
   147    294    147    +.765 878 317    +.433 708 489      1.834
   266    532    266    +.765 878 320    +.433 710 229      3.297

SOURCE TEXT(S):
"CODE" 33120;