NUMAL Section 5.2.1.1.1.1.C

BEGIN SECTION : 5.2.1.1.1.1.C (February, 1979)

PROCEDURE : RK4A.

AUTHOR:J.A.ZONNEVELD.

CONTRIBUTORS:M.BAKKER AND I.BRINK.

INSTITUTE: MATHEMATICAL CENTRE,

RECEIVED: 730715.

BRIEF DESCRIPTION:

    RK4A SOLVES AN INITIAL VALUE PROBLEM FOR A SINGLE FIRST ORDER
    ORDINARY DIFFERENTIAL EQUATION   DY / DX = F(X,Y), WHERE F(X,Y)
    MAY BECOME LARGE, E.G. IN THE NEIGHBOURHOOD OF A SINGULARITY.
    RK4A INTERCHANGES THE DEPENDENT AND INDEPENDENT VARIABLE.
    THE EQUATION IS UPPOSED TO BE NON-STIFF.

KEYWORDS:

    INITIAL VALUE PROBLEM,
    SINGLE FIRST ORDER ORDINARY DIFFERENTIAL EQUATION
    LARGE DERIVATIVE VALUES.

CALLING SEQUENCE:

    THE HEADING OF THE PROCEDURE READS:
    "PROCEDURE" RK4A(X,XA,B,Y,YA,FXY,E,D,FI,XDIR,POS);
    "VALUE" FI,XDIR,POS;
    "BOOLEAN" FI,XDIR,POS;
    "REAL" X,XA,B,Y,YA,FXY;
    "ARRAY" E,D;
    "CODE" 33016;

    THE MEANING OF THE FORMAL PARAMETERS IS:

    X:    <VARIABLE>;
          THE INDEPENDENT VARIABLE ;
          UPON COMPLETION OF A CALL X IS EQUAL TO THE MOST RECENT
          VALUE OF THE INDEPENDENT VARIABLE;
    XA:   <ARITHMETIC EXPRESSION>;
          ENTRY: THE INITIAL VALUE OF X;
    B:    <ARITHMETIC EXPRESSION>;
          B DEPENDS ON X AND Y.
          THE EQUATION B=0 FULFILLED WITHIN THE TOLERANCES E[4] AND
          E[5]  SPECIFIES THE END OF THE INTEGRATION INTERVAL.
          AT THE END OF EACH INTEGRATION STEP B IS EVALUATED AND IS
          TESTED FOR CHANGE OF SIGN;
    Y:    <VARIABLE>;
          THE DEPENDENT VARIABLE;
    YA:   <ARITHMETIC EXPRESSION>;
          ENTRY: THE VALUE OF Y AT X=XA;
    FXY:  <ARITHMETIC EXPRESSION>;
          FXY GIVES THE VALUE OF DY/DX;
          FXY USES X AND Y AS JENSEN PARAMETERS.
    E:    <ARRAY IDENTIFIER>;
          "ARRAY" E[0:5];
          ENTRY:
          E[0], E[2] : RELATIVE TOLERANCES, FOR X AND Y RESPECTIVELY;
          E[1], E[3] : ABSOLUTE TOLERANCES, FOR X AND Y RESPECTIVELY;
          E[4], E[5] : TOLERANCES USED IN THE DETERMINATION OF
          THE ZERO OF B;
    D:    <ARRAY IDENTIFIER>;
          "ARRAY" D[0:4];
          AFTER COMPLETION OF EACH STEP WE HAVE :
          IF D[0]>0 THEN X IS THE INTEGRATION VARIABLE;
          IF D[0]<0 THEN Y IS THE INTEGRATION VARIABLE;
          D[1] IS THE NUMBER OF STEPS SKIPPED;
          D[2] IS THE STEPSIZE;
          D[3] IS EQUAL TO THE LAST VALUE OF X;
          D[4] IS EQUAL TO THE LAST VALUE OF Y;
    FI:   <BOOLEAN EXPRESSION>;
          IF FI="TRUE" THEN THE INTEGRATION IS STARTED WITH INITIAL
          VALUES X=XA,Y=YA.
          IF FI="FALSE" THEN THE INTEGRATION IS STARTED WITH X=D[3],
          Y=D[4];
    XDIR, POS : <BOOLEAN EXPRESSION>;
          IF FI="TRUE" THEN THE INTEGRATION STARTS IN SUCH A WAY THAT
          IF POS="TRUE" AND XDIR="TRUE" THEN X INCREASES,
          IF POS="TRUE" AND XDIR="FALSE" THEN Y INCREASES,
          IF POS="FALSE" AND XDIR="TRUE" THEN X DECREASES,
          IF POS="FALSE" AND XDIR="FALSE" THEN Y DECREASES.

PROCEDURES USED : ZEROIN = CP34150.

METHOD AND PERFORMANCE:

    RK4A IS BASED ON AN EXPLICIT, 5-TH ORDER RUNGE-KUTTA METHOD AND
    INTERCHANGES THE DEPENDENT AND INDEPENDENT INTEGRATION VARIABLE.
    IF  ABS(F(X,Y)) < 1, X IS USED AS INTEGRATION VARIABLE,OTHERWISE Y.
    THE PROCEDURE IS PROVIDED WITH STEP SIZE AND ERROR CONTROL.
    FOR DETAILS,E.G. HOW TO USE ARRAY E AND HOW TO SPECIFY THE ENDPOINT
    SEE [1] ( RK4A IS A SLIGHTLY ADAPTED VERSION OF RK4).

REFERENCES:
    [1]J.A.ZONNEVELD,
       AUTOMATIC NUMERICAL INTEGRATION.
       MATHEMATICAL CENTRE TRACT 8(1970).

EXAMPLE OF USE:

    THE SOLUTION OF THE DIFFERENTIAL EQUATION
    DY/DX=1-2*(X**2+Y), X>=0,
    Y=0               , X =0,
    IS REPRESENTED BY THE PARABOLA Y=X*(1-X);
    WE WOULD LIKE TO FIND THE VALUE OF X FOR WHICH THE CURVE
    OF THE SOLUTION INTERSECTS THE LINE Y+X=0.
    THE SOLUTION CAN BE OBTAINED BY THE FOLLOWING PROGRAM:

    "BEGIN" "COMMENT" INTEGRATION OF DY/DX=1-2*(Y+X**2), X>=0,
        Y(0)=0,UNTIL THE CONDITION Y+X=0 IS SATISFIED;
        "REAL" X,Y;"ARRAY" D[0:4],E[0:5];
        E[0]:=E[1]:=E[2]:=E[3]:=E[4]:=E[5]:="-6;
        RK4A(X,0.0,X+Y,Y,0.0,1-2*(X*X+Y),E,D,"TRUE","TRUE","TRUE");
        OUTPUT(61,"("10B"("X=")"+D.9D,10B"("EXACTLY:")"2B+D.9D/,
        10B"("Y=")"+D.9D/,10B"("Y-X*(1-X)=")"+.10D")",X,2,
        Y,Y-X*(1-X))
    "END"

    THE PROGRAM PRINTS THE FOLLOWING RESULTS:

    X=1.9999998554     EXACTLY: 2.0000000000

    Y=-1.9999995347427

    Y-X*(1-X)=0.0000000313

SOURCE TEXT(S):
"CODE" 33016;