NUMAL Section 5.2.1.1.1.1.D

BEGIN SECTION : 5.2.1.1.1.1.D (February, 1979)

PROCEDURE : RK4NA.

AUTHOR:J.A.ZONNEVELD.

CONTRIBUTORS:M.BAKKER AND I.BRINK.

INSTITUTE: MATHEMATICAL CENTRE.

RECEIVED: 730715.

BRIEF DESCRIPTION:

    RK4NA SOLVES AN INITIAL VALUE PROBLEM FOR A SYSTEM OF FIRST
    ORDER ORDINARY DIFFERENTIAL EQUATIONS  DY / DX = F(X,Y), OF WHICH
    THE DERIVATIVE COMPONENTS ARE SUPPOSED TO BECOME LARGE, E.G. IN
    THE NEIGHBOURHOOD OF SINGULARITIES. RK4NA INTERCHANGES THE
    INDEPENDENT VARIABLE AND ONE DEPENDENT VARIABLE. THE SYSTEM
    IS SUPPOSED TO BE NON-STIFF.

KEYWORDS:

    INITIAL VALUE PROBLEM,
    SYSTEM OF FIRST ORDER ORDINARY DIFFERENTIAL EQUATIONS,
    LARGE DERIVATIVE COMPONENTS.

CALLING SEQUENCE:

    THE HEADING OF THE PROCEDURE READS:
    "PROCEDURE" RK4NA(X,XA,B,FXJ,J,E,D,FI,N,L,POS);
    "VALUE" FI,N,L,POS;
    "INTEGER" J,N,L;
    "BOOLEAN" FI,POS;
    "REAL" B,FXJ;
    "ARRAY" X,XA,E,D;
    "CODE" 33017;

    THE MEANING OF THE FORMAL PARAMETERS IS:

    X:    <ARRAY IDENTIFIER>;
          "ARRAY" X[0:N];
          X[0] IS THE INDEPENDENT VARIABLE,
          X[1],...,X[N] ARE THE DEPENDENT VARIABLES;
          EXIT:THE SOLUTION AT B=0;
    XA:   <ARRAY IDENTIFIER>;
          "ARRAY" XA[0:N];
          ENTRY: THE INITIAL VALUES OF X[0],...,X[N];
    B:    <ARITHMETIC EXPRESSION>;
          B DEPENDS ON X[0],...,X[N];
          IF THE EQUATION B=0 IS SATISFIED WITHIN A
          CERTAIN TOLERANCE,THE INTEGRATION IS TERMINATED;
          B IS EVALUATED AND TESTED FOR CHANGE OF SIGN AT THE
          END OF EACH STEP;
          FOR THE TOLERANCE SEE PARAMETER E.
    FXJ:  <ARITHMETIC EXPRESSION>;
          FXJ DEPENDS ON X[0],...,X[N] AND J, DEFINING THE RIGHT
          HAND SIDE OF THE DIFFERENTIAL EQUATION;
          AT EACH CALL IT DELIVERS :DX[J]/DX[0];
    J:    <VARIABLE>;
          J IS USED AS A JENSEN PARAMETER FOR FXJ;
    E:    <ARRAY IDENTIFIER>;
          "ARRAY" E[0:2*N+3];
          ENTRY: E[2*J] AND E[2*J+1] , 0<=J<=N ,
          ARE THE RELATIVE AND THE ABSOLUTE TOLERANCE ,
          RESPECTIVELY, ASSOCIATED WITH X[J];
          E[2*N+2] AND E[2*N+3] ARE THE RELATIVE AND ABSOLUTE
          TOLERANCE  USED IN THE DETERMINATION OF THE ZERO OF B;
    D:    <ARRAY IDENTIFIER>;
          "ARRAY" D[0:N+3];
          AFTER COMPLETION OF EACH STEP WE HAVE :
          ENTIER(D[0]+.5) IS THE NUMBER OF STEPS SKIPPED;
          D[2] IS THE STEP LENGTH;
          D[J+3] IS THE LAST VALUE OF X[J], J=0,...,J=N;

    FI:   <BOOLEAN EXPRESSION>
          IF FI="TRUE" THEN THE INTEGRATION IS STARTED WITH INITIAL
          CONDITIONS X[J]=XA[J];IF FI="FALSE" THEN THE INTEGRATION IS
          CONTINUED WITH X[J]=D[J+3];
    N:    <ARITHMETIC EXPRESSION>;
          THE NUMBER OF EQUATIONS;
    L:    <ARITHMETIC EXPRESSION>;
          AN INTEGER TO BE SUPPLIED BY THE USER, 0<=L<=N (SEE POS);
    POS:  <BOOLEAN EXPRESSION>
          IF FI="TRUE" THEN THE INTEGRATION STARTS IN SUCH A WAY
          THAT X[L] INCREASES IF POS="TRUE" AND X[L] DECREASES IF
          POS="FALSE";
          IF FI="FALSE" THEN POS IS OF NO SIGNIFICANCE.

PROCEDURES USED : ZEROIN = CP34150.

REQUIRED CENTRAL MEMORY : CIRCA 9*N MEMORY PLACES.

METHOD AND PERFORMANCE :
    RK4NA IS BASED ON AN EXPLICIT, 5-TH ORDER RUNGE-KUTTA METHOD
    AND INTERCHANGES VARIABLES. THE PROCEDURE IS PROVIDED WITH STEPSIZE
    AND ERROR CONTROL. FOR DETAILS, E.G. HOW TO USE ARRAY E, HOW TO
    SPECIFY THE ENDPOINT, HOW TO USE L AND POS, SEE [1] ( RK4NA IS A
    SLIGHTLY ADAPTED VERSION OF RK4N ).

REFERENCES:

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

EXAMPLE OF USE:

    THE PERIOD OF THE SOLUTION OF THE VAN DER POL EQUATION
    DX[1]/DT = X[2]
    DX[2]/DT = 10*(1-X[1]**2)*X[2]-X[1], T>=0,
    CAN BE OBTAINED BY THE FOLLOWING PROGRAM:

    "BEGIN" "COMMENT" VAN DER POL;
        "REAL" X0;

        "PROCEDURE" PRINT(X);"ARRAY" X;
        OUTPUT(61,"("/,4(+ZD.10D3B)")",X[0],X[1],X[2],X0);

        "INTEGER" J,K;"BOOLEAN" FIRST;
        "ARRAY" E[0:7],XA,X[0:2],D[0:5];
        "FOR" K:=0,1,2,3,4,5 "DO" E[K]:=.1"-6; E[6]:=E[7]:="-8 ;
        OUTPUT(61,"(""("VAN DER POL")",/BB"("EPS=")"D.10D,/BB
        "("THE VALUES OF X[0],X[1],X[2],P,RESPECTIVELY")"")",E[0]);
        X0:=XA[0]:=XA[2]:=0;XA[1]:=2;J:=0;PRINT(XA);
        FIRST:="TRUE";
        "FOR" J:=1,2,3,4 "DO"
        "BEGIN" RK4NA(X,XA,X[2],"IF" K=1 "THEN" X[2] "ELSE"
            10*(1-X[1]**2)*X[2]-X[1],K,E,D,FIRST,2,0,"TRUE");
            X0:=X[0]-X0;PRINT(X);FIRST:="FALSE"; X0:=X[0]
        "END"
    "END"

    THE PROGRAM PRINTS THE FOLLOWING RESULTS:

    VAN DER POL

    EPS=0.0000001000

    THE VALUES OF X[0],X[1],X[2],P,RESPECTIVELY:

    +0.00000000  +2.00000000  +0.00000000  +0.00000000

    +9.32386570  -2.01428560  +0.00000000  +9.32386570

    +18.86305411 +2.01428557  +0.00000000  +9.53918840

    +28.40224194  -2.01428858  -0.00000000  +9.53918783

    +37.94143003  +2.01428558  +0.00000000  +9.53918809

SOURCE TEXT(S):
"CODE" 33017;