NUMAL Section 5.2.1.1.2.1.C

BEGIN SECTION : 5.2.1.1.2.1.C (February, 1979)

PROCEDURE : RK3

AUTHOR:J.A.ZONNEVELD.

CONTRIBUTORS: M.BAKKER AND I.BRINK.

INSTITUTE: MATHEMATICAL CENTRE.

RECEIVED: 730715.

BRIEF DESCRIPTION:

    RK3 INTEGRATES THE SCALAR INITIAL VALUE PROBLEM
    (D/DX) (D/DX) Y = F(X,Y)   (WITHOUT THE DERIVATIVE (D/DX) Y IN F),
    A <= X <= B OR B <= X <= A,  Y(A) AND (D/DX) Y(A) PRESCRIBED.

KEYWORDS:

    INITIAL VALUE PROBLEM,
    SECOND ORDER DIFFERENTIAL EQUATION.

CALLING SEQUENCE:

    THE HEADING OF THE PROCEDURE READS:
    "PROCEDURE" RK3(X,A,B,Y,YA,Z,ZA,FXY,E,D,FI);
    "VALUE" B,FI;
    "REAL" X,A,B,Y,YA,Z,ZA,FXY;
    "BOOLEAN" FI;
    "ARRAY" E,D;
    "CODE" 33014;

    THE MEANING OF THE FORMAL PARAMETERS IS:
    X:    <VARIABLE>;
          THE INDEPENDENT VARIABLE.
          UPON COMPLETION OF A CALL OF RK3 ,
          IT IS EQUAL TO B;
    A:    <ARITHMETIC EXPRESSION>;
          THE STARTING VALUE OF X;
    B:    <ARITHMETIC EXPRESSION>;
          A VALUE PARAMETER, GIVING THE END VALUE OF X;
          B <= A IS ALLOWED;
    Y:    <VARIABLE>;
          THE DEPENDENT VARIABLE;
          EXIT : THE VALUE OF Y(X) AT X = B;
    YA:   <ARITHMETIC EXPRESSION>;
          ENTRY : THE VALUE OF Y AT X=A;
    Z:    <VARIABLE>;
          THE DERIVATIVE  DY/DX;
          EXIT : THE VALUE OF DY/DX AT X = B;
    ZA:   <ARITHMETIC EXPRESSION>;
          ENTRY : THE VALUE OF DY/DX AT X=A;
    FXY:  <ARITHMETIC EXPRESSION>;
          AN EXPRESSION,DEPENDING ON X AND Y ,GIVING THE VALUE OF
          (D/DX)(D/DX)Y;
    E:   <ARRAY IDENTIFIER>;
          "ARRAY" E[1:4];
          E[1] AND E[3] ARE USED AS RELATIVE TOLERANCES,
          E[2] AND E[4] ARE USED AS ABSOLUTE TOLERANCES
          FOR Y AND DY/DX, RESPECTIVELY;
    D:    <ARRAY IDENTIFIER>;
          "ARRAY" D[1:5];
          EXIT:
          ENTIER(D[1]+.5) IS THE NUMBER OF STEPS SKIPPED;
          D[2] IS THE LAST STEP LENGTH USED;
          D[3] IS EQUAL TO B;
          D[4] IS EQUAL TO Y(B);
          D[5] IS EQUAL TO DY/DX FOR X=B;
    FI:   <BOOLEAN EXPRESSION>;
          IF FI="TRUE" THEN THE INTEGRATION STARTS AT X=A WITH A TRIAL
          STEP B-A;IF FI="FALSE" THEN THE INTEGRATION IS CONTINUED
          VIZ.  WITH THE INITIAL CONDITIONS X=D[3], Y=D[4], Z=D[5] AND
          STEP LENGTH  H=D[2]*SIGN(B-D[3]); A,YA,ZA ARE IGNORED.

PROCEDURES USED : NONE.

METHOD AND PERFORMANCE :
    RK3 INTEGRATES (D/DX)(D/DX)Y = F(X,Y) FROM X TO B,WITH IF FI="TRUE"
    THEN X=A, Y=YA,DY/DX=ZA ELSE X=D[3], Y=D[4], Z=D[5].
    A 5-TH ORDER RUNGE-KUTTA METHOD IS USED.
    UPON COMPLETION OF A CALL OF RK3 WE HAVE  X=D[3]=B, Y=D[4]=Y[B],
    Z=D[5], I.E. THE VALUE OF DY/DX FOR X=B.
    RK3 USES AS ITS MINIMAL ABSOLUTE STEP LENGTH
    HMIN=MIN (E[2*J-1]*INT+E[2*J]) WITH 1<=J<=2 AND INT=
    ABS(B-("IF" FI "THEN" A "ELSE" D[3])).
    IF A STEP OF LENGTH ABS(H)<=HMIN IS REJECTED ,A STEP SIGN(H)*HMIN
    IS SKIPPED.  A STEP IS REJECTED IF THE ABSOLUTE VALUE OF THE LAST
    TERM TAKEN INTO ACCOUNT IS GREATER THEN (ABS(DY/DX)*E[1]+E[2])*
    ABS(H)/INT OR IF THAT TERM IS GREATER THEN (ABS(FXY)*E[3]+E[4])*
    ABS(H)/INT ( INT = ABS(B - A) ).
    SEE REF[1].

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

EXAMPLE OF USE:

    "BEGIN" "COMMENT" SOLUTION OF Y"=X*Y,Y(0)=0,Y'(0)=1;

    "REAL" "PROCEDURE" YEXACT(X);"VALUE" X;"REAL" X;
    "BEGIN" "INTEGER" N;"REAL" X3,S,TERM;
        X3:=X**3;TERM:=X;S:=0;
        "FOR" N:=3,N+3 "WHILE" ABS(TERM)>"-14 "DO"
        "BEGIN" S:=S+TERM;TERM:=TERM*X3/N/(N+1)
        "END";
        YEXACT:=S
    "END";

    "REAL" X,B,Y,Z;"BOOLEAN" FI;"ARRAY" D,E[1:5];
    E[1]:=E[3]:="-8;E[2]:=E[4]:="-12;
    "FOR" B:=.25,.50,.75,1.00 "DO"
    "BEGIN" FI:=B<.30;
        RK3(X,0.0,B,Y,0.0,Z,1.0,X*Y,E,D,FI);
        OUTPUT(61,"("10B"("Y-YEXACT=")".10D,5B"("X=")"Z.2D,
        5B"("Y=")"2D.10D//")",Y-YEXACT(X),X,Y)
    "END"
    "END"

    DELIVERS:

    Y-YEXACT=0.0000000000     X= .25     Y=00.2503256420

    Y-YEXACT=0.0000000000     X= .50     Y=00.5052238559

    Y-YEXACT=0.0000000000     X= .75     Y=00.7766332813

    Y-YEXACT=0.0000000000     X=1.00     Y=01.0853396481

SOURCE TEXT(S):
"CODE" 33014;