NUMAL Section 5.2.1.2.1.3
BEGIN SECTION : 5.2.1.2.1.3 (December, 1979)
AUTHOR: M. BAKKER.
INSTITUTE: MATHEMATICAL CENTRE, AMSTERDAM.
RECEIVED: 791231.
BRIEF DESCRIPTION:
THE PROCEDURE NONLIN FEMLAGSKEW SOLVES A NONLINEAR TWO POINT
BOUNDARY VALUE PROBLEM WITH SPHERICAL COORDINATES.
IT SOLVES THE DIFFERENTIAL EQUATION
(X**NC*Y')'/X**NC = F(X, Y, Y'), 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,
BOUNDARY VALUE PROBLEMS,
RITZ-GALERKIN METHOD,
SPHERICAL COORDINATES,
GLOBAL METHODS.
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 27, 1976.
[3] BABUSKA, I.,
NUMERICAL STABILITY IN PROBLEMS OF LINEAR ALGEBRA,
S.I.A.M. J. NUM. ANAL., VOL.9, P. 53-77 (1972).
[4] BAKKER, M.,
GALERKIN METHODS IN SPHERICAL REGIONS, TO APPEAR.
SUBSECTION: NONLIN FEM LAG SKEW.
CALLING SEQUENCE:
THE HEADING OF THE PROCEDURE READS:
"PROCEDURE" NONLIN FEM LAG SKEW(X, Y, N, F, FY, FZ, NC, E);
"INTEGER" N, NC;
"REAL" "PROCEDURE" F, FY, FZ;
"ARRAY" X, Y, E;
"CODE" 33314;
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 SEGMENT [A,B];
Y: <ARRAY IDENTIFIER>;
"ARRAY" Y[0:N];
ENTRY: Y[I] (I = 0, 1, ... , N) IS AN INITIAL APPROXIMATE
SOLUTION AT X[I] OF THE DIFFERENTIAL EQUATION
(3) (Y'*X**NC)'/X**NC = F(X,Y,Y') , 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];
EXIT: Y[I] (I = 0, 1, ... , N) IS THE GALERKIN
SOLUTION AT X[I] OF THE (3)-(4);
F: <PROCEDURE IDENTIFIER>;
THE HEADING OF F READS:
"REAL" "PROCEDURE" F(X,Y,Z); "VALUE" X,Y,Z; "REAL" X,Y,Z;
F(X,Y,Z) IS THE RIGHT HAND SIDE OF (3);
FY: <PROCEDURE IDENTIFIER>;
THE HEADING OF FY READS:
"REAL" "PROCEDURE" FY(X,Y,Z); "VALUE" X,Y,Z; "REAL" X,Y,Z;
FY(X,Y,Z) IS THE DERIVATIVE OF F WITH RESPECT TO Y;
FZ: <PROCEDURE IDENTIFIER>;
THE HEADING OF FZ READS:
"REAL" "PROCEDURE" FZ(X,Y,Z); "VALUE" X,Y,Z; "REAL" X,Y,Z;
FZ(X,Y,Z) IS THE DERIVATIVE OF F WITH RESPECT TO Z;
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: DUPVEC CP 31030.
REQUIRED CENTRAL MEMORY:
FIVE AUXILIARY ARRAYS OF N REALS ARE USED.
RUNNING TIME:
LET IT BE THE NUMBER OF NEWTON ITERATIONS; THEN
IT*N EVALUATIONS OF F, FY, FZ ARE NEEDED;
DATA AND RESULTS:
THE FUNCTIONS F, FY AND FZ ARE REQUIRED TO BE SUFFICIENTLY
SMOOTH IN THEIR VARIABLES ON THE INTERIOR OF EVERY SEGMENT
<X[I],X[I+1]> (I = 0, ..., N - 1);
METHOD AND PERFORMANCE:
LET Y[0](X) BE SOME INITIAL APPROXIMATION OF Y(X); THEN
THE NONLINEAR PROBLEM IS SOLVED BY SUCCESIVELY SOLVING
- (D[K]'*X**NC)'/X**NC
+ FY(X,Y[K](X),Y[K]'(X))*D[K](X)
+ FZ(X,Y[K](X),Y[K]'(X))*D[K]'(X)
= (Y[K]'*X**NC)'/X**NC
- F(X,Y[K],Y[K]'(X)), X[0] < X < X[N],
E[1]*D[K](X[0]) + E[2]*D[K]'(X[0]) = 0;
E[4]*D[K](X[N]) + E[5]*D[K]'(X[N]) = 0;
WITH GALERKIN'S METHOD (SEE PREVIOUS SECTION) AND PUTTING
Y[K+1](X) = Y[K](X) + D[K](X), K = 0,1,...
THIS IS THE SO-CALLED NEWTON-KANTOROWITCH METHOD;
EXAMPLE OF USE:
WE SOLVE THE BOUNDARY VALUE PROBLEM
(Y'*X**2)'/X**2 = EXP(Y)+EXP(Y')-EXP(1-X**2)-EXP(2*X)-6;
0 < X < 1, Y'(0) = Y(1) = 0;
FOR THE BOUNDARY CONDITIONS THIS MEANS THAT
E[2] = E[4] = 1; E[1] = E[3] = E[5] = E[6] = 0;
THE ANALYTIC SOLUTION IS Y(X) = 1 - X**2; WE APPROXIMATE
THE SOLUTION ON A UNIFORM GRID, I.E. X[I] = I/N,
I = 0, ..., N; WE CHOOSE N=10,20 AND COMPUTE THE MAXIMUM ERROR;
THE PROGRAM READS AS FOLLOWS:
"BEGIN" "INTEGER" N, NC;
"FOR" NC:= 0,1,2 "DO"
"FOR" N:= 25, 50 "DO"
"BEGIN" "INTEGER" I;"ARRAY" X, Y[0:N], E[1:6]; "REAL" RHO, D;
"REAL" "PROCEDURE" F(X,Y,Z); "VALUE" X,Y,Z; "REAL" X,Y,Z;
F:= EXP(Y)+EXP(Z)-EXP(1-X**2)-EXP(-2*X)-2-2*NC;
"REAL" "PROCEDURE" FY(X,Y,Z); "VALUE" X,Y,Z; "REAL" X,Y,Z;
FY:= EXP(Y);
"REAL" "PROCEDURE" FZ(X,Y,Z); "VALUE" X,Y,Z; "REAL" X,Y,Z;
FZ:= EXP(Z);
E[2]:= E[4]:= 1; E[1]:= E[3]:= E[5]:= E[6]:= 0;
"FOR" I:= 0 "STEP" 1 "UNTIL" N "DO"
"BEGIN" X[I]:= I/N; Y[I]: = 0 "END";
OUTPUT(61,"("//,4B"("N = ")"ZD,4B"("NC = ")"ZD")",N,NC);
NONLIN FEM LAG SKEW(X, Y, N, F, FY, FZ, NC, E);
RHO:= 0;
"FOR" I:= 0 "STEP" 1 "UNTIL" N "DO"
"BEGIN" D:= ABS(Y[I] - 1 + X[I]**2);
"IF" RHO < D "THEN" RHO:= D
"END";
OUTPUT(61,"("24B"("MAX.ERROR= ")",D.DD"+ZD")",RHO)
"END"
"END"
RESULTS:
N = 25 NC = 0 MAX.ERROR= 2.47" -4
N = 50 NC = 0 MAX.ERROR= 6.19" -5
N = 25 NC = 1 MAX.ERROR= 1.41" -3
N = 50 NC = 1 MAX.ERROR= 3.99" -4
N = 25 NC = 2 MAX.ERROR= 2.44" -3
N = 50 NC = 2 MAX.ERROR= 7.02" -4
ONE OBSERVES THAT THE MAXIMUM ERROR DECREASES BY ABOUT
0.25 WHEN THE MESH SIZE IS HALVED.
SOURCE TEXT(S):
"CODE" 33314;