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;