NUMAL Section 5.2.1.1.3
BEGIN SECTION : 5.2.1.1.3 (November, 1976)
AUTHORS: P.A. BEENTJES, H.G.J. ROZENHART.
INSTITUTE: MATHEMATICAL CENTRE.
RECEIVED: 760201
BRIEF DESCRIPTION:
ARKMAT SOLVES AN INITIAL VALUE PROBLEM, GIVEN AS A SYSTEM OF FIRST
ORDER (NON-LINEAR) DIFFERENTIAL EQUATIONS BY MEANS OF A STABILIZED
RUNGE KUTTA METHOD;
IN PARTICULAR THIS PROCEDURE IS SUITABLE FOR THE INTEGRATION OF
SYSTEMS WHERE THE DEPENDENT VARIABLE AND THE RIGHTHAND SIDE ARE
STORED IN A RECTANGULAR ARRAY INSTEAD OF A VECTOR , I.E.
DU / DT = F( T, U), WHERE U AND F ARE (N * M) MATRICES ( SEE METHOD
AND PERFORMANCE).
KEYWORDS:
MATRIX DIFFERENTIAL EQUATIONS,
INITIAL VALUE PROBLEMS,
EXPLICIT ONE-STEP METHODS,
STABILIZED RUNGE KUTTA METHODS.
CALLING SEQUENCE:
THE DECLARATION OF THE PROCEDURE IN THE CALLING PROGRAM READS:
"PROCEDURE" ARKMAT(T, TE, M, N, U, DER, TYPE, ORDER, SPR, OUT);
"VALUE" M, N, TYPE, ORDER; "INTEGER" M, N, TYPE, ORDER;
"REAL" T, TE, SPR; "ARRAY" U; "PROCEDURE" DER, OUT;
"CODE" 33066;
THE MEANING OF THE FORMAL PARAMETERS IS
T: <VARIABLE>;
THE INDEPENDENT VARIABLE T;
ENTRY: THE INITIAL VALUE T0;
EXIT : THE FINAL VALUE TE;
TE: <ARITHMETIC EXPRESSION>;
ENTRY: THE FINAL VALUE OF T;
M: <ARITHMETIC EXPRESSION>;
NUMBER OF COLUMNS OF U;
N: <ARITHMETIC EXPRESSION>;
NUMBER OF ROWS OF U;
U: <ARRAY IDENTIFIER>;
"ARRAY" U[1:N,1:M];
ENTRY: THE INITIAL VALUES OF THE SOLUTION OF THE SYSTEM OF
DIFFERENTIAL EQUATIONS AT T=T0;
EXIT : THE VALUES OF THE SOLUTION AT T=TE;
DER: <PROCEDURE IDENTIFIER>;
THE HEADING OF THIS PROCEDURE READS:
"PROCEDURE" DER(T, V, FTV); "VALUE" T;
"REAL" T; "ARRAY" V, FTV;
THIS PROCEDURE MUST BE GIVEN BY THE USER AND PERFORMS
AN EVALUATION OF THE RIGHTHAND SIDE F( T, V) OF THE
SYSTEM; UPON COMPLETION OF DER,THE RIGHTHAND SIDE SHOULD
BE STORED IN FTV[1:N,1:M];
TYPE: <VARIABLE>;
ENTRY: THE TYPE OF THE SYSTEM OF DIFFERENTIAL EQUATIONS TO
BE SOLVED;
THE USER SHOULD SUPPLY ONE OF THE FOLLOWING VALUES;
1: IF NO SPECIFICATION OF THE TYPE CAN BE MADE;
2: IF THE EIGENVALUES OF THE JACOBIAN MATRIX OF THE
RIGHTHAND SIDE ARE NEGATIVE REAL;
3: IF THE EIGENVALUES OF THE JACOBIAN MATRIX OF THE
RIGHTHAND SIDE ARE PURELY IMAGINARY;
ORDER: <VARIABLE>;
THE ORDER OF THE RUNGE KUTTA METHOD USED;
ENTRY: FOR TYPE=2 THE USER MAY CHOOSE ORDER=1 OR ORDER=2;
ORDER SHOULD BE 2 FOR THE OTHER TYPES;
SPR: <ARITHMETIC EXPRESSION>;
ENTRY: THE SPECTRAL RADIUS OF THE JACOBIAN MATRIX OF THE
RIGHTHAND SIDE, WHEN THE SYSTEM IS WRITTEN IN ONE
DIMENSIONAL FORM (I.E. VECTORFORM);
THE INTEGRATION STEP WILL EQUAL CONSTANT/SPR (SEE DATA AND
RESULTS);
IF NECESSARY SPR CAN BE UPDATED (AFTER EACH STEP) BY MEANS
OF THE PROCEDURE OUT;
OUT: <PROCEDURE IDENTIFIER>
THE HEADING OF THIS PROCEDURE READS:
"PROCEDURE" OUT;
AFTER EACH INTEGRATION STEP PERFORMED, INFORMATION CAN BE
OBTAINED OR UPDATED BY THIS PROCEDURE, E.G. THE VALUES OF
T, U[1:N,1:M] AND SPR.
DATA AND RESULTS:
IF THE USER WANTS TO PERFORM THE INTEGRATION WITH A PRESCRIBED STEP
H, HE HAS TO GIVE SPR THE VALUE CONSTANT/H, WHERE CONSTANT HAS THE
FOLLOWING VALUES:
CONSTANT= 4.3 IF TYPE=1 AND ORDER=2;
CONSTANT= 156 IF TYPE=2 AND ORDER=1;
CONSTANT= 64 IF TYPE=2 AND ORDER=2;
CONSTANT= 8 IF TYPE=3 AND ORDER=2;
PROCEDURES USED:
ELMCOL = CP34023,
DUPMAT = CP31035.
REQUIRED CENTRAL MEMORY:
TWO AUXILIARY ARRAYS OF ORDER N*M ARE DECLARED.
METHOD AND PERFORMANCE:
ARKMAT IS AN IMPLEMENTATION OF LOW ORDER STABILIZED RUNGE KUTTA
METHODS (SEE REFERENCE[1]);
THE INTEGRATION STEPSIZE USED WILL DEPEND ON:
1. THE TYPE OF SYSTEM TO BE SOLVED (I.E. HYPERBOLIC OR PARABOLIC);
2. THE SPECTRAL RADIUS OF THE JACOBIAN MATRIX OF THE SYSTEM;
3. THE INDICATED ORDER OF THE PARTICULAR RUNGE KUTTA METHOD;
THE PROCEDURE ARKMAT IS ESPECIALLY INTENDED FOR SYSTEMS OF
DIFFERENTIAL EQUATIONS ARISING FROM INITIAL BOUNDARY VALUE PROBLEMS
IN TWO DIMENSIONS, E.G. WHEN THE METHOD OF LINES IS APPLIED TO THIS
KIND OF PROBLEMS,THE RIGHTHAND SIDE OF THE RESULTING SYSTEM IS MUCH
EASIER TO DESCRIBE IN MATRIX THAN IN VECTOR FORM; BECAUSE OF THIS
FACT THE ARRAY OF DEPENDENT VARIABLES U IS A MATRIX, RATHER THAN A
VECTOR.
REFERENCE:
[1]. P.J. VAN DER HOUWEN.
STABILIZED RUNGE KUTTA METHOD WITH LIMITED
STORAGE REQUIREMENTS.
MATH. CENTR. REPORT TW 124/71.
EXAMPLE OF USE:
GIVEN THE FOLLOWING SYSTEM OF EQUATIONS:
DU / DT = V( T, X, Y),
(1)
DV / DT = D( DU / DX) / DX + D( DU / DY) / DY,
( ORIGINATING FROM THE INITIAL BOUNDARY VALUE PROBLEM
D( DU / DT) / DT = D( DU / DX) / DX + D( DU / DY) / DY,
ON THE DOMAIN 0 <= X <= PI , 0 <= Y <= 1 ),
WITH THE FOLLOWING BOUNDARY CONDITIONS:
U( T, 0, Y) = U( T, PI, Y) = U( T, X, 1) = 0,
U( T, X, 0) = SIN( X ) * COS( SQRT( 1 + PI * PI / 4) * T),
AND THE INITIAL VALUES:
U( 0, X, Y) = SIN( X ) * COS( PI * Y / 2),
V( 0, X, Y) = 0;
BY APPLYING THE METHOD OF LINES TO PROBLEM (1), USING A TEN BY TEN
GRID ON THE INDICATED DOMAIN, THE SYSTEM IS TRANSFORMED TO A MATRIX
-DIFFERENTIAL EQUATION; THE SOLUTION OF THE LATTER PROBLEM AT T=1
IS COMPUTED BY THE FOLLOWING PROGRAM, USING A CONSTANT STEPSIZE .1;
"BEGIN" "REAL" HPI,H1,H2,H1K,H2K,T,TE;
"INTEGER" I,J,N,M,TYP,ORDE,TEL;"ARRAY" U[1:20,1:10];
"PROCEDURE" DERIV(T,U,DU); "VALUE" T; "REAL" T;"ARRAY" U,DU;
"BEGIN" "FOR" I:=2 "STEP" 1 "UNTIL" N-1 "DO"
"FOR" J:=2 "STEP" 1 "UNTIL" M-1 "DO"
"BEGIN" DU[I,J]:=U[I+N,J];
DU[I+N,J]:=(U[I,J+1]-2*U[I,J]+U[I,J-1])/H1K+
(U[I+1,J]-2*U[I,J]+U[I-1,J])/H2K
"END";
"FOR" J:=1,M "DO"
"BEGIN" INIMAT(N+1,N+N,J,J,DU,0);
"FOR" I:=1 "STEP" 1 "UNTIL" N "DO" DU[I,J]:=U[N+1,J]
"END";
"FOR" I:=1,N "DO"
"FOR" J:=2 "STEP" 1 "UNTIL" M-1 "DO"
"BEGIN" DU[I,J]:=U[I+N,J];
"IF" I=1 "THEN" DU[N+1,J]:=(U[1,J+1]-2*U[1,J]+U[1,J-1])/H1K+
(2*U[2,J]-2*U[1,J])/H2K
"ELSE" DU[2*N,J]:=0
"END"
"END" DERIV;
"PROCEDURE" OUT;
"BEGIN" TEL:=TEL+1;
"IF" T=TE "THEN"
"BEGIN" OUTPUT(61,"("//,3B,"("X")",7B,"("Y")",4B,
"("U(1,X,Y)")",7B,"("U(1,X,Y)")",/,16B,"("COMPUTED")",7B,
"("EXACT")",//")");
"FOR" I:= 1 "STEP" 1 "UNTIL" 10 "DO"
OUTPUT(61,"("2(-D.3D2B),2(-D.6D6B),/")",
(I-1)*H1,(I-1)*H2,U[I,I],SIN(H1*(I-1))*COS(HPI*H2*(I-1))*
COS(T*SQRT(1+HPI*HPI)));
OUTPUT(61,"("/,"("NUMBER OF INTEGRATION STEPS: ")"
,ZZZD")",TEL);
OUTPUT(61,"("//,"(" TYPE IS:")",ZD,"(" ORDER IS:")",
ZD")",TYP,ORDE);
"END";
"END" OUT;
"PROCEDURE" START;
"BEGIN" "FOR" J:=1 "STEP" 1 "UNTIL" M "DO" U[N,J]:=SIN(H1*(J-1));
"FOR" I:=1 "STEP" 1 "UNTIL" N "DO"
"BEGIN" "REAL" COS1; COS1:=COS(H2*HPI*(I-1));
"FOR" J:=1 "STEP" 1 "UNTIL" M "DO" U[I,J]:=U[N,J]*COS1
"END";
INIMAT(N+1,N+N,1,M,U,0)
"END" START;
HPI:=2*ARCTAN(1);H2:=1/9;H1:=(2*HPI)/9;N:=M:=10;
H1K:=H1*H1;H2K:=H2*H2;TEL:=0;
T:=0; TE:=1 ; START; TYP:=3; ORDE:=2;
ARKMAT(T,TE,M,N+N,U,DERIV,TYP,ORDE,80.0,OUT)
"END"
THIS PROGRAM DELIVERS:
X Y U(1,X,Y) U(1,X,Y)
COMPUTED EXACT
0.000 0.000 0.000000 0.000000
0.349 0.111 -0.095201 -0.096735
0.698 0.222 -0.170723 -0.173474
1.047 0.333 -0.211983 -0.215398
1.396 0.444 -0.213228 -0.216663
1.745 0.556 -0.178920 -0.181802
2.094 0.667 -0.122388 -0.124360
2.443 0.778 -0.062138 -0.063139
2.793 0.889 -0.016787 -0.017057
3.142 1.000 0.000000 -0.000000
NUMBER OF INTEGRATION STEPS: 10
TYPE IS: 3 ORDER IS: 2
SOURCE TEXT(S):
"CODE" 33066;