NUMAL Section 5.2.1.2.1.2.1.2
BEGIN SECTION : 5.2.1.2.1.2.1.2 (January, 1976)
AUTHOR: M. BAKKER.
INSTITUTE: MATHEMATICAL CENTRE, AMSTERDAM.
RECEIVED: 751231.
BRIEF DESCRIPTION:
THIS SECTION CONTAINS A PROCEDURE FOR THE SOLUTION
OF SECOND ORDER SKEW-ADJOINT LINEAR TWO POINT
BOUNDARY VALUE PROBLEMS;
FEM LAG SKEW;
THIS PROCEDURE SOLVES THE DIFFERENTIAL EQUATION
- Y'' + Q(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].
KEY WORDS AND PHRASES:
SECOND ORDER DIFFERENTIAL EQUATIONS,
TWO POINT BOUNDARY VALUE PROBLEMS,
SKEW-ADJOINT BOUNDARY VALUE PROBLEMS,
GALERKIN'S METHOD,
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 SKEW.
CALLING SEQUENCE:
THE HEADING OF THE PROCEDURE READS:
"PROCEDURE" FEM LAG SKEW(X, Y, N, Q, R, F, ORDER, E);
"VALUE" N, ORDER; "INTEGER" N, ORDER;
"ARRAY" X, Y, E;
"REAL" "PROCEDURE" Q, R, F;
"CODE" 33302;
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] TO THE DIFFERENTIAL EQUATION
(1) - Y'' + Q(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];
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);
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 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 Q(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 SKEW HAS SOME RESTRICTIONS IN ITS USE:
(I) Q(X) IS NOT ALLOWED TO HAVE VERY LARGE VALUES IN SOME SENSE:
THE PRODUCT Q(X)*(X[J] - X[J-1]) SHOULD NOT BE TOO LARGE
ON THE CLOSED INTERVAL <X[J-1],X[J]>, OTHERWISE
THE BOUNDARY VALUE PROBLEM MAY DEGENERATE TO A SINGULAR
PERTURBATION OR BOUNDARY LAYER PROBLEM, FOR WHICH EITHER
SPECIAL METHODS OR A SUITABLY CHOSEN GRID ARE NEEDED;
(II) Q(X), R(X) AND F(X) ARE REQUIRED TO BE SUFFICIENTLY
DIFFERENTIABLE ON THE DOMAIN OF THE BOUNDARY VALUE PROBLEM;
THEY ARE, HOWEVER, THE DERIVATIVES ARE ALLOWED TO HAVE
DISCONTINUITIES AT THE GRID POINTS, IN WHICH CASE THE ORDER OF
ACCURACY (2, 4 OR 6) IS PRESERVED;
(III) IF Q(X) AND R(X) SATISFY THE INEQUALITY R(X) >= Q'(X)/2,
THE EXISTENCE OF A UNIQUE SOLUTION IS GUARANTEED, OTHERWISE
THIS 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-DIGITS ACCURACY CAN BE OBTAINED
WITH A MODERATE MESH SIZE (E.G. < 0.1) ALREADY, PROVIDED
A SIXTH ORDER METHOD IS USED.
METHOD AND PERFORMANCE:
PROBLEM (1)-(2) IS SOLVED BY MEANS OF GALERKIN'S METHOD WITH
CONTINUOUS PIECEWISE POLYNOMIAL FUNCTIONS (SEE [1], [2]);
THE SOLUTION IS APPROXIMATED BY A FUNCTION WHICH IS CONTINUOUS
ON THE 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]); SINCE THE FINAL TRIDIAGONAL SYSTEM
IS OF FINITE DIFFERENCE TYPE, IT IS SOLVED BY MEANS
OF BABUSKA'S METHOD (SEE [4]).
EXAMPLE OF USE:
WE SOLVE THE BOUNDARY VALUE PROBLEM
- Y'' + Y'*COS(X) + Y*EXP(X) = SIN(X)*(1 + EXP(X)) + COS(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" 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)) + COS(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 SKEW(X, Y, N, Q, 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= 2.95" -3
ORDER=4 MAX. ERROR= 2.56" -5
ORDER=6 MAX. ERROR= 4.26" -8
N=20
ORDER=2 MAX. ERROR= 7.55" -4
ORDER=4 MAX. ERROR= 1.68" -6
ORDER=6 MAX. ERROR= 6.76"-10
NOTICE THAT THE MAXIMUM ERROR DECREASES BY ABOUT
2**(-ORDER) WHEN THE MESH SIZE IS HALVED.
SOURCE TEXT(S):
"CODE" 33302;