NUMAL Section 5.2.1.2.1.2.2.1
BEGIN SECTION : 5.2.1.2.1.2.2.1 (January, 1976)
AUTHOR: M. BAKKER.
INSTITUTE: MATHEMATICAL CENTRE, AMSTERDAM.
RECEIVED: 751231.
BRIEF DESCRIPTION:
THIS SECTION CONTAINS A PROCEDURE FOR THE SOLUTION
OF FOURTH ORDER SELF-ADJOINT LINEAR TWO POINT
BOUNDARY VALUE PROBLEMS;
FEM HERM SYM;
THIS PROCEDURE SOLVES THE DIFFERENTIAL EQUATION
(P(X)*Y'')'' - (Q(X)*Y')' + R(X)*Y = F(X), A < X < B,
WITH BOUNDARY CONDITIONS
Y(A) = E[1], Y'(A) = E[2],
Y(B) = E[3], Y'(B) = E[4].
KEY WORDS AND PHRASES:
FOURTH ORDER DIFFERENTIAL EQUATIONS,
TWO POINT BOUNDARY VALUE PROBLEMS,
SELF-ADJOINT BOUNDARY VALUE PROBLEMS,
GALERKIN'S METHOD,
DIRICHLET BOUNDARY CONDITIONS,
GLOBAL METHODS.
LANGUAGE: ALGOL 60.
REFERENCES:
[1] STRANG, G. AND G.J. FIX,
AN ANALYSIS OF THE FINITE ELEMENT METHOD,
PRENTICE-HALL, ENGLE WOOD 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).
CALLING SEQUENCE:
THE HEADING OF THE PROCEDURE READS:
"PROCEDURE" FEM HERM SYM(X, Y, N, P, Q, R, F, ORDER, E);
"VALUE" N, ORDER; "INTEGER" N, ORDER;
"ARRAY" X, Y, E;
"REAL" "PROCEDURE" P, Q, R, F;
"CODE" 33303;
THE MEANING OF THE FORMAL PARAMETERS IS:
N: <ARITHMETIC EXPRESSION>;
THE UPPER BOUND OF THE ARRAY X; 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[1:2*N-2];
EXIT: Y[2*I-1] IS AN APPROXIMATION TO Y(X[I]),
Y[2*I] IS AN APPROXIMATION TO Y'(X[I]),
WHERE Y(X) IS THE SOLUTION OF THE DIFFERENTIAL EQUATION
(1) (P(X)*Y'')'' - (Q(X)*Y')' + R(X)*Y = F(X) , A< X < B,
WITH BOUNDARY CONDITIONS
Y(A) = E[1], Y'(A) = E[2],
(2)
Y(B) = E[3], Y'(B) = E[4];
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);
P(X) SHOULD BE STRICTLY POSITIVE;
Q: <PROCEDURE IDENTIFIER>;
THE HEADING OF Q READS:
"REAL" "PROCEDURE" Q(X); "VALUE" X; "REAL" X;
Q(X) IS THE COEFFICIENT OF Y' IN (1);
Q(X) SHOULD BE NONNEGATIVE;
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);
R(X) SHOULD BE NONNEGATIVE;
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[2*I-1]-Y(X[I])) <= C1 * H**ORDER,
ABS(Y[2*I]-Y'(X[I]) <= C2 * H**ORDER, I = 1,...,N-1;
ORDER CAN ONLY BE CHOSEN EQUAL TO 4, 6, 8;
E: <ARRAY IDENTIFIER>;
"ARRAY" E[1:4];
E[1], ... , E[4] DESCRIBE THE BOUNDARY CONDITIONS (2).
PROCEDURES USED: CHLDECSOLBND = CP 34333
REQUIRED CENTRAL MEMORY:
ONE AUXILIARY ARRAY OF 8*(N-1) REALS IS USED.
RUNNING TIME:
LET K = ORDER/2; THEN
(A) K*N + 1 EVALUATIONS OF P(X), Q(X), R(X) AND F(X) ARE NEEDED;
(B) ABOUT (ORDER-3)*50*N MULTIPLICATIONS/DIVISIONS ARE NEEDED;
(C) ONE CALL OF CHLDECSOLBND IS DONE.
DATA AND RESULTS:
THE PROCEDURE FEM HERM SYM HAS SOME RESTRICTIONS:
(I) P(X) SHOULD BE POSITIVE ON THE CLOSED
INTERVAL <X[0],X[N]> AND Q(X) AND R(X) SHOULD BE
NONNEGATIVE THERE;
(II) P(X), Q(X), R(X) AND F(X) ARE REQUIRED TO BE SUFFICIENTLY
SMOOTH ON THE INTERVAL <X[0],X[N]> EXCEPT AT THE KNOTS,
WHERE DISCONTINUITIES OF THE DERIVATIVES ARE ALLOWED;
IN THAT CASE THE ORDER OF ACCURACY IS PRESERVED;
(III) THE USER SHOULD NOT EXPECT HIGHER ACCURACY THAN 12
DECIMALS DUE TO THE LOSS OF DIGITS DURING THE EVALUATION OF THE
MATRIX AND VECTOR AND DURING THE REDUCTION TO A PENTADIAGONAL
SYSTEM; THIS ACCURACY CAN BE REACHED VERY EASILY WHEN AN EIGTH
ORDER METHOD IS USED
METHOD AND PERFORMANCE:
PROBLEM (1)-(2) IS SOLVED BY MEANS OF GALERKIN'S METHOD WITH
CONTINUOUSLY DIFFERENTIABLE PIECEWISE POLYNOMIAL FUNCTIONS
(SEE [1], [2]) : THE SOLUTION IS APPROXIMATED BY A FUNCTION
WHICH IS CONTINUOUSLY DIFFERENTIABLE ON THE CLOSED INTERVAL
<X[0],X[N]> AND A POLYNOMIAL OF DEGREE LESS THAN OR EQUAL TO
K (K = 1 + ORDER//2) ON EACH CLOSED SEGMENT <X[J-1],X[J]>
(J = 1, ..., N);
THIS FUNCTION IS ENTIRELY DETERMINED BY THE VALUES OF
THE ZEROETH AND FIRST DERIVATIVE AT THE KNOTS X[J] AND BY
THE VALUES IT HAS AT (K-3) INTERIOR KNOTS ON EACH CLOSED
SEGMENT <X[J-1],X[J]>; THE VALUES OF THE FUNCTION AND ITS
DERIVATIVE AT THE KNOTS ARE OBTAINED BY THE SOLUTION OF AN
(ORDER + 1)-DIAGONAL LINEAR SYSTEM OF (K-1)*N - 2 UNKNOWNS;
THE ENTRIES OF THE MATRIX AND THE VECTOR ARE INNER PRODUCTS
WHICH ARE APPROXIMATED BY PIECEWISE K-POINT LOBATTO
QUADRATURE (SEE [3]); THE EVALUATION OF THE MATRIX AND
VECTOR IS PERFORMED SEGMENT BY SEGMENT;
IF K > 3 THE RESULTING LINEAR SYSTEM CAN BE REDUCED
TO A PENTADIAGONAL SYSTEM BY MEANS OF STATIC
CONDENSATION; THIS IS POSSIBLE BECAUSE THE FUNCTION
VALUES AT THE INTERIOR KNOTS ON EACH SEGMENT <X[J-1],X[J]>
DO NOT DEPEND ON FUNCTION VALUES OUTSIDE THAT SEGMENT;
THE FINAL PENTADIAGONAL SYSTEM, SINCE THE MATRIX IS POSITIVE
DEFINITE AND SYMMETRIC, IS SOLVED BY MEANS OF CHOLESKY'S
DECOMPOSITION METHOD (SEE SECTION 3.1.2.1.1.2.1.3).
EXAMPLE OF USE:
WE SOLVE THE BOUNDARY VALUE PROBLEM
WE SOLVE THE BOUNDARY VALUE PROBLEM
Y'''' - (Y'*COS(X))' + Y*EXP(X) = SIN(X)*(1 + EXP(X) + COS(X)*2),
0 < X < PI;
Y(0) = Y(PI) = 0; Y'(0) = 1; Y'(PI) = -1;
PI = 3.14159265358979;
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 = 5, 10 AND WE COMPUTE
THE MAXIMUM DEVIATIONS FROM Y(X[I]) AND Y'(X[I])
FOR ORDER = 4, 6, 8;
THE PROGRAM READS AS FOLLOWS:
"BEGIN" "INTEGER" N; "FOR" N:= 5, 10 "DO"
"BEGIN" "INTEGER" I, ORDER; "REAL" PI; "ARRAY" X[0:N],
Y[1:2*N-2], E[1:4];
"REAL" "PROCEDURE" P(X); "VALUE" X; "REAL" X; P:= 1;
"REAL" "PROCEDURE" Q(X); "VALUE" X; "REAL" X;
Q:= COS(X);
"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)+ 2*COS(X));
E[1]:= E[3]:= 0; E[2]:= 1; E[4]:= - 1;
PI:= 3.14159265358979;
"FOR" I:= 0 "STEP" 1 "UNTIL" N "DO" X[I]:= PI*I/N;
OUTPUT(61,"("//,6B"("N=")"ZD")",N);
"FOR" ORDER:= 4, 6, 8 "DO"
"BEGIN" "REAL" RHO1, RHO2, D1, D2;
FEM HERM SYM(X, Y, N, P, Q, R, F, ORDER, E);
RHO1:= RHO2:= 0;
"FOR" I:= 1 "STEP" 1 "UNTIL" N - 1 "DO"
"BEGIN" D1:= ABS(Y[2*I-1] - SIN(X[I]));
"IF" RHO1 < D1 "THEN" RHO1:= D1;
D2:= ABS(Y[2*I] - COS(X[I]));
"IF" RHO2 < D2 "THEN" RHO2:= D2
"END";
OUTPUT(61,"("/,16B"("ORDER=")"D,/,
24B"("MAX ABS(Y[2*I-1]-Y(X[I])) = ")",D.3D"+ZD,
/,24B"("MAX ABS(Y[2*I]-Y'(X[I])) = ")",D.3D"+ZD")",
ORDER,RHO1,RHO2)
"END"
"END"
"END"
RESULTS:
N= 5
ORDER=4
MAX ABS(Y[2*I-1]-Y(X[I])) = 4.822" -4
MAX ABS(Y[2*I]-Y'(X[I])) = 4.548" -4
ORDER=6
MAX ABS(Y[2*I-1]-Y(X[I])) = 5.651" -6
MAX ABS(Y[2*I]-Y'(X[I])) = 2.035" -6
ORDER=8
MAX ABS(Y[2*I-1]-Y(X[I])) = 2.264" -8
MAX ABS(Y[2*I]-Y'(X[I])) = 1.600" -8
N=10
ORDER=4
MAX ABS(Y[2*I-1]-Y(X[I])) = 2.657" -5
MAX ABS(Y[2*I]-Y'(X[I])) = 2.870" -5
ORDER=6
MAX ABS(Y[2*I-1]-Y(X[I])) = 8.398" -8
MAX ABS(Y[2*I]-Y'(X[I])) = 3.572" -8
ORDER=8
MAX ABS(Y[2*I-1]-Y(X[I])) = 7.981"-11
MAX ABS(Y[2*I]-Y'(X[I])) = 6.796"-11
NOTICE THAT THE MAXIMUM ERROR IS DIVIDED BY
2**ORDER, WHEN THE MESH SIZE IS HALVED.
SOURCE TEXT(S):
"CODE" 33303;