NUMAL Section 5.2.1.2.1.2.1.1
BEGIN SECTION : 5.2.1.2.1.2.1.1 (December, 1979)
AUTHOR: M. BAKKER.
INSTITUTE: MATHEMATICAL CENTRE, AMSTERDAM.
RECEIVED: 751231/ REVISED 791231.
BRIEF DESCRIPTION:
THIS SECTION CONTAINS THREE PROCEDURES FOR THE SOLUTION
OF SECOND ORDER SELF-ADJOINT LINEAR TWO POINT
BOUNDARY VALUE PROBLEMS;
(1) FEM LAG SYM;
THIS PROCEDURE SOLVES THE DIFFERENTIAL EQUATION
- (P(X)*Y')' + R(X)*Y = F(X), A < X < B,
WITH BOUNDARY CONDITIONS
E[1]*Y(A) + E[2]*Y'(A) = E[3],
E[4]*Y(B) + E[5]*Y'(B) = E[6].
(2) FEM LAG;
THIS PROCEDURE SOLVES THE DIFFERENTIAL EQUATION
- Y'' + R(X)*Y = F(X), A < X < B,
WITH BOUNDARY CONDITIONS
E[1]*Y(A) + E[2]*Y'(A) = E[3],
E[4]*Y(B) + E[5]*Y'(B) = E[6].
(3) FEM LAG SPHER:
THIS PROCEDURE SOLVES THE DIFFERENTIAL EQUATION
WITH SPHERICAL COORDINATES
- (X**NC*Y')'/X**NC + R(X)*Y = F(X), A < X < B,
WITH BOUNDARY CONDITIONS
E[1]*Y(A) + E[2]*Y'(A) = E[3],
E[4]*Y(B) + E[5]*Y'(B) = E[6].
KEY WORDS AND PHRASES:
SECOND ORDER DIFFERENTIAL EQUATIONS,
TWO POINT BOUNDARY VALUE PROBLEMS,
SELF-ADJOINT BOUNDARY VALUE PROBLEMS,
RITZ-GALERKIN METHOD,
SPHERICAL COORDINATES,
GLOBAL METHODS.
LANGUAGE: ALGOL 60.
REFERENCES:
[1] STRANG, G. AND G.J. FIX,
AN ANALYSIS OF THE FINITE ELEMENT METHOD,
PRENTICE-HALL, ENGLEWOOD CLIFFS, NEW JERSEY, 1973.
[2] BAKKER, M., EDITOR,
COLLOQUIUM ON DISCRETIZATION METHODS, CHAPTER 3 (DUTCH),
MATHEMATISCH CENTRUM, MC-SYLLABUS, 1976.
[3] HEMKER, P.W.,
GALERKIN'S METHOD AND LOBATTO POINTS,
MATHEMATISCH CENTRUM, REPORT 24/75 (1975).
[4] BABUSKA, I.,
NUMERICAL STABILITY IN PROBLEMS OF LINEAR ALGEBRA,
S.I.A.M. J. NUM. ANAL., VOL.9, P. 53-77 (1972).
SUBSECTION: FEM LAG SYM.
CALLING SEQUENCE:
THE HEADING OF THE PROCEDURE READS:
"PROCEDURE" FEM LAG SYM(X, Y, N, P, R, F, ORDER, E);
"VALUE" N, ORDER; "INTEGER" N, ORDER;
"ARRAY" X, Y, E;
"REAL" "PROCEDURE" P, R, F;
"CODE" 33300;
THE MEANING OF THE FORMAL PARAMETERS IS:
N: <ARITHMETIC EXPRESSION>;
THE UPPER BOUND OF THE ARRAYS X AND Y; N > 1;
X: <ARRAY IDENTIFIER>;
"ARRAY" X[0:N];
ENTRY: A = X[0] < X[1] < ... < X[N] = B
IS A PARTITION OF THE INTERVAL [A,B];
Y: <ARRAY IDENTIFIER>;
"ARRAY" Y[0:N];
EXIT: Y[I] (I = 0, 1, ... , N) IS THE APPROXIMATE
SOLUTION AT X[I] OF THE DIFFERENTIAL EQUATION
(1) - (P(X)*Y')' + R(X)*Y = F(X), A < X < B,
WITH BOUNDARY CONDITIONS
E[1]*Y(A) + E[2]*Y'(A) = E[3],
(2)
E[4]*Y(B) + E[5]*Y'(B) = E[6];
P: <PROCEDURE IDENTIFIER>;
THE HEADING OF P READS:
"REAL" "PROCEDURE" P(X); "VALUE" X; "REAL" X;
P(X) IS THE COEFFICIENT OF Y' IN (1);
R: <PROCEDURE IDENTIFIER>;
THE HEADING OF R READS:
"REAL" "PROCEDURE" R(X); "VALUE" X; "REAL" X;
R(X) IS THE COEFFICIENT OF Y IN (1);
F: <PROCEDURE IDENTIFIER>;
THE HEADING OF F READS:
"REAL" "PROCEDURE" F(X); "VALUE" X; "REAL" X;
F(X) IS THE RIGHT HAND SIDE OF (1);
ORDER: <ARITHMETIC EXPRESSION>;
ENTRY: ORDER DENOTES THE ORDER OF ACCURACY REQUIRED FOR THE
APPROXIMATE SOLUTION OF (1)-(2); LET H = MAX(X[I] - X[I-1]);
THEN ABS(Y[I] - Y(X[I])) <= C*H**ORDER, I = 0, ... , N;
ORDER CAN BE CHOSEN EQUAL TO 2, 4 OR 6 ONLY;
E: <ARRAY IDENTIFIER>;
"ARRAY" E[1:6];
E[1], ... , E[6] DESCRIBE THE BOUNDARY CONDITIONS (2);
E[1] AND E[4] ARE NOT ALLOWED TO VANISH BOTH.
PROCEDURES USED: NONE.
REQUIRED CENTRAL MEMORY:
FOUR AUXILIARY ARRAYS OF N REALS ARE USED.
RUNNING TIME:
LET K = ORDER/2; THEN
(A) K*N + 1 EVALUATIONS OF P(X), R(X) AND F(X) ARE NEEDED;
(B) ABOUT 17*2**(K-1)*N MULTIPLICATIONS/DIVISIONS ARE NEEDED.
DATA AND RESULTS:
THE PROCEDURE FEM LAG SYM HAS SOME RESTRICTIONS IN ITS USE:
(I) P(X) SHOULD BE POSITIVE ON THE CLOSED INTERVAL <X[J-1],X[J]>;
(II) P(X), R(X) AND F(X) ARE REQUIRED TO BE SUFFICIENTLY SMOOTH
ON <X[0],X[N]> EXCEPT AT THE GRID POINTS WHERE P(X) SHOULD BE
AT LEAST CONTINUOUS;
IN THAT CASE THE ORDER OF ACCURACY (2, 4, OR 6) IS PRESERVED;
(III) R(X) SHOULD BE NONNEGATIVE ON <X[0],X[N]>;
IF, HOWEVER, THE PROBLEM HAS PURE DIRICHLET BOUNDARY CONDITIONS
(I.E. E[2] = E[5] = 0) THIS CONDITION CAN BE WEAKENED TO THE
REQUIREMENT THAT
R(X) > - P0*(PI/(X[N] - X[0]))**2,
WHERE P0 IS THE MINIMUM OF P(X) ON <X[0],X[N]> AND PI HAS
THE VALUE 3.14159...; HOWEVER, ONE SHOULD NOTE THAT THE
PROBLEM MAY BE ILL-CONDITIONED WHEN R(X) IS QUITE NEAR THAT
LOWER BOUND; FOR OTHER NEGATIVE VALUES OF R(X) THE EXISTENCE
OF A SOLUTION REMAINS AN OPEN QUESTION;
(IV) THE USER SHOULD NOT EXPECT GREATER ACCURACY THAN 12 DECIMALS
DUE TO THE LOSS OF DIGITS DURING THE EVALUATION OF THE MATRIX
AND THE VECTOR OF THE LINEAR SYSTEM TO BE SOLVED AND DURING ITS
REDUCTION TO A TRIDIAGONAL SYSTEM; WHEN THE SOLUTION OF THE
PROBLEM IS NOT TOO WILD, THIS 12-DIGIT ACCURACY CAN ALREADY BE
OBTAINED WITH A MODERATE MESH SIZE (E.G. < 0.1), PROVIDED THAT
A SIXTH ORDER METHOD IS USED.
METHOD AND PERFORMANCE:
PROBLEM (1)-(2) IS SOLVED BY MEANS OF GALERKIN'S METHOD WITH
CONTINUOUS PIECEWISE POLYNOMIALS (SEE [1], [2]);
THE SOLUTION IS APPROXIMATED BY A FUNCTION WHICH IS CONTINUOUS ON
THE CLOSED INTERVAL <X[0],X[N]> AND A POLYNOMIAL OF DEGREE LESS
THAN OR EQUAL TO K (K = ORDER//2) ON EACH SEGMENT <X[J-1],X[J]>
(J = 1, ..., N); THIS PIECEWISE POLYNOMIAL IS ENTIRELY
DETERMINED BY THE VALUES IT HAS
AT THE KNOTS X[J] AND ON (K-1) INTERIOR KNOTS ON EACH SEGMENT
<X[J-1],X[J]>; THESE VALUES ARE OBTAINED BY THE SOLUTION OF AN
(ORDER + 1)-DIAGONAL LINEAR SYSTEM WITH A SPECIALLY STRUCTURED
MATRIX (SEE [2]); THE ENTRIES OF THE MATRIX AND THE VECTOR ARE
INNER PRODUCTS WHICH ARE APPROXIMATED BY PIECEWISE (K+1)-POINT
LOBATTO QUADRATURE (SEE [3]); THE EVALUATION OF THE MATRIX AND
THE VECTOR IS DONE SEGMENT BY SEGMENT: ON EACH SEGMENT
THE CONTRIBUTIONS TO THE ENTRIES OF THE MATRIX AND THE
VECTOR ARE COMPUTED AND EMBEDDED IN THE GLOBAL MATRIX AND
VECTOR; SINCE THE FUNCTION VALUES ON THE INTERIOR
POINTS OF EACH SEGMENT ARE NOT COUPLED WITH THE FUNCTION
VALUES OUTSIDE THAT SEGMENT, THE RESULTING LINEAR SYSTEM
CAN BE REDUCED TO A TRIDIAGONAL SYSTEM BY MEANS OF STATIC
CONDENSATION (SEE [2]); THE FINAL TRIDIAGONAL SYSTEM,
SINCE IT IS OF FINITE DIFFERENCE TYPE, IS SOLVED BY
MEANS OF BABUSKA'S METHOD (SEE [4]).
EXAMPLE OF USE:
WE SOLVE THE BOUNDARY VALUE PROBLEM
-(Y'*EXP(X))'+Y*COS(X)=EXP(X)*(SIN(X)-COS(X))+SIN(2*X)/2,
0 < X < PI = 3.14159265358979,
Y(0) = Y(PI) = 0;
FOR THE BOUNDARY CONDITIONS THIS MEANS THAT
E[1] = E[4] = 1; E[2] = E[3] = E[5] = E[6] = 0;
THE ANALYTIC SOLUTION IS Y(X) = SIN(X); WE APPROXIMATE
THE SOLUTION ON A UNIFORM GRID, I.E. X[I] = I*PI/N,
I = 0, ..., N; WE CHOOSE N=10,20 AND COMPUTE FOR ORDER = 2,4,6
THE MAXIMUM ERROR; THE PROGRAM READS AS FOLLOWS:
"BEGIN" "INTEGER" N; "FOR" N:= 10, 20 "DO"
"BEGIN" "INTEGER" I, ORDER; "REAL" PI; "ARRAY" X, Y[0:N], E[1:6];
"REAL" "PROCEDURE" R(X); "VALUE" X; "REAL" X;
R:= COS(X);
"REAL" "PROCEDURE" P(X); "VALUE" X; "REAL" X;
P:= EXP(X);
"REAL" "PROCEDURE" F(X); "VALUE" X; "REAL" X;
F:= EXP(X)*(SIN(X)-COS(X)) + SIN(2*X)/2;
E[1]:= E[4]:= 1; E[2]:= E[3]:= E[5]:= E[6]:= 0;
PI:= 3.14159265358979;
"FOR" I:= 0 "STEP" 1 "UNTIL" N "DO" X[I]:= PI*I/N;
OUTPUT(61,"("//,6B"("N=")"ZD")",N);
"FOR" ORDER:= 2, 4, 6 "DO"
"BEGIN" "REAL" RHO, D;
FEM LAG SYM(X, Y, N, P, R, F, ORDER, E);
RHO:= 0;
"FOR" I:= 0 "STEP" 1 "UNTIL" N "DO"
"BEGIN" D:= ABS(Y[I] - SIN(X[I]));
"IF" RHO < D "THEN" RHO:= D
"END";
OUTPUT(61,"("/,16B"("ORDER=")"D,4B"("MAX.ERROR= ")",
D.DD"+ZD")",ORDER,RHO)
"END"
"END"
"END"
RESULTS:
N=10
ORDER=2 MAX. ERROR= 1.36" -2
ORDER=4 MAX. ERROR= 7.55" -5
ORDER=6 MAX. ERROR= 3.48" -8
N=20
ORDER=2 MAX. ERROR= 3.41" -3
ORDER=4 MAX. ERROR= 4.79" -6
ORDER=6 MAX. ERROR= 5.51"-10
ONE OBSERVES THAT THE MAXIMUM ERROR DECREASES BY ABOUT
2**(-ORDER) WHEN THE MESH SIZE IS HALVED.
SUBSECTION: FEM LAG.
CALLING SEQUENCE:
THE HEADING OF THE PROCEDURE READS:
"PROCEDURE" FEM LAG(X, Y, N, R, F, ORDER, E);
"VALUE" N, ORDER; "INTEGER" N, ORDER;
"ARRAY" X, Y, E;
"REAL" "PROCEDURE" R, F;
"CODE" 33301;
THE MEANING OF THE FORMAL PARAMETERS IS:
N: <ARITHMETIC EXPRESSION>;
THE UPPER BOUND OF THE ARRAYS X AND Y; N > 1;
X: <ARRAY IDENTIFIER>;
"ARRAY" X[0:N];
ENTRY: A = X[0] < X[1] < ... < X[N] = B IS A
PARTITION OF THE SEGMENT [A,B];
Y: <ARRAY IDENTIFIER>;
"ARRAY" Y[0:N];
EXIT: Y[I] (I = 0, 1, ... , N) IS THE APPROXIMATE
SOLUTION AT X[I] OF THE DIFFERENTIAL EQUATION
(3) - Y''+ R(X)*Y = F(X), A < X < B,
WITH BOUNDARY CONDITIONS
(4) E[1]*Y(A) + E[2]*Y'(A) = E[3],
E[4]*Y(B) + E[5]*Y'(B) = E[6];
R: <PROCEDURE IDENTIFIER>;
THE HEADING OF R READS:
"REAL" "PROCEDURE" R(X); "VALUE" X; "REAL" X;
R(X) IS THE COEFFICIENT OF Y IN (3);
F: <PROCEDURE IDENTIFIER>;
THE HEADING OF F READS:
"REAL" "PROCEDURE" F(X); "VALUE" X; "REAL" X;
F(X) IS THE RIGHT HAND SIDE OF (3);
ORDER: <ARITHMETIC <EXPRESSION>;
ENTRY: ORDER DENOTES THE ORDER OF ACCURACY REQUIRED FOR THE
APPROXIMATE SOLUTION OF (3)-(4); LET H = MAX(X[I] - X[I-1]);
THEN ABS(Y[I] - Y(X[I])) <= C*H**ORDER, I = 0, ... , N;
ORDER CAN CAN BE CHOSEN EQUAL TO 2, 4 OR 6 ONLY;
E: <ARRAY IDENTIFIER>;
"ARRAY" E[1:6];
E[1], ... , E[6] DESCRIBE THE BOUNDARY CONDITIONS (4);
E[1] AND E[4] ARE NOT ALLOWED TO VANISH BOTH.
PROCEDURES USED: NONE.
REQUIRED CENTRAL MEMORY:
FOUR AUXILIARY ARRAYS OF N REALS ARE USED.
RUNNING TIME:
LET K = ORDER/2; THEN
(A) K*N + 1 EVALUATIONS OF R(X) AND F(X) ARE NEEDED;
(B) ABOUT 12*2**(K-1)*N MULTIPLICATIONS/DIVISIONS ARE NEEDED.
DATA AND RESULTS: SEE PREVIOUS SUBSECTION.
METHOD AND PERFORMANCE: SEE PREVIOUS SUBSECTION.
EXAMPLE OF USE:
WE SOLVE THE BOUNDARY VALUE PROBLEM
- Y'' + Y*EXP(X) = SIN(X)*(1+EXP(X),
0 < X < PI = 3.14159265358979,
Y(0) = Y(PI) = 0;
FOR THE BOUNDARY CONDITIONS THIS MEANS THAT
E[1] = E[4] = 1; E[2] = E[3] = E[5] = E[6] = 0;
THE ANALYTIC SOLUTION IS Y(X) = SIN(X); WE APPROXIMATE
THE SOLUTION ON A UNIFORM GRID, I.E. X[I] = I*PI/N,
I = 0, ..., N; WE CHOOSE N=10,20 AND COMPUTE FOR ORDER = 2,4,6
THE MAXIMUM ERROR; THE PROGRAM READS AS FOLLOWS:
"BEGIN" "INTEGER" N; "FOR" N:= 10, 20 "DO"
"BEGIN" "INTEGER" I, ORDER; "REAL" PI; "ARRAY" X, Y[0:N], E[1:6];
"REAL" "PROCEDURE" R(X); "VALUE" X; "REAL" X;
R:= EXP(X);
"REAL" "PROCEDURE" F(X); "VALUE" X; "REAL" X;
F:= SIN(X)*(1 + EXP(X));
E[1]:= E[4]:= 1; E[2]:= E[3]:= E[5]:= E[6]:= 0;
PI:= 3.14159265358979;
"FOR" I:= 0 "STEP" 1 "UNTIL" N "DO" X[I]:= PI*I/N;
OUTPUT(61,"("//,6B"("N=")"ZD")",N);
"FOR" ORDER:= 2, 4, 6 "DO"
"BEGIN" "REAL" RHO, D;
FEM LAG(X, Y, N, R, F, ORDER, E);
RHO:= 0;
"FOR" I:= 0 "STEP" 1 "UNTIL" N "DO"
"BEGIN" D:= ABS(Y[I] - SIN(X[I]));
"IF" RHO < D "THEN" RHO:= D
"END";
OUTPUT(61,"("/,16B"("ORDER=")"D,4B"("MAX.ERROR= ")",
D.DD"+ZD")",ORDER,RHO)
"END"
"END"
"END"
RESULTS:
N=10
ORDER=2 MAX. ERROR= 1.60" -3
ORDER=4 MAX. ERROR= 1.55" -5
ORDER=6 MAX. ERROR= 7.28"-10
N=20
ORDER=2 MAX. ERROR= 4.01" -4
ORDER=4 MAX. ERROR= 9.80" -7
ORDER=6 MAX. ERROR= 9.38"-12
NOTICE THAT THE MAXIMUM ERROR DECREASES BY ABOUT
2**(-ORDER) WHEN THE MESH SIZE IS HALVED.
SUBSECTION: FEM LAG SPHER.
CALLING SEQUENCE:
THE HEADING OF THE PROCEDURE READS:
"PROCEDURE" FEM LAG SPHER(X, Y, N, NC, R, F, ORDER, E);
"VALUE" N, NC, ORDER; "INTEGER" N, NC, ORDER;
"ARRAY" X, Y, E;
"REAL" "PROCEDURE" R, F;
"CODE" 33308;
THE MEANING OF THE FORMAL PARAMETERS IS:
N: <ARITHMETIC EXPRESSION>;
THE UPPER BOUND OF THE ARRAYS X AND Y; N > 1;
NC: <EXPRESSION>;
IF NC = 0, CARTESIAN COORDINATES ARE USED;
IF NC = 1, POLAR COORDINATES ARE USED;
IF NC = 2, SPHERICAL COORDINATES ARE USED;
X: <ARRAY IDENTIFIER>;
"ARRAY" X[0:N];
ENTRY: A = X[0] < X[1] < ... < X[N] = B IS A
PARTITION OF THE INTERVAL [A,B];
Y: <ARRAY IDENTIFIER>;
"ARRAY" Y[0:N];
EXIT: Y[I] (I = 0, 1, ... , N) IS THE APPROXIMATE
SOLUTION AT X[I] OF THE DIFFERENTIAL EQUATION
(1) - (X**NC*Y')'/X**NC + R(X)*Y = F(X), A < X < B,
WITH BOUNDARY CONDITIONS
E[1]*Y(A) + E[2]*Y'(A) = E[3],
(2)
E[4]*Y(B) + E[5]*Y'(B) = E[6];
R: <PROCEDURE IDENTIFIER>;
THE HEADING OF R READS:
"REAL" "PROCEDURE" R(X); "VALUE" X; "REAL" X;
R(X) IS THE COEFFICIENT OF Y IN (1);
F: <PROCEDURE IDENTIFIER>;
THE HEADING OF F READS:
"REAL" "PROCEDURE" F(X); "VALUE" X; "REAL" X;
F(X) IS THE RIGHT HAND SIDE OF (1);
ORDER: <ARITHMETIC EXPRESSION>;
ENTRY: ORDER DENOTES THE ORDER OF ACCURACY REQUIRED FOR THE
APPROXIMATE SOLUTION OF (1)-(2); LET H = MAX(X[I] - X[I-1]);
THEN ABS(Y[I] - Y(X[I])) <= C*H**ORDER, I = 0, ... , N;
ORDER CAN BE CHOSEN EQUAL TO 2 OR 4 ONLY;
E: <ARRAY IDENTIFIER>;
"ARRAY" E[1:6];
E[1], ... , E[6] DESCRIBE THE BOUNDARY CONDITIONS (2);
E[1] AND E[4] ARE NOT ALLOWED TO VANISH BOTH.
PROCEDURES USED: NONE.
REQUIRED CENTRAL MEMORY:
FOUR AUXILIARY ARRAYS OF N REALS ARE USED.
RUNNING TIME:
LET K = ORDER/2; THEN
(A) K*N EVALUATIONS OF R(X) AND F(X) ARE NEEDED;
(B) IF NC > 0 AND ORDER=4, THEN N SQUARE ROOTS ARE EVALUATED;
DATA AND RESULTS:
THE PROCEDURE FEM LAG SPHER HAS SOME RESTRICTIONS IN ITS USE:
R(X) AND F(X) ARE REQUIRED TO BE SUFFICIENTLY SMOOTH
ON <X[0],X[N]> EXCEPT AT THE GRID POINTS; FURTHERMORE R(X)
SHOULD BE NONNEGATIVE.
METHOD AND PERFORMANCE:
PROBLEM (1)-(2) IS SOLVED BY MEANS OF GALERKIN'S METHOD WITH
CONTINUOUS PIECEWISE POLYNOMIALS (SEE [1], [2]);
THE SOLUTION IS APPROXIMATED BY A FUNCTION WHICH IS CONTINUOUS ON
THE CLOSED INTERVAL <X[0],X[N]> AND A POLYNOMIAL OF DEGREE LESS
THAN OR EQUAL TO K (K = ORDER//2) ON EACH SEGMENT <X[J-1],X[J]>
(J = 1, ..., N); THIS PIECEWISE POLYNOMIAL IS ENTIRELY
DETERMINED BY THE VALUES IT HAS
AT THE KNOTS X[J] AND ON (K-1) INTERIOR KNOTS ON EACH SEGMENT
<X[J-1],X[J]>; THESE VALUES ARE OBTAINED BY THE SOLUTION OF AN
(ORDER + 1)-DIAGONAL LINEAR SYSTEM WITH A SPECIALLY STRUCTURED
MATRIX (SEE [2]); THE ENTRIES OF THE MATRIX AND THE VECTOR ARE
INNER PRODUCTS WHICH ARE APPROXIMATED BY SOME PIECEWISE K-POINT
GAUSSIAN QUADRATURE (SEE [4]); THE EVALUATION OF THE MATRIX AND
THE VECTOR IS DONE SEGMENT BY SEGMENT: ON EACH SEGMENT
THE CONTRIBUTIONS TO THE ENTRIES OF THE MATRIX AND THE VECTOR
ARE COMPUTED AND EMBEDDED IN THE GLOBAL MATRIX AND VECTOR;
SINCE THE FUNCTION VALUES ON THE INTERIOR
POINTS OF EACH SEGMENT ARE NOT COUPLED WITH THE FUNCTION
VALUES OUTSIDE THAT SEGMENT, THE RESULTING LINEAR SYSTEM
CAN BE REDUCED TO A TRIDIAGONAL SYSTEM BY MEANS OF STATIC
CONDENSATION (SEE [2]); THE FINAL TRIDIAGONAL SYSTEM,
SINCE IT IS OF FINITE DIFFERENCE TYPE, IS SOLVED BY
MEANS OF BABUSKA'S METHOD (SEE [3]).
EXAMPLE OF USE:
WE SOLVE THE BOUNDARY VALUE PROBLEM
-(Y'*X**NC)'/X**NC + Y = 1 - X**4 + (12 + 4*NC)*X**2,
0 < X < 1; Y'(0) = Y(1) = 0;
FOR THE BOUNDARY CONDITIONS THIS IMPLIES THAT
E[2] = E[4] = 1; E[1] = E[3] = E[5] = E[6] = 0;
THE ANALYTIC SOLUTION IS Y(X) = 1 - X**4; WE APPROXIMATE
THE SOLUTION ON A UNIFORM GRID, I.E. X[I] = I/N, I = 0, ..., N;
I = 0, ..., N; WE CHOOSE N=10,20 AND COMPUTE FOR ORDER = 2,4
THE MAXIMUM ERROR; THE PROGRAM READS AS FOLLOWS:
"BEGIN" "INTEGER" N, NC;
"FOR" N:= 10, 20 "DO" "FOR" NC:= 0, 1, 2 "DO"
"BEGIN" "INTEGER" I, ORDER; "ARRAY" X, Y[0:N], E[1:6];
"REAL" "PROCEDURE" R(X); "VALUE" X; "REAL" X;
R:= 1;
"REAL" "PROCEDURE" F(X); "VALUE" X; "REAL" X;
F:= (12 + 4*NC)*X**2 + 1 - X**4;
E[2]:= E[4]:= 1; E[1]:= E[3]:= E[5]:= E[6]:= 0;
"FOR" I:= 0 "STEP" 1 "UNTIL" N "DO" X[I]:= I/N;
OUTPUT(61,"("//,6B"("N=")"ZZD,6B"("NC=")"ZD")",N,NC);
"FOR" ORDER:= 2, 4 "DO"
"BEGIN" "REAL" RHO, D;
FEM LAG SPHER(X, Y, N, NC, R, F, ORDER, E);
RHO:= 0;
"FOR" I:= 0 "STEP" 1 "UNTIL" N "DO"
"BEGIN" D:= ABS(Y[I] - 1 + X[I]**4);
"IF" RHO < D "THEN" RHO:= D
"END";
OUTPUT(61,"("/,16B"(" ORDER=")"ZD,4B"("MAX.ERROR= ")",
D.DD"+ZD")",ORDER,RHO)
"END"
"END"
"END"
RESULTS:
N= 10 NC= 0
ORDER= 2 MAX.ERROR= 4.37" -3
ORDER= 4 MAX.ERROR= 2.93" -6
N= 10 NC= 1
ORDER= 2 MAX.ERROR= 1.42" -2
ORDER= 4 MAX.ERROR= 5.49" -5
N= 10 NC= 2
ORDER= 2 MAX.ERROR= 2.46" -2
ORDER= 4 MAX.ERROR= 1.27" -4
N= 20 NC= 0
ORDER= 2 MAX.ERROR= 1.09" -3
ORDER= 4 MAX.ERROR= 1.83" -7
N= 20 NC= 1
ORDER= 2 MAX.ERROR= 3.53" -3
ORDER= 4 MAX.ERROR= 3.91" -6
N= 20 NC= 2
ORDER= 2 MAX.ERROR= 6.10" -3
ORDER= 4 MAX.ERROR= 9.26" -6
ONE OBSERVES THAT THE MAXIMUM ERROR DECREASES BY ABOUT
2**(-ORDER) WHEN THE MESH SIZE IS HALVED.
SOURCE TEXT(S):
"CODE" 33300;
"CODE" 33301;
"CODE" 33308;