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;