NUMAL Section 5.2.1.1.2.1.D

BEGIN SECTION : 5.2.1.1.2.1.D (February, 1979)

PROCEDURE : RK3N.

AUTHOR:J.A.ZONNEVELD.

CONTRIBUTORS: M.BAKKER AND I.BRINK.

INSTITUTE:MATHEMATICAL CENTRE.

RECEIVED: 730715.

BRIEF DESCRIPTION:

    RK3N INTEGRATES THE VECTOR INITIAL VALUE  PROBLEM
    (D/DX) (D/DX) Y = F(X,Y), A <= X <= B OR B <= X <= A,
    Y[J] (A)  AND  (D/DX) Y[J] (A)  PRESCRIBED.

KEYWORDS:

    INITIAL VALUE PROBLEM,
    SECOND ORDER DIFFERENTIAL EQUATION.

CALLING SEQUENCE:

    THE HEADING OF THE PROCEDURE READS:
    "PROCEDURE" RK3N(X,A,B,Y,YA,Z,ZA,FXYJ,J,E,D,FI,N);
    "VALUE" B,FI,N;
    "INTEGER" J,N;
    "REAL" X,A,B,FXYJ;
    "BOOLEAN" FI;
    "ARRAY" Y,YA,Z,ZA,E,D;
    "CODE" 33015;

    THE MEANING OF THE FORMAL PARAMETERS IS:
    X:    <VARIABLE>;
          THE INDEPENDENT VARIABLE.
          UPON COMPLETION OF A CALL OF RK3N,
          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:    <ARRAY IDENTIFIER>;
          "ARRAY" Y[1:N];
          THE VECTOR OF DEPENDENT VARIABLES;
          EXIT : THE VALUE OF Y[J](X) AT X = B;
    YA:   <ARRAY IDENTIFIER>;
          "ARRAY" YA[1:N];
          ENTRY : THE STARTING VALUES OF Y[J],I.E. THE VALUES AT X=A;
    Z:    <ARRAY IDENTIFIER>;
          "ARRAY" Z[1:N];
          THE DERIVATIVES OF THE DEPENDENT VARIABLES, Z[J] = DY[J]/DX;
          EXIT : THE VALUE OF Z[J](X) AT X = B;
    ZA:   <ARRAY IDENTIFIER>;
          "ARRAY" ZA[1:N];
          ENTRY : THE STARTING VALUES OF Z[J],I.E. THE VALUES AT X=A;
    FXYJ: <ARITHMETIC EXPRESSION>;
          AN EXPRESSION DEPENDING ON X,Y[1],...,Y[N],J,
          GIVING THE VALUE OF (D/DX)(D/DX)Y[J];
    J:    <VARIABLE>;
          A VARIABLE OF TYPE INTEGER,USED IN THE ACTUAL PARAMETER
          CORRESPONDING TO FXYJ,TO DENOTE THE NUMBER OF THE EQUATION
          REQUIRED (JENSEN'S DEVICE);
    E:    <ARRAY IDENTIFIER>;
          "ARRAY" E[1:4*N];
          THE ELEMENT E[2*J-1] IS A RELATIVE AND E[2*J] IS AN ABSOLUTE
          TOLERANCE ASSOCIATED WITH Y[J];
          E[2*(N+J)-1] IS A RELATIVE AND E[2*(N+J)] IS AN ABSOLUTE
          TOLERANCE ASSOCIATED WITH Z[J];
    D:    <ARRAY IDENTIFIER>;
          "ARRAY" D[1:2*N+3];
          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],...,D[N+3] ARE EQUAL TO Y[1],...,Y[N] FOR X=B;
          D[N+4],...,D[2*N+3] ARE EQUAL TO THE DERIVATIVES
          Z[1],...,Z[N] FOR X=B;
    FI:   <BOOLEAN EXPRESSION>;
          IF FI="TRUE" THEN THE INTEGRATION STARTS AT A ,WITH A TRIAL
          STEP B-A;IF FI="FALSE" THEN THE INTEGRATION IS CONTINUED VIZ.
          WITH THE INITIAL CONDITIONS:X=D[3],Y[J]=D[J+3],Z[J]=D[N+J+3],
          AND STEP LENGTH H=D[2]*SIGN(B-D[3]); A,YA,ZA ARE IGNORED;
    N:    <ARITHMETIC EXPRESSION>;
          THE NUMBER OF EQUATIONS.

PROCEDURES USED: NONE.

REQUIRED CENTRAL MEMORY:
    EIGHT ARAYS OF ORDER N AND ONE OF ORDER 4 * N ARE USED.

METHOD AND PERFORMANCE :
    RK3N INTEGRATES (D/DX)(D/DX)Y=F(X,Y) FROM X TO B,WITH,IF FI="TRUE"
    THEN X=A, Y[J]=YA[J], Z[J]=ZA[J].IF FI="FALSE" THEN X=D[3],
    Y[J]=D[J+3], Z[J]=D[N+3+J], USING A 5-TH ORDER RUNGE-KUTTA METHOD.
    UPON COMPLETION OF A CALL OF RK3N WE HAVE X=D[3]=B, Y[J]=D[J+3]
    THE VALUE OF THE DEPENDENT VARIABLES FOR X=B, Z[J]= D[N+3+J],
    THE VALUE OF THE DERIVATIVES OF Y[J] AT X=B.
    RK3N USES AS ITS MINIMAL ABSOLUTE STEP LENGTH:
    HMIN=MIN (E[2*J-1]*INT+E[2*J]) ,WITH 1<=J<=2*N 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(Z[J])*E[2*J-1]+E[2*J])*
    ABS(H)/INT OR IF THAT TERM IS GREATER THEN (ABS(FXYJ)*E[2*(J+N)-1]
    +E[2*(J+N)])*ABS(H)/INT FOR ANY VALUE OF J, 1<=J<=N (INT=ABS(B-A)).
    SEE REF[1].

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

EXAMPLE OF USE:
    THE SECOND ORDER (VECTOR) DIFFERENTIAL EQUATION

        (D/DX)(D/DX)Y[1] = +Y[2],

        (D/DX)(D/DX)Y[2] = -Y[1], X>=0,

        Y[1] = Y[2] = 1,

        (D/DX)Y[1] = (D/DX)Y[2] = 0, X = 0,

    WHOSE EXACT SOLUTION IS GIVEN BY

     Y[1]=COSH(X/SQRT(2))*COS(X/SQRT(2))+SINH(X/SQRT(2))*SIN(X/SQRT(2))
     Y[2]=COSH(X/SQRT(2))*COS(X/SQRT(2))-SINH(X/SQRT(2))*SIN(X/SQRT(2))

    CAN BE INTEGRATED BY RK3N BECAUSE THE SECOND DERIVATIVE IS NOT
    EXPRESSED IN THE FIRST. THE PROGRAM READS AS FOLLOWS:

    "BEGIN" "INTEGER" K,B; "REAL" X; "BOOLEAN" FI;
        "ARRAY" Y,YA,Z[1:2],E[1:8],D[0:7];
        "INTEGER" "PROCEDURE" EVEN(N); "VALUE" N; "INTEGER" N;
        EVEN:= "IF" N//2 = N/2 "THEN" +1 "ELSE" -1;
        "PROCEDURE" EXACT(X,Y); "VALUE" X; "REAL" X; "ARRAY" Y;
        "BEGIN" "INTEGER" I,N; "REAL" X2,TERM;
            Y[1]:=Y[2]:=0; TERM:=1; X2:= X*X*.5;
            "FOR" N:=1, N+1 "WHILE" ABS(TERM)>"-14 "DO"
            "BEGIN" "FOR" I:=1,2 "DO"
                Y[I]:=Y[I] + TERM*EVEN((I+N-2)//2);
                TERM:= TERM*X2  /N/(N*2-1)
            "END"
        "END";
        "FOR" K:=1,2,3,4,5,6,7,8 "DO" E[K]:="-7; FI:= "TRUE";
        Y[1]:=Y[2]:=1; Z[1]:=Z[2]:=0; B:=0; AA: B:= B+1;
      RK3N(X,0.0,B,Y,Y,Z,Z,"IF"K=1"THEN"Y[2]"ELSE"-Y[1],K,E,D,FI,2);
        EXACT(X,YA); OUTPUT(61,"("//10B
        "("ABS(YEXACT[1]-Y[1])+ABS(YEXACT[2]-Y[2])=")".10D"(""00")"
           ")", ABS(Y[1]-YA[1])+ABS(YA[2]-Y[2]) );
        FI:="FALSE" ; "IF" B<5 "THEN" "GO TO" AA
    "END"
    RESULTS:
    FOR X=1,2,3,4,5 THE FOLLOWING ERRORS ARE NOTICED (E[K]="-7,
    K=1,...,8):
          ABS(YEXACT[1]-Y[1])+ABS(YEXACT[2]-Y[2])=.0000000005"00
          ABS(YEXACT[1]-Y[1])+ABS(YEXACT[2]-Y[2])=.0000000018"00
          ABS(YEXACT[1]-Y[1])+ABS(YEXACT[2]-Y[2])=.0000000046"00
          ABS(YEXACT[1]-Y[1])+ABS(YEXACT[2]-Y[2])=.0000000126"00
          ABS(YEXACT[1]-Y[1])+ABS(YEXACT[2]-Y[2])=.0000000293"00

SOURCE TEXT(S):
"CODE" 33015;