NUMAL Section 5.2.1.1.1.1.E
BEGIN SECTION : 5.2.1.1.1.1.E (February, 1979)
PROCEDURE : RK5NA.
AUTHOR: J.A.ZONNEVELD.
CONTRIBUTORS: M.BAKKER AND I.BRINK.
INSTITUTE: MATHEMATICAL CENTRE.
RECEIVED: 730715.
BRIEF DESCRIPTION:
RK5NA SOLVES AN INITIAL VALUE PROBLEM FOR A SYSTEM OF FIRST ORDER
ORDINARY DIFFERENTIAL EQUATIONS DY / DX = F(X,Y), OF WHICH THE
DERIVATIVE COMPONENTS ARE SUPPOSED TO BECOME LARGE, E.G. IN THE
NEIGHBOURHOOD OF SINGULARITIES. RK5NA INTRODUCES THE ARC LENGTH
AS INTEGRATION VARIABLE. THE SYSTEM IS SUPPOSED TO BE NON-STIFF.
KEYWORDS:
INITIAL VALUE PROBLEM.
SYSTEM OF FIRST ORDER ORDINARY DIFFERENTIAL EQUATIONS.
LARGE DERIVATIVE COMPONENTS.
CALLING SEQUENCE:
THE HEADING OF THE PROCEDURE READS:
"PROCEDURE" RK5NA(X, XA, B, FXJ, J, E, D, FI, N, L, POS);
"VALUE" FI, N, L, POS; "INTEGER" J, N, L; "BOOLEAN" FI, POS;
"REAL" B, FXJ; "ARRAY" X, XA, E, D;
"CODE" 33018;
THE MEANING OF THE FORMAL PARAMETERS IS:
X: <ARRAY IDENTIFIER>;
"ARRAY" X[0 : N];
THE DEPENDENT VARIABLES;
X[0], ..., X[N] CAN BE USED AS JENSEN PARAMETERS;
XA: <ARRAY IDENTIFIER>;
"ARRAY" XA[0 : N];
ENTRY: THE INITIAL VALUES OF X[J], J = 0, ..., N;
B: <ARITHMETIC EXPRESSION>;
B DEPENDS ON X[0], ..., X[N];
IF,WITHIN SOME TOLERANCE,B = 0 THEN THE INTEGRATION IS TER-
MINATED; SEE ALSO THE EXPLANATION OF THE PARAMETER E;
FXJ: <ARITHMETIC EXPRESSION>;
THE RIGHT HAND SIDE OF THE DIFFERENTIAL EQUATION;
FXJ DEPENDS ON X[0], ..., X[N], J, GIVING THE VALUE OF
DX[J] / DX[0];
J: <VARIABLE>;
J IS USED AS A JENSEN PARAMETER TO DENOTE,IN THE ACTUAL
PARAMETER CORRESPONDING TO FXJ,THE NUMBER OF THE FUNCTION
REQUIRED;
E: <ARRAY IDENTIFIER>;
"ARRAY" E[0 : 2 * N + 3];
ENTRY:
E[2 * J] AND E[2 * J + 1] ARE THE RELATIVE AND THE ABSOLUTE
TOLERANCE,RESPECTIVELY,ASSOCIATED WITH X[J],J = 0,..., N, WHILE
E[2 * N + 2] AND E[2 * N + 3] ARE THE ONES ASSOCIATED WITH B;
D: <ARRAY IDENTIFIER>;
"ARRAY" D[1 : N + 3];
AFTER COMPLETION OF EACH STEP WE HAVE:
ABS(D[1]) THE ARC LENGTH,
D[2] THE STEP LENGTH,
D[I + 3] THE LATEST VALUE OF X[I],
I = 0, ..., N;
FI: <BOOLEAN EXPRESSION>;
IF FI = "TRUE" THEN THE INTEGRATION IS STARTED WITH INITIAL
CONDITIONS X[I] = XA[I], I = 0, ..., N;
IF FI = "FALSE" THEN THE INTEGRATION IS CONTINUED WITH
X[I] = D[I + 3];
N: <ARITHMETIC EXPRESSION>;
THE NUMBER OF EQUATIONS;
L: <ARITHMETIC EXPRESSION>;
1 <= L <= N; SEE THE EXPLANATION OF POS;
POS: <BOOLEAN EXPRESSION>;
IF FI = "TRUE" THEN THE INTEGRATION STARTS IN SUCH A WAY THAT
X[L] INCREASES IF POS = "TRUE" AND DECREASES IF POS = "FALSE".
IF FI = "FALSE" THEN POS IS IGNORED.
PROCEDURES USED : ZEROIN = CP34150.
REQUIRED CENTRAL MEMORY : CIRCA 9 * N MEMORY PLACES.
METHOD AND PERFORMANCE :
RK5NA IS BASED ON A 5-TH ORDER RUNGE-KUTTA METHOD AND INTEGRATES
THE SYSTEM DX[J] / DX[0] = F(J,X[0],..,X[N]) / F(0,X[0],..,X[N]).
THE ARC LENGTH IS INTRODUCED AS INTEGRATION VARIABLE.
THE INTEGRATION PROCESS IS TERMINATED IF SOME CONDITION ON
X[0], .. ,X[N], TO BE SUPPLIED BY THE USER, IS SATISFIED.
RK5NA USES STEPLENGTH AND ERROR CONTROL. DETAILS ABOUT THE
PROCEDURE AND THE UNDERLYING THEORY ARE GIVEN IN [1] ( RK5NA IS
A SLIGHTLY ADAPTED VERSION OF RK5N).
REFERENCES:
[1]. J.A.ZONNEVELD,
AUTOMATIC NUMERICAL INTEGRATION,
MATH.CENTRE TRACT 8 (1970) .
EXAMPLE OF USE:
THE VAN DER POL EQUATION IN THE PHASE PLANE
DX[1] / DX[0] = (10*(1-X[0]**2)*X[1]-X[0])/X[1]
CAN BE INTEGRATED BY THE PROCEDURE RK5NA; THE STARTING VALUES ARE
X[0] = 2, X[1] = 0. THE INTEGRATION PROCEEDS UNTIL THE NEXT ZERO OF
X[1],THEN IT CONTINUES UNTIL THE NEXT ZERO AND SO ON UNTIL THE
FOURTH ZERO IS ENCOUNTERED.IN THE EXAMPLE THE OUTPUT IS GIVEN FOR
THE TOLERANCES E[0]=E[1]=E[2]=E[3]= "-6, E[4]=E[5]= "-10.
THE PROGRAM READS AS FOLLOWS:
"BEGIN" "COMMENT" VAN DER POL IN THE PHASE PLANE;
"INTEGER" J,K; "BOOLEAN" FIRST;
"ARRAY" E[0:5],XA,X[0:1],D[1:4];
"PROCEDURE" PRINT(X); "ARRAY" X;
"BEGIN" OUTPUT(61,"("/B"("X[0]=")"+D.10D,
10B"("X[1]=")"+D.10D,10B"("S=")"3D.10D")",X[0],
X[1],ABS(D[1]))
"END";
"FOR" K:=0,1,2,3 "DO" E[K]:="-6 ; E[4]:=E[5]:="-10 ; D[1]:=0;
XA[0]:=2; XA[1]:=0; J:=0; PRINT(XA); AA:
FIRST:=J=0;
RK5NA(X,XA,X[1],"IF" K=0 "THEN" X[1] "ELSE"
10*(1-X[0]**2)*X[1]-X[0],K,E,D,FIRST,1,1,"FALSE");;J:=J+1;
PRINT(X); "IF" J<4 "THEN" "GO TO" AA
"END"
RESULTS:
X[0]=+2.0000000000 X[1]=+0.0000000000 S=000.0000000000
X[0]=-2.0142853657 X[1]=-0.0000000012 S=029.3873834087
X[0]=+2.0142853659 X[1]=+0.0000000001 S=058.7884331939
X[0]=-2.0142853659 X[1]=-0.0000000000 S=088.1894829781
X[0]=+2.0142853659 X[1]=+0.0000000000 S=117.5905327623
SOURCE TEXT(S):
"CODE" 33018;