NUMAL Section 5.2.1.2.2.1.2

BEGIN SECTION : 5.2.1.2.2.1.2 (December, 1979)

AUTHORS: T.M.T.COOLEN AND R.PLOEGER.

INSTITUTE: MATHEMATICAL CENTRE, AMSTERDAM.

RECEIVED: 740301.

BRIEF DESCRIPTION:

    THIS SECTION CONTAINS TWO PROCEDURES :

    RICHARDSON SOLVES A SYSTEM OF LINEAR EQUATIONS WITH A COEFFICIENT
    MATRIX HAVING POSITIVE REAL EIGENVALUES BY MEANS OF A NON-
    STATIONARY SECOND ORDER ITERATIVE METHOD: RICHARDSON'S METHOD.
    SINCE RICHARDSON'S METHOD IS PARTICULARLY SUITABLE FOR SOLVING
    A SYSTEM OF LINEAR EQUATIONS THAT IS OBTAINED BY DISCRETIZING A
    TWO-DIMENSIONAL ELLIPTIC BOUNDARY VALUE PROBLEM, THE PROCEDURE
    RICHARDSON IS PROGRAMMED IN SUCH A WAY THAT THE SOLUTION VECTOR
    IS GIVEN AS A TWO-DIMENSIONAL ARRAY U[J,L], LJ<=J<=UJ, LL<=L<=UL.
    THE COEFFICIENT MATRIX IS NOT STORED, BUT EACH ROW CORRESPONDING
    TO A PAIR (J,L) IS GENERATED WHEN NEEDED.
    RICHARSON CAN ALSO BE USED TO DETERMINE THE EIGENVALUE OF THE
    COEFFICIENT MATRIX CORRESPONDING TO THE DOMINANT EIGENFUNCTION.

    ELIMINATION, USED IN CONNECTION WITH THE PROCEDURE RICHARDSON,
    (THIS  SECTION)  SOLVES  A  SYSTEM  OF  LINEAR   EQUATIONS  WITH
    A  COEFFICIENT   MATRIX   HAVING   POSITIVE    REAL
    EIGENVALUES BY MEANS OF A NON-STATIONARY SECOND ORDER ITERATIVE
    METHOD, WHICH IS AN ACCELERATION OF RICHARDSON'S METHOD. IN
    GENERAL, ELIMINATION CANNOT BE USED BY ITSELF IN A SENSIBLE WAY.
    SINCE RICHARDSON'S METHOD AND ITS ACCELERATION ARE PARTICULARLY
    SUITABLE FOR SOLVING A SYSTEM OF LINEAR EQUATIONS THAT IS OBTAINED
    BY DISCRETIZING A TWO-DIMENSIONAL ELLIPTIC BOUNDARY VALUE PROBLEM,
    THE PROCEDURES RICHARDSON AND ELIMINATION ARE PROGRAMMED IN SUCH A
    WAY THAT THE SOLUTION VECTOR IS GIVEN AS A TWO-DIMENSIONAL ARRAY
    U[J,L], LJ<=J<=UJ, LL<=L<=UL. THE COEFFICIENT MATRIX IS NOT STORED,
    BUT EACH ROW CORRESPONDING TO A PAIR(J,L) IS GENERATED WHEN NEEDED.

KEYWORDS:

    DIFFERENTIAL EQUATION,
    TWO-DIMENSIONAL BOUNDARY VALUE PROBLEM,
    SYSTEM OF LINEAR EQUATIONS,
    COEFFICIENT MATRIX HAVING POSITIVE REAL EIGENVALUES,
    NON-STATIONARY SECOND ORDER ITERATIVE METHOD,
    RICHARDSON'S METHOD.
    ACCELERATION OF RICHARDSON'S METHOD.


SUBSECTION : RICHARDSON.

CALLING SEQUENCE:

    THE HEADING OF THE PROCEDURE READS:
    "PROCEDURE" RICHARDSON(U,LJ,UJ,LL,UL,INAP,RESIDUAL,A,B,N,DISCR,K,
    RATECONV,DOMEIGVAL,OUT);
    "VALUE" LJ,UJ,LL,UL,A,B;
    "INTEGER" N,K,LJ,UJ,LL,UL;
    "REAL" A,B,RATECONV,DOMEIGVAL;
    "BOOLEAN" INAP;
    "ARRAY" U,DISCR;
    "PROCEDURE" RESIDUAL, OUT;
    "CODE" 33170;

    THE MEANING OF THE FORMAL PARAMETERS IS:
    U:      <ARRAY IDENTIFIER>;
            "ARRAY" U[LJ:UJ,LL:UL];
            AFTER EACH ITERATION THE APPROXIMATE SOLUTION CALCULATED BY
            THE PROCEDURE RICHARDSON IS STORED INTO U.
            ENTRY: IF INAP IS CHOSEN TO BE "TRUE" THEN AN INITIAL
                   APPROXIMATION OF THE SOLUTION, OTHERWISE ARBITRARY;
            EXIT: THE FINAL APPROXIMATION OF THE SOLUTION;
    LJ,UJ:  <ARITHMETIC EXPRESSION>;
            LOWER AND UPPER BOUND FOR THE FIRST SUBSCRIPT OF U;
    LL,UL:  <ARITHMETIC EXPRESSION>;
            LOWER AND UPPER BOUND FOR THE SECOND SUBSCRIPT OF U;
    INAP:   <BOOLEAN EXPRESSION>;
            IF THE USER WISHES TO INTRODUCE AN INITIAL APPROXIMATION
            INAP="TRUE" SHOULD BE CHOSEN; THE CHOICE INAP="FALSE" HAS
            THE EFFECT THAT ALL COMPONENTS OF U ARE SET EQUAL TO 1
            BEFORE THE FIRST ITERATION IS PERFORMED;
    RESIDUAL: <PROCEDURE IDENTIFIER>;
            THE HEADING OF THIS PROCEDURE READS :
            "PROCEDURE" RESIDUAL(U); "ARRAY" U;
            SUPPOSE THAT THE SYSTEM OF EQUATIONS AT HAND IS AU= F;
            FOR ANY ENTRY U THE PROCEDURE RESIDUAL SHOULD CALCULATE
            THE  RESIDUAL  AU - F  IN  EACH  POINT  J,L,  WHERE
            LJ<=J<=UJ, LL<=L<=UL, AND SUBSTITUTE THESE VALUES IN THE
            ARRAY U;
    A,B:    <ARITHMETIC EXPRESSION>;
            IF ONE WISHES TO FIND THE SOLUTION OF THE BOUNDARY VALUE
            PROBLEM, IN  A  AND  B  THE USER SHOULD GIVE A LOWER AND
            UPPER BOUND FOR THE EIGENVALUES FOR WHICH THE CORRESPONDING
            EIGENFUNCTIONS IN THE EIGENFUNCTION EXPANSION OF THE RESIDU
            AL AU - F, WITH U = THE INITIAL APPROXIMATION, SHOULD BE
            REDUCED; IF THE DOMINANT EIGENVALUE IS TO BE FOUND, ONE
            SHOULD CHOOSE  A  GREATER THAN THIS EIGENVALUE
            (SEE METHOD AND PERFORMANCE);

    N:      <ARITHMETIC EXPRESSION>;
            N GIVES THE TOTAL NUMBER OF ITERATIONS TO BE PERFORMED; THE
            VALUE OF N SHOULD EITHER BE GIVEN, OR MADE DEPENDENT OF
            SOME JENSEN PARAMETER; E.G. K AND RATECONV CAN SERVE
            FOR THIS PURPOSE;
    DISCR:  <ARRAY IDENTIFIER>;
            "ARRAY" DISCR[1:2];
            AFTER EACH ITERATION THE PROCEDURE RICHARDSON DELIVERS
            IN DISCR[1] THE EUCLIDEAN NORM OF THE RESIDUAL, AND
            IN DISCR[2] THE MAXIMUM NORM OF THE RESIDUAL;
    K:      <VARIABLE>
            K COUNTS THE NUMBER OF ITERATIONS RICHARDSON IS PERFORMING;
            IT CAN SERVE AS A JENSEN PARAMETER FOR N AND OUT;
    RATECONV: <VARIABLE>;
            AFTER EACH ITERATION THE AVERAGE RATE OF CONVERGENCE IS
            ASSIGNED TO RATECONV;
    DOMEIGVAL: <VARIABLE>;
            AFTER EACH ITERATION THE VALUE OF THE DOMINANT EIGENVALUE,
            IF PRESENT, IS ASSIGNED TO DOMEIGVAL; IF THERE IS NO
            DOMINANT EIGENVALUE, THE VALUE OF DOMEIGVAL IS MEANINGLESS,
            WHICH MANIFESTS ITSELF BY SHOWING NO CONVERGENCE TO A
            FIXED VALUE;
    OUT:    <PROCEDURE IDENTIFIER>;
            THE HEADING OF THIS PROCEDURE, TO BE WRITTEN BY THE USER,
            READS :
            "PROCEDURE" OUT(K); "VALUE" K; "INTEGER"K;
            BY THIS PROCEDURE ONE HAS ACCESS TO THE FOLLOWING
            QUANTITIES:
            FOR 0<=K<=N THE K-TH ITERAND IN U,THE EUCLIDEAN AND
            MAXIMUM NORM OF THE K-TH RESIDUAL IN DISCR[1] AND DISCR[2],
            RESPECTIVELY;
            FOR 0<K<=N ALSO THE AVERAGE RATE OF CONVERGENCE AND THE
            APPROXIMATION TO THE DOMINANT EIGENVALUE, BOTH WITH RESPECT
            TO THE K-TH ITERAND U, IN RATECONV AND DOMEIGVAL,
            RESPECTIVELY;
            MOREOVER, OUT CAN BE USED TO LET N BE DEPENDENT ON THE
            ACCURACY REACHED IN APPROXIMATING THE DOMINANT EIGENVALUE.

DATA AND RESULTS: SEE REF[1],[2].

PROCEDURES USED: NONE.

REQUIRED CENTRAL MEMORY:

    APPROXIMATELY  3*(UJ - LJ + 1) * (UL - LL + 1) WORDS.

LANGUAGE: ALGOL 60.

METHOD AND PERFORMANCE:

    SUPPOSE THE SYSTEM OF EQUATIONS TO BE SOLVED READS AU = F, WHERE
    A IS A MATRIX HAVING POSITIVE REAL EIGENVALUES. DENOTING THE
    K-TH ITERATE BY  U(K), U(K) BEING THE VECTOR  U(K)[J,L], LJ<=J<=UJ,
    LL<=L<=UL, THE SO-CALLED RESIDUAL WITH RESPECT TO THE K-TH ITERATE
    IS DEFINED BY
            R(K) = AU(K) - F.
    A SECOND ORDER NON-STATIONARY ITERATIVE METHOD IS GIVEN BY
            U(K+1) = BETA K * U(K) + (1 - BETA K) * U(K-1)
                     - OMEGA K * R(K),
    OR, EQUIVALENTLY, IF U IS THE (UNKNOWN) EXACT SOLUTION OF AU = F,
            U(K) - U = PK(A) (U(0) - U),
    WHERE PK DENOTES A POLYNOMIAL OF DEGREE K. RICHARDSON'S METHOD
    CONSISTS OF CHOOSING THIS POLYNOMIAL IN SUCH A WAY THAT AMONGST ALL
    POLYNOMIALS  PK(X) OF DEGREE K WITH PK(0)= 1 IT HAS MINIMAL MAXIMUM
    NORM OVER THE INTERVAL [C,D], WHERE  C > 0 SHOULD BE CHOSEN TO BE A
    LOWER BOUND, AND D AN UPPER BOUND FOR THE EIGENVALUES OF  A.
    APPLICATION OF THIS POLYNOMIAL TO THE INITIAL ERROR  U(0) - U  HAS
    THE EFFECT THAT EACH COMPONENT OF THE INITIAL ERROR IN ITS EIGEN-
    FUNCTION EXPANSION IS REDUCED BY A FACTOR LESS OR EQUAL TO THE NORM
    OF THE POLYNOMIAL.
    THE POLYNOMIALS
            PK(X) = CK((A+B-2*X)/(A-B)) / CK((A+B)/(A-B))
    WHERE CK(Y) DENOTES THE K-TH CHEBYSHEV POLYNOMIAL, HAVE THE
    DESIRED PROPERTIES. THUS, THE VALUES OF THE PARAMETERS BETA K
    AND OMEGA K MAY BE DETERMINED FROM THE RECURRENCE RELATIONS FOR
    CHEBESHEV POLYNOMIALS.
    IN COMPUTATION U(K) - U  IS NOT AVAILABLE, SO ONE USES R(K) AS
    A MEASURE FOR THE ERROR.

    THE  ELEMENTS  OF  THE  MATRIX  A  ARE  NOT  STORED,  BUT
    GENERATED WHEN NEEDED. MORE PRECISELY, THIS MEANS THAT THE
    (UJ-LJ+1) * (UL-LL+1) COMPONENTS OF  AU(K) - F  ARE CALCULATED FOR
    EACH PAIR (J,L)  LJ<J<UJ, LL<L<UL. THE USER SHOULD INTRODUCE THE
    EQUATION TO BE SOLVED IN THIS MANNER BY MEANS OF THE PROCEDURE
    RESIDUAL.
    CLEARLY, THE METHOD IS PARTICULARLY SUITABLE FOR SPARSE MATRICES,
    FOR EXAMPLE MATRICES THAT ARE OBTAINED BY DISCRETIZING ELLIPTIC
    PARTIAL DIFFERENTIAL EQUATIONS.
    THE SHARPER THE BOUNDS  C  AND  D  FOR THE EIGENVALUES 0F  A ARE,
    THE BETTER APPROXIMATE SOLUTION ONE GETS FOR A GIVEN VALUE OF  K,
    SINCE THE ASYMPTOTIC RATE OF CONVERGENCE (K TO INFINITY) IS
    2 * SQRT(C/D).

    NOW LET ALPHA1 BE THE SMALLEST EIGENVALUE OF  A. IF ONE CHOOSES
    C > ALPHA1, THEN, STARTING WITH ANY INITIAL APPROXIMATION, FOR A
    SUFFICIENTLY LARGE NUMBER OF ITERATIONS THE PROCEDURE RICHARDSON
    WILL DELIVER AN APPROXIMATE VALUE FOR THIS EIGENVALUE.

    LET US EXPLAIN THIS FACT FOR THE CASE  ALPHA1 < C < ALPHA2, WHERE
    ALPHA2  IS THE SECOND SMALLEST EIGENVALUE OF  A. THE POLYNOMIAL
    PK(X)  HAS SMALL MAXIMUM VALUE OVER THE INTERVAL [C,D] (WHICH, OF
    COURSE, DEPENDS ON  K), BUT BECOMES LARGE FOR  X < A. SO, IF ONE
    APPLIES PK(A) TO AN EIGENFUNCTION OF A, THIS EIGENFUNCTION WILL
    ONLY BE REDUCED CONSIDERABLY IF IT CORRESPONDS TO AN EIGENVALUE
    > C.  CONSEQUENTLY, THE EIGENFUNCTION CORRESPONDING TO  ALPHA1 WILL
    BECOME DOMINANT IN THE EIGENFUNCTION EXPANSION OF
            PK(A) (U(0) - U)
    FOR SUFFICIENTLY LARGE K.

    SEE REF[1],[2] FOR DETAILS.

REFERENCES:

    [1].T.M.T.COOLEN, P.W.HEMKER, P.J.VAN DER HOUWEN  AND
        E.SLAGT.
        ALGOL 60 PROCEDURES FOR INITIAL AND BOUNDARY VALUE PROBLEMS
        (DUTCH).
        MC-SYLLABUS 20, MATHEMATICAL CENTRE, 1973, AMSTERDAM.

    [2].P.J.VAN DER HOUWEN.
        FINITE DIFFERENCE METHODS FOR SOLVING PARTIAL DIFFERENTIAL
        EQUATIONS.
        MATHEMATICAL CENTRE TRACT NO. 20, 1968.

EXAMPLE OF USE:

    THE APPROXIMATE SOLUTION OF THE BOUNDARY VALUE PROBLEM
        - ((D/DX)**2 + (D/DY)**2) U(X,Y) = -2*(X*X+Y*Y), O<X,Y<PI,
        U(X,0) = 0, U(X,PI) = PI*PI*X*X,  0 < X < PI,
        U(0,Y) = 0, U(PI,Y) = PI*PI*X*X,  0 < Y < PI,
    WHICH HAS THE ANALYTICAL SOLUTION X*X*Y*Y, MAY BE OBTAINED BY THE
    FOLLOWING PROGRAM:

"BEGIN" "COMMENT" DIRICHLET PROBLEM FOR LAPLACE'S EQUATION;

  "PROCEDURE" RESIDUAL(U); "ARRAY" U;
  "BEGIN" "INTEGER" UJMIN1,ULMIN1,LJPLUS1;
  "REAL" U2; "REAL" "ARRAY" U1[LJ:UJ];
    UJMIN1:= UJ - 1; ULMIN1 := UL - 1; LJPLUS1:= LJ + 1;
    "FOR" J:= LJ "STEP" 1 "UNTIL" UJ "DO"
    "BEGIN" U1[J]:= U[J,LL]; U[J,LL]:= 0; "END";

    "FOR" L:= LL + 1 "STEP" 1 "UNTIL" ULMIN1 "DO"
    "BEGIN"  U1[LJ]:= U[LJ,L]; U[LJ,L]:= 0;
      "FOR"  J:= LJPLUS1"STEP" 1 "UNTIL" UJMIN1 "DO"
      "BEGIN" U2:= U[J,L];
        U[J,L]:=(4 * U2 - U1[J-1]  - U1[J] - U[J+1,L] - U[J,L+1])
        - F(J*H,L*H)*H2;
        U1[J]:= U2
      "END";
      U[UJ,L]:= 0;
    "END";
    "FOR" J:= LJ "STEP" 1 "UNTIL" UJ "DO" U[J,UL]:= 0
  "END" RESIDUAL;

  "REAL" "PROCEDURE" F(X,Y); "VALUE" X,Y; "REAL" X,Y;
  F:= -2*(X*X + Y*Y);

  "REAL" "PROCEDURE" ANALSOL(X,Y); "VALUE" X,Y; "REAL" X,Y;
  ANALSOL:= X*X*Y*Y;

  "PROCEDURE" INITAPPR(U,J,L,G); "INTEGER" J,L; "ARRAY" U; "REAL" G;
  "FOR" J:= LJ "STEP" 1 "UNTIL" UJ "DO"
  "FOR" L:= LL "STEP" 1 "UNTIL" UL "DO"
  U[J,L]:= "IF" J=LJ "OR" J=UJ "OR" L=LL "OR" L=UL "THEN"G "ELSE" 1;

  "PROCEDURE"OUT1(K); "VALUE" K; "INTEGER" K;
  "IF" K = N "THEN" OUTPUT(61,"("//"("  K   DISCR[1]       DISCR[2]
  RATECONV")",//,+ZDB,3(+.7D"+ZDB)")",K,DISCR[1],DISCR[2],RATECONV);
  "INTEGER" J,L,LJ,UJ,LL,UL,N,K;
  "REAL" H,PI,D1,D2,H2,DOMEIGVAL,RATECONV,A,B;
  "REAL" "ARRAY" DISCR[1:2];
  OUTPUT(61,"("/"("GIVE LJ,UJ,LL,UL,N,A,B")"/")");
  ININTEGER(70,LJ); ININTEGER(70,UJ);
  ININTEGER(70,LL); ININTEGER(70,UL);
  ININTEGER(70,N); INREAL(70,A); INREAL(70,B);

  "BEGIN" "REAL" "ARRAY" U[LJ:UJ,LL:UL];
    PI:=3.1415 92653 58979; H:= PI/(UJ - LJ); H2:= H * H;
    INITAPPR(U,J,L,ANALSOL(J*H,L*H));
    RICHARDSON(U,LJ,UJ,LL,UL,"TRUE",RESIDUAL,A,B,N,DISCR,K,
    RATECONV ,DOMEIGVAL,OUT1);
  "END"
"END"

    IT DELIVERS WITH
    LJ = 0, UJ = 11, LL = 0, UL = 11, N = 50, A = .163, B = 7.83
    THE FOLLOWING RESULTS:

   K   DISCR[1]      DISCR[2]      RATECONV

 +50 +.1401828" -3 +.4666866" -4 +.2921718" +0  .


SUBSECTION : ELIMINATION.

CALLING SEQUENCE:

    THE HEADING OF THE PROCEDURE READS:
    "PROCEDURE" ELIMINATION(U,LJ,UJ,LL,UL,RESIDUAL,A,B,N,DISCR,K,
    RATECONV,DOMEIGVAL,OUT);
    "VALUE" LJ,UJ,LL,UL,A,B;
    "INTEGER" N,K,LJ,UJ,LL,UL;
    "REAL" A,B,RATECONV,DOMEIGVAL;
    "ARRAY" U,DISCR;
    "PROCEDURE" RESIDUAL, OUT; "CODE" 33171;

    THE MEANING OF THE FORMAL PARAMETERS IS:
    U:      <ARRAY IDENTIFIER>;
            "ARRAY" U[LJ:UJ,LL:UL];
            AFTER EACH ITERATION THE APPROXIMATE SOLUTION CALCULATED BY
            THE PROCEDURE ELIMINATION IS STORED INTO U;
            ENTRY: AN INITIAL APPROXIMATION OF THE SOLUTION, WHICH
                   IS OBTAINED BY USE OF RICHARDSON;
            EXIT: THE FINAL APPROXIMATION OF THE SOLUTION;
    LJ,UJ:  <ARITHMETIC EXPRESSION>;
            LOWER AND UPPER BOUND FOR THE FIRST SUBSCRIPT OF U;
    LL,UL:  <ARITHMETIC EXPRESSION>;
            LOWER AND UPPER BOUND FOR THE SECOND SUBSCRIPT OF U;
    RESIDUAL: <PROCEDURE IDENTIFIER>;
            THE HEADING OF THIS PROCEDURE READS :
            "PROCEDURE" RESIDUAL(U); "ARRAY" U;
            SUPPOSE THAT THE SYSTEM OF EQUATIONS AT HAND IS AU= F;
            FOR ANY ENTRY U THE PROCEDURE RESIDUAL SHOULD CALCULATE
            THE SO-CALLED RESIDUAL AU - F IN EACH POINT J,L, WHERE
            LJ<=J<=UJ, LL<=L<=UL, AND SUBSTITUTE THESE VALUES IN THE
            ARRAY U;
    A,B:    <ARITHMETIC EXPRESSION>;
            A AND B SHOULD HAVE THE SAME VALUES AS IN THE PRECEDING
            CALL OF RICHARDSON (SEE DESCRIPTION OF RICHARDSON);
    N:      <VARIABLE>;
            THE NUMBER OF ITERATIONS THE PROCEDURE ELIMINATION NEEDS
            TO ELIMINATE THE EIGENFUNCTION BELONGING TO THE DOMINANT
            EIGENVALUE, IS ASSIGNED TO N;
    DISCR:  <ARRAY IDENTIFIER>;
            "ARRAY" DISCR[1:2];
            AFTER EACH ITERATION THE PROCEDURE  ELIMINATION DELIVERS
            IN DISCR[1] THE EUCLIDEAN NORM OF THE RESIDUAL, AND
            IN DISCR[2] THE MAXIMUM NORM OF THE RESIDUAL;
    K:      <VARIABLE>
            K COUNTS THE NUMBER OF ITERATIONS ELIMINATION IS PERFORMING
            IT CAN SERVE AS A JENSEN PARAMETER FOR OUT;
    RATECONV: <VARIABLE>;
            AFTER EACH ITERATION THE AVERAGE RATE OF CONVERGENCE IS
            ASSIGNED TO RATECONV;

    DOMEIGVAL: <ARITHMETIC EXPRESSION>;
            BEFORE A CALL OF ELIMINATION THE VALUE OF THE EIGENVALUE
            FOR WHICH THE CORRESPONDING EIGENFUNCTION HAS TO BE
            ELIMINATED, SHOULD BE ASSIGNED TO DOMEIGVAL; IF AFTER
            APPLICATION OF ELIMINATION THERE IS A NEW DOMINANT EIGEN-
            FUNCTION, THEN DOMEIGVAL WILL BE EQUAL TO THE CORRESPOND-
            ING EIGENVALUE; OTHERWISE, THE VALUE OF DOMEIGVAL BECOMES
            MEANINGLESS;
    OUT:    <PROCEDURE IDENTIFIER>;
            THE HEADING OF THIS PROCEDURE, TO BE WRITTEN BY THE USER,
            READS :
            "PROCEDURE" OUT(K); "VALUE" K; "INTEGER"K;
            BY THIS PROCEDURE ONE HAS ACCESS TO THE FOLLOWING
            QUANTITIES:
            FOR 0<=K<=N THE K-TH ITERAND IN U,THE EUCLIDEAN AND
            MAXIMUM NORM OF THE K-TH RESIDUAL IN DISCR[1] AND DISCR[2],
            RESPECTIVELY;
            FOR 0<K<=N ALSO THE AVERAGE RATE OF CONVERGENCE WITH
            RESPECT TO THE K-TH ITERAND U, IN RATECONV;
            FOR K = N, POSSIBLY THE DOMINANT EIGENVALUE OF THE
            COEFFICIENT MATRIX OF THE EQUATION AU= F, IN DOMEIGVAL.

DATA AND RESULTS: SEE REF[1],[2].

PROCEDURES USED:

    RICHARDSON = CP33170,
    TAN = CP35120,
    TANH = CP35113,
    ARCCOS = CP35122,
    ZEROIN = CP34150.

REQUIRED CENTRAL MEMORY:

    APPROXIMATELY  3*(UJ - LJ + 1) * (UL - LL + 1) WORDS.

METHOD AND PERFORMANCE:

    SEE THIS HEADING IN THE DESCRIPTION OF THE PROCEDURE RICHARDSON.
    SOME ADDITIONAL REMARKS WILL BE MADE HERE.
    IN ORDER TO USE ELIMINATION THE INITIAL APPROXIMATION OF THE
    SOLUTION OF
            AU = F
    IS FIRST TREATED BY MEANS OF RICHARDSON'S METHOD, WHERE  C  IS
    CHOSEN GREATER THAN THE SMALLEST EIGENVALUE. AFTER APPLICATION OF
    RICHARDSON, THE EIGENFUNCTION CORRESPONDING TO THIS EIGENVALUE HAS
    BECOME DOMINANT IN THE QUANTITY
            PK(A) (U(0) - U),
    WITH
            PK(X) = CK((C+D-2*X)/(C-D)) / CK((C+D)/(C-D)),
    WHEREAS THE CONTRIBUTION OF THE OTHER EIGENFUNCTIONS TO THE ERROR
    U(K) - U  AND TO R(K) HAS BEEN REDUCED CONSIDERABLY. CONSEQUENTLY
    THE ERROR U(K) - U HAS VERY SMALL COMPONENTS IN THE SUBSPACE
    SPANNED BY ALL EIGENVECTORS BUT THE "FIRST", IN WHICH DIRECTION IT
    HAS A VERY LARGE COMPONENT.
    THE CONTRIBUTION OF THE "FIRST" EIGENFUNCTION TO  R(K)  IS NOW
    "ELIMINATED" BY APPLICATION OF A POLYNOMIAL OPERATOR E(A) SUCH
    THAT E(X) HAS A ZERO IN THE FIRST EIGENVALUE.
    THE POLYNOMIAL IS CHOSEN IN SUCH A WAY THAT A MAXIMAL RATE OF CON-
    VERGENCE WITH RESPECT TO THE INITIAL APPROXIMATION USED IN
    RICHARDSON IS OBTAINED.

    FOR DETAILS SEE REF[1],[2].

REFERENCES:

    [1].T.M.T.COOLEN, P.W.HEMKER, P.J.VAN DER HOUWEN  AND
        E.SLAGT.
        ALGOL 60 PROCEDURES FOR INITIAL AND BOUNDARY VALUE PROBLEMS
        (DUTCH).
        MC-SYLLABUS 20, MATHEMATICAL CENTRE, 1976, AMSTERDAM.

    [2].P.J.VAN DER HOUWEN.
        FINITE DIFFERENCE METHODS FOR SOLVING PARTIAL DIFFERENTIAL
        EQUATIONS.
        MATHEMATICAL CENTRE TRACT NO. 20, 1968.

EXAMPLE OF USE:

    THE APPROXIMATE SOLUTION OF THE BOUNDARY VALUE PROBLEM
        - ((D/DX)**2 + (D/DY)**2) U(X,Y) = -2*(X*X+Y*Y), O<X,Y<PI,
        U(X,0) = 0, U(X,PI) = PI*PI*X*X,  0 < X < PI,
        U(0,Y) = 0, U(PI,Y) = PI*PI*X*X,  0 < Y < PI,
    WHICH HAS THE ANALYTICAL SOLUTION X*X*Y*Y, MAY BE OBTAINED BY THE
    FOLLOWING PROGRAM:

"BEGIN" "COMMENT" DIRICHLET PROBLEM FOR LAPLACE'S EQUATION;

  "PROCEDURE" RESIDUAL(U); "ARRAY" U;
  "BEGIN" "INTEGER" UJMIN1,ULMIN1,LJPLUS1;
  "REAL" U2; "REAL" "ARRAY" U1[LJ:UJ];
    UJMIN1:= UJ - 1; ULMIN1 := UL - 1; LJPLUS1:= LJ + 1;
    "FOR" J:= LJ "STEP" 1 "UNTIL" UJ "DO"
    "BEGIN" U1[J]:= U[J,LL]; U[J,LL]:= 0; "END";
    "FOR" L:= LL + 1 "STEP" 1 "UNTIL" ULMIN1 "DO"
    "BEGIN"  U1[LJ]:= U[LJ,L]; U[LJ,L]:= 0;
      "FOR"  J:= LJPLUS1"STEP" 1 "UNTIL" UJMIN1 "DO"
      "BEGIN" U2:= U[J,L];
        U[J,L]:=(4 * U2 - U1[J-1]  - U1[J] - U[J+1,L] - U[J,L+1])
        - F(J*H,L*H)*H2;
        U1[J]:= U2
      "END";
      U[UJ,L]:= 0;
    "END";
    "FOR" J:= LJ "STEP" 1 "UNTIL" UJ "DO" U[J,UL]:= 0
  "END" RESIDUAL;

  "REAL" "PROCEDURE" F(X,Y); "VALUE" X,Y; "REAL" X,Y;
  F:= -2*(X*X + Y*Y);

  "REAL" "PROCEDURE" ANALSOL(X,Y); "VALUE" X,Y; "REAL" X,Y;
  ANALSOL:= X*X*Y*Y;

  "PROCEDURE" INITAPPR(U,J,L,G); "INTEGER" J,L; "ARRAY" U; "REAL" G;
  "FOR" J:= LJ "STEP" 1 "UNTIL" UJ "DO"
  "FOR" L:= LL "STEP" 1 "UNTIL" UL "DO"
  U[J,L]:= "IF" J=LJ "OR" J=UJ "OR" L=LL "OR" L=UL "THEN"G "ELSE" 1;

  "PROCEDURE"OUT3(K); "VALUE" K; "INTEGER" K;
  "IF" K=P "THEN" OUTPUT(61,"("//,+ZDB,3(+.7D"+ZDB)")",K,DISCR[1],
  DISCR[2],RATECONV);

  "PROCEDURE"OUT1(K); "VALUE" K; "INTEGER" K;
  "IF" K=N "THEN" OUTPUT(61,"("//"("  K   DISCR[1]      DISCR[2]")",
  "("      RATECONV")",//,+ZDB,3(+.7D"+ZDB)")",
     K,DISCR[1],DISCR[2],RATECONV);

  "PROCEDURE" OUT2(K); "VALUE" K; "INTEGER" K;
  "BEGIN"
    "IF" K = 0 "THEN" D1:= D2:= 1 "ELSE"
    "BEGIN" D2:= D1; D1:= DOMEIGVAL;
      N:= "IF" ABS((D1 - D2)/D2) < 10.0**(-Q) "THEN" K "ELSE" NN;
      OUT1(K)
    "END"
  "END" OUT2;

  "INTEGER" J,L,LJ,UJ,LL,UL,NN,N,P,K,Q;
  "REAL" H,PI,D1,D2,H2,RATECONVR,RATECONVE,DOMEIGVAL,RATECONV,A,B,VAR;
  "REAL" "ARRAY" DISCR[1:2];
  OUTPUT(61,"("/"("GIVE LJ,UJ,LL,UL,N,Q,A,B")"/")");
  ININTEGER(70,LJ); ININTEGER(70,UJ);
  ININTEGER(70,LL); ININTEGER(70,UL);
  ININTEGER(70, N); ININTEGER(70, Q);
  INREAL(70, A); INREAL(70, B);

  "BEGIN" "REAL" "ARRAY" U[LJ:UJ,LL:UL];
    PI:=3.1415 92653 58979; H:= PI/(UJ - LJ); H2:= H * H;
    INITAPPR(U,J,L,ANALSOL(J*H,L*H));
    NN:= N;
    RICHARDSON(U,LJ,UJ,LL,UL,"TRUE",RESIDUAL,A,B,N,DISCR,K,
    RATECONV ,DOMEIGVAL,OUT2); RATECONVR:= RATECONV;
    OUTPUT(61,"("//+.7D"+ZD4B"("DOMINANT EIGENVALUE")"")",DOMEIGVAL);
    ELIMINATION(U,LJ,UJ,LL,UL,RESIDUAL,A ,B,P,DISCR,K,
    RATECONV ,DOMEIGVAL,OUT3); RATECONVE:= RATECONV;
    NN:= N + P; OUTPUT(61,"("//+Z2D13B"("TOTAL NUMBER OF ITERATIONS")"
    ")",NN);
    OUTPUT(61,"("/+.7D"+ZD4B"("RATE OF CONVERGENCE WITH RESPECT TO")",
    /17B"("THE ZEROTH ITERAND OF RICHARDSON")"")",
    (N * RATECONVR + P * RATECONVE)/NN);
  "END"
"END"

    IT DELIVERS WITH
    LJ = 0, UJ = 11, LL = 0, UL = 11, N = 50, Q = 4, A = .326, B = 7.83
    THE FOLLOWING RESULTS:

   K   DISCR[1]      DISCR[2]      RATECONV

 +45 +.4998463" -1 +.8903863" -2 +.2009943" +0

 +.1620445" +0    DOMINANT EIGENVALUE

  +7 +.3563865" -5 +.6714375" -6 +.1360086" +1

  +52             TOTAL NUMBER OF ITERATIONS
 +.3570259" +0    RATE OF CONVERGENCE WITH RESPECT TO
                  THE ZEROTH ITERAND OF RICHARDSON

SOURCE TEXT(S):
"CODE" 33170;
"CODE" 33171;