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;