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;