NUMAL Section 5.2.1.1.2.1.B

BEGIN SECTION : 5.2.1.1.2.1.B (February, 1979)

PROCEDURE : RK2N.

AUTHOR:J.A.ZONNEVELD.

CONTRIBUTORS: M.BAKKER AND I.BRINK.

INSTITUTE : MATHEMATICAL CENTRE.

RECEIVED: 730715.

BRIEF DESCRIPTION:

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

KEYWORDS :

    INITIAL VALUE PROBLEM,
    SECOND ORDER DIFFERENTIAL EQUATION.

CALLING SEQUENCE:

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

    THE MEANING OF THE FORMAL PARAMETERS IS:
    X:    <VARIABLE>;
          THE INDEPENDENT VARIABLE.
          UPON COMPLETION OF A CALL OF RK2N,
          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;
    Y:    <ARRAY IDENTIFIER>;
          "ARRAY" Y[1:N];
          THE VECTOR OF DEPENDENT VARIABLES;
          EXIT: THE VALUE OF Y[J] (B), (J = 1, .. ,N);
    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 FIRST DERIVATIVES OF THE DEPENDENT VARIABLES;
          EXIT : THE VALUE OF (D/DX)Y[J](B) (J = 1, .. ,N);
    ZA:   <ARRAY IDENTIFIER>;
          "ARRAY" ZA[1:N];
          ENTRY : THE STARTING VALUES OF Z[J],I.E. THE VALUES AT X=A;
    FXYZJ:<ARITHMETIC EXPRESSION>;
          AN EXPRESSION DEPENDING ON X,J,Y[J],Z[J]  (J=1,...,N),
          GIVING THE VALUE OF (D/DX)(D/DX)Y[J];
    J:    <VARIABLE>;
          A VARIABLE OF TYPE INTEGER,USED IN THE ACTUAL PARAMETER
          CORRESPONDING TO FXYZJ,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 INITIAL CONDITIONS:X=D[3],Y[J]=D[J+3],Z[J]=
          D[N+3+J] AND STEP LENGTH H=D[2]*SIGN(B-D[3]), AND
          A, YA, ZA ARE IGNORED;
    N:    <ARITHMETIC EXPRESSION>;
          THE NUMBER OF EQUATIONS.

PROCEDURES USED: NONE.

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

METHOD AND PERFORMANCE :
    RK2N INTEGRATES (D/DX)(D/DX)Y = F(X,Y,Z) FROM X TO B,WITH, EITHER
    (IF FI = "TRUE") X=A, Y[J]=YA[J], Z[J]=ZA[J], OR  (IF FI="FALSE")
    X = D[3], Y[J]=D[J+3], Z[J]=D[N+J+3], J=1,...,N, USING A 5-TH ORDER
    RUNGE-KUTTA METHOD.
    UPON COMPLETION OF A CALL OF RK2N WE HAVE:X=D[3]=B, Y[J]=D[J+3]
    THE VALUE OF THE DEPENDENT VARIABLES FOR X=B, Z[J]=D[N+J+3], THE
    VALUE OF THE DERIVATIVES OF Y[J] AT X=B, J=1,...,N.
    RK2N 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
    COMPUTED DISCRETIZATION ERROR IS GREATER THAN
    ( ABS(Z[J]) * E[2 * J - 1] + E[2 * J] ) * ABS(H) / INT
    OR IF THAT TERM IS GREATER THEN (ABS(FXYZJ)*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].

EXAMPLE OF USE:

    THE SECOND ORDER (VECTOR) DIFFERENTIAL EQUATION

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

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

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

    WITH ANALYTIC SOLUTION

      Y[1] = -EXP(-X)*(EXP(-X)*(EXP(-X)*(EXP(-X)/3+.5)-1)-5/6),
      Y[2] = -EXP(-X)*(EXP(-X)*(EXP(-X)*(EXP(-X)/3-.5)+1)-5/6)

    CAN BE INTEGRATED BY RK2N FROM 0 TO 5 WITH  1,2,3,4 AS REFERENCE
    POINTS. THE PROGRAM READS AS FOLLOWS:

    "BEGIN" "REAL" B, X, EXPX; "INTEGER" K; "BOOLEAN" FI;
        "ARRAY" Y,YA,Z,ZA[0:2],E[1:8],D[0:7];
        "FOR"  K:=1,2,3,4,5,6,7,8 "DO" E[K]:="-7;
        YA[1]:=ZA[2]:=1; YA[2]:=ZA[1]:=0; B:=1; AA: FI:=B=1;
      RK2N(X,0.0,B,Y,YA,Z,ZA,-5*(Y[K]+Z[K])+("IF"K=1"THEN"Y[2]"ELSE"
        Y[1]),K,E,D,FI,2);
        "COMMENT" COMPUTATION  OF THE EXACT VALUES OF Y AND DY/DX;
        EXPX:=EXP(-X);
        YA[1]:=-EXPX*(EXPX*(EXPX*(EXPX/3+.5)-1)-5/6);
        YA[2]:=-EXPX*(EXPX*(EXPX*(EXPX/3-.5)+1)-5/6);
        ZA[1]:=+EXPX*(EXPX*(EXPX*(EXPX/.75+1.5)-2)-5/6);
        ZA[2]:=+EXPX*(EXPX*(EXPX*(EXPX/.75-1.5)+2)-5/6);
        OUTPUT(61,"("/20B"("X=")"D.4D/,
        10B"("Y[1]-YEXACT[1]=")"+.14D ,10B"("Y[2]-YEXACT[2]=")"+.14D4/,
        10B"("Z[1]-ZEXACT[1]=")"+.14D ,10B"("Z[2]-ZEXACT[2]=")"+.14D
        5/")",X,Y[1]-YA[1],Y[2]-YA[2],Z[1]-ZA[1],Z[2]-ZA[2]);
        B:=B+1; "IF" B<5 "THEN" "GO TO" AA
    "END"
    RESULTS:

          X=1.0000
     Y[1]-YEXACT[1]=+.00000000002955     Y[2]-YEXACT[2]=+.0000000000567
     Z[1]-ZEXACT[1]=-.00000000013770     Z[2]-ZEXACT[2]=-.0000000002422

          X=2.0000
     Y[1]-YEXACT[1]=-.00000000085294     Y[2]-YEXACT[2]=+.0000000001486
     Z[1]-ZEXACT[1]=+.00000000378800     Z[2]-ZEXACT[2]=-.0000000006509

          X=3.0000
     Y[1]-YEXACT[1]=-.00000000162707     Y[2]-YEXACT[2]=-.0000000004796
     Z[1]-ZEXACT[1]=+.00000000803265     Z[2]-ZEXACT[2]=+.0000000019380

          X=4.0000
     Y[1]-YEXACT[1]=-.00000000117993    Y[2]-YEXACT[2]=-.0000000008505
     Z[1]-ZEXACT[1]=+.00000000633393     Z[2]-ZEXACT[2]=+.0000000039114

SOURCE TEXT(S):
"CODE" 33013;