function fem1d_adaptive ( )
%*****************************************************************************80
%
%% MAIN is the main program for FEM1D_ADAPTIVE.
%
% Discussion:
%
% FEM1D_ADAPTIVE solves a 1D problem using an adaptive finite element method.
%
% The equation to be treated is:
%
% -d/dx ( P(x) * dU(x)/dx ) + Q(x) * U(x) = F(x)
%
% by the finite-element method using piecewise linear basis
% functions.
%
% An adaptive method is used to try to reduce the maximum
% error by refining the mesh in certain places.
%
% Here U is an unknown scalar function of X defined on the
% interval [XL,XR], and P, Q and F are given functions of X.
%
% The values of U at XL and XR are also specified.
%
% The interval [XL,XR] is "meshed" with N+1 points,
%
% XN(0) = XL, XN(1) = XL+H, XN(2) = XL+2*H, ..., XN(N) = XR.
%
% This creates N subintervals, with interval I having endpoints
% XN(I-1) and XN(I).
%
% The algorithm tries to guarantee a certain amount
% of accuracy by examining the current solution, estimating the error
% in each subinterval, and, if necessary, subdividing one or more
% subintervals and repeating the calculation.
%
% We can think of the adaptive part of the algorithm as a refined
% problem. The program re-solves the problem on the pair of
% intervals J and J+1, which extend from node J-1 to node J+1.
% The values of U that were just computed at nodes J-1 and J+1
% will be used as the boundary values for this refined problem.
% The intervals J and J+1 will each be evenly divided into NY
% smaller subintervals. This boundary value problem is solved,
% and the derivatives of the original and refined solutions are
% then compared to get an estimate of the error.
%
% Licensing:
%
% This code is distributed under the GNU LGPL license.
%
% Modified:
%
% 06 November 2006
%
% Author:
%
% MATLAB version by John Burkardt
%
% Parameters:
%
% real ADIAG(NU).
% ADIAG(I) is the "diagonal" coefficient of the I-th
% equation in the linear system. That is, ADIAG(I) is
% the coefficient of the I-th unknown in the I-th equation.
%
% real ALEFT(NU).
% ALEFT(I) is the "left hand" coefficient of the I-th
% equation in the linear system. That is, ALEFT(I) is the
% coefficient of the (I-1)-th unknown in the I-th equation.
% There is no value in ALEFT(1), since the first equation
% does not refer to a "0-th" unknown.
%
% real ARITE(NU).
% ARITE(I) is the "right hand" coefficient of the I-th
% equation in the linear system. ARITE(I) is the coefficient
% of the (I+1)-th unknown in the I-th equation. There is
% no value in ARITE(NU) because the NU-th equation does not
% refer to an "NU+1"-th unknown.
%
% real ETA(N).
% ETA(I) is the error estimate for interval I. It is computed
% as the sum of two quantities, one associated with the left
% and one with the right node of the interval.
%
% real F(NU).
% ASSEMBLE stores into F the right hand side of the linear
% equations.
% SOLVE replaces those values of F by the solution of the
% linear equations.
%
% real FY(M).
% FY is the right hand side of the linear system of the refined
% problem.
%
% real H(N)
% H(I) is the length of subinterval I. This code uses
% equal spacing for all the subintervals.
%
% real HY(M).
% HY(I) is the length of subinterval I in the refined problem.
%
% integer IBC.
% IBC declares what the boundary conditions are.
% 1, at the left endpoint, U has the value UL,
% at the right endpoint, U' has the value UR.
% 2, at the left endpoint, U' has the value UL,
% at the right endpoint, U has the value UR.
% 3, at the left endpoint, U has the value UL,
% and at the right endpoint, U has the value UR.
% 4, at the left endpoint, U' has the value UL,
% at the right endpoint U' has the value UR.
%
% integer IBCY.
% IBCY declares the boundary conditions for the refined problem
% which should always be that the value of U is specified at
% both the left and right endpoints. This corresponds to a
% value of IBCY = 3.
%
% integer INDX(0:N).
% For a node I, INDX(I) is the index of the unknown
% associated with node I.
% If INDX(I) is equal to -1, then no unknown is associated
% with the node, because a boundary condition fixing the
% value of U has been applied at the node instead.
% Unknowns are numbered beginning with 1.
% If IBC is 2 or 4, then there is an unknown value of U
% at node 0, which will be unknown number 1. Otherwise,
% unknown number 1 will be associated with node 1.
% If IBC is 1 or 4, then there is an unknown value of U
% at node N, which will be unknown N or N+1,
% depending on whether there was an unknown at node 0.
%
% integer INDY(0:M).
% INDY(I) records the index of the unknown associated with
% node I for the refined problem.
%
% integer JADD(N).
% JADD(I) is 1 if the error estimates show that interval I
% should be subdivided.
%
% integer KOUNT, the number of adaptive steps that have been taken.
%
% integer M.
% M is the number of subintervals used in the refined problem.
% M is equal to NY for computations centered at node 0 or node N,
% and otherwise, M is equal to 2*NY.
%
% integer N
% The number of subintervals into which the interval
% [XL,XR] is broken.
%
% integer NL.
% The number of basis functions used in a single
% subinterval. (NL-1) is the degree of the polynomials
% used. For this code, NL is fixed at 2, meaning that
% piecewise linear functions are used as the basis.
%
% integer NMAX, the maximum number of unknowns that can be handled.
%
% integer NODE(NL,N).
% For each subinterval I:
% NODE(1,I) is the number of the left node, and
% NODE(2,I) is the number of the right node.
%
% integer NODEY(NL,M).
% NODEY performs the same function for the refined problem that
% NODE performs for the full problem, recording the node numbers
% associated with a particular subinterval.
%
% integer NQUAD
% The number of quadrature points used in a subinterval.
%
% integer NU.
% NU is the number of unknowns in the linear system.
% Depending on the value of IBC, there will be N-1,
% N, or N+1 unknown values, which are the coefficients
% of basis functions.
%
% integer NUY.
% The number of unknowns in the refined problem.
%
% integer NY.
% NY is the number of subintervals into which a given interval
% will be subdivided, before solving the refined probelm.
%
% integer PROBLEM, chooses the problem to be solved.
% The user must choose this value by setting it in routine GETPRB.
% * 1, u = x, p = 1, q = 0, f = 0, ibc = 3, ul = 0, ur = 1.
% The program should find the solution exactly, and the
% adaptive code should find that there is no reason to
% subdivide any interval.
% * 2, u = x*x, p = 1, q = 0, f = -2, ibc = 3, ul = 0, ur = 1.
% This problem should find the solution exactly, and
% the adaptive code should again find there is nothing
% to do.
% *3, u = sin(pi*x/2), p = 1, q = 0, ibc = 3, f = 0.25*pi*pi*sin(pi*x/2),
% ul = 0, ur = 1.
% *4, u = cos(pi*x/2), p = 1, q = 0, ibc = 3, f = 0.25*pi*pi*cos(pi*x/2),
% ul = 1, ur = 0.
% *5: u = x**(beta+2)/((beta+2)*(beta+1)), p = 1, q = 1, ibc = 3,
% f = -x**beta + (x**(beta+2))/((beta+2)*(beta+1)),
% ul = 0, ur = 1/((beta+2)*(beta+1))
% (beta must be greater than -2, and not equal to -1)
% *6: u = atan((x-0.5)/alpha), p = 1, q = 0, ibc = 3,
% f = 2*alpha*(x-0.5) / (alpha^2 + (x-0.5)^2) ^2,
% ul = u(0), ur = u(1)
%
% integer STATUS, reports status of subdivision.
% 0, a new subdivision was carried out.
% 1, no more subdivisions are needed.
% -1, no more subdivisions can be carried out.
%
% real TOL.
% A tolerance that is used to determine whether the estimated
% error in an interval is so large that it should be subdivided
% and the problem solved again.
%
% real UL.
% If IBC is 1 or 3, UL is the value that U is required
% to have at X = XL.
% If IBC is 2 or 4, UL is the value that U' is required
% to have at X = XL.
%
% real UR.
% If IBC is 2 or 3, UR is the value that U is required
% to have at X = XR.
% If IBC is 1 or 4, UR is the value that U' is required
% to have at X = XR.
%
% real WQUAD(NQUAD).
% WQUAD(I) is the weight associated with the I-th point
% of an NQUAD point Gaussian quadrature rule.
%
% real XL.
% XL is the left endpoint of the interval over which the
% differential equation is being solved.
%
% real XN(0:N).
% XN(I) is the location of the I-th node. XN(0) is XL,
% and XN(N) is XR.
%
% real XQUAD(NQUAD,NMAX), the I-th quadrature point
% in interval J.
%
% real XQUADY(NQUAD,NMAY ), the I-th quadrature point
% in subinterval J of the refined problem.
%
% real XR.
% XR is the right endpoint of the interval over which the
% differential equation is being solved.
%
% Workspace, double precision XT(0:NMAX), used to compute a new
% set of nodes.
%
% real YN(0:M).
% YN(I) is the location of the I-th node in the refined
% problem.
%
nl = 2;
nmax = 30;
nquad = 2;
timestamp ( );
fprintf ( 1, '\n' );
fprintf ( 1, 'FEM1D_ADAPTIVE\n' );
fprintf ( 1, ' MATLAB version\n' );
fprintf ( 1, '\n' );
fprintf ( 1, 'Solve the two-point boundary value problem:\n' );
fprintf ( 1, '\n' );
fprintf ( 1, ' -d/dx ( P(x) * dU(x)/dx ) + Q(x) * U(x) = F(x)\n' );
fprintf ( 1, '\n' );
fprintf ( 1, 'on the interval [0,1], specifying the value\n' );
fprintf ( 1, 'of U at each endpoint.\n' );
fprintf ( 1, '\n' );
fprintf ( 1, 'The number of basis functions per element is %d\n', nl );
fprintf ( 1, '\n' );
fprintf ( 1, 'The number of quadrature points per element is %d\n', nquad );
problem = get_problem ( );
fprintf ( 1, '\n' );
fprintf ( 1, ' Problem index = %d\n', problem );
fprintf ( 1, '\n' );
if ( problem == 1 )
fprintf ( 1, '\n' );
fprintf ( 1, ' "Linear" problem:\n' );
fprintf ( 1, ' (No refinement needed)\n' );
fprintf ( 1, '\n' );
fprintf ( 1, ' U(X) = X\n' );
fprintf ( 1, ' P(X) = 1.0\n' );
fprintf ( 1, ' Q(X) = 0.0\n' );
fprintf ( 1, ' F(X) = 0.0\n' );
fprintf ( 1, ' IBC = 3\n' );
fprintf ( 1, ' UL = 0.0\n' );
fprintf ( 1, ' UR = 1.0\n' );
elseif ( problem == 2 )
fprintf ( 1, '\n' );
fprintf ( 1, ' "Quadratic" problem:\n' );
fprintf ( 1, ' (No refinement needed)\n' );
fprintf ( 1, '\n' );
fprintf ( 1, ' U(X) = X*X\n' );
fprintf ( 1, ' P(X) = 1.0\n' );
fprintf ( 1, ' Q(X) = 0.0\n' );
fprintf ( 1, ' F(X) = -2.0\n' );
fprintf ( 1, ' IBC = 3\n' );
fprintf ( 1, ' UL = 0.0\n' );
fprintf ( 1, ' UR = 1.0\n' );
elseif ( problem == 3 )
fprintf ( 1, '\n' );
fprintf ( 1, ' "SINE" problem:\n' );
fprintf ( 1, '\n' );
fprintf ( 1, ' U(X) = SIN(PI*X/2)\n' );
fprintf ( 1, ' P(X) = 1.0\n' );
fprintf ( 1, ' Q(X) = 0.0\n' );
fprintf ( 1, ' F(X) = PI*PI*SIN(PI*X/2)/4\n' );
fprintf ( 1, ' IBC = 3\n' );
fprintf ( 1, ' UL = 0.0\n' );
fprintf ( 1, ' UR = 1.0\n' );
elseif ( problem == 4 )
fprintf ( 1, '\n' );
fprintf ( 1, ' "COSINE" problem:\n' );
fprintf ( 1, '\n' );
fprintf ( 1, ' U(X) = COS(PI*X/2)\n' );
fprintf ( 1, ' P(X) = 1.0\n' );
fprintf ( 1, ' Q(X) = 0.0\n' );
fprintf ( 1, ' F(X) = PI*PI*COS(PI*X/2)/4\n' );
fprintf ( 1, ' IBC = 3\n' );
fprintf ( 1, ' UL = 0.0\n' );
fprintf ( 1, ' UR = 1.0\n' );
elseif ( problem == 5 )
beta = get_beta ( );
fprintf ( 1, '\n' );
fprintf ( 1, ' "RHEINBOLDT" problem:\n' );
fprintf ( 1, '\n' );
fprintf ( 1, ' U(X) = X**(B+2)/((B+2)*(B+1))\n' );
fprintf ( 1, ' P(X) = 1.0\n' );
fprintf ( 1, ' Q(X) = 1.0\n' );
fprintf ( 1, ' F(X) = -X**B+(X**B+2))/((B+2)*(B+1))\n' );
fprintf ( 1, ' IBC = 3\n' );
fprintf ( 1, ' UL = 0.0\n' );
fprintf ( 1, ' UR = 1/((B+2)*(B+1))\n' );
fprintf ( 1, ' B = %f\n', beta );
elseif ( problem == 6 )
alpha = get_alpha ( );
fprintf ( 1, '\n' );
fprintf ( 1, ' "ARCTAN" problem:\n' );
fprintf ( 1, '\n' );
fprintf ( 1, ' U(X) = ATAN((X-0.5)/A)\n' );
fprintf ( 1, ' P(X) = 1.0\n' );
fprintf ( 1, ' Q(X) = 0.0\n' );
fprintf ( 1, ' F(X) = 2*A*(X-0.5)/(A^2+(X-0.5)^2)^2\n' );
fprintf ( 1, ' IBC = 3\n' );
fprintf ( 1, ' UL = ATAN(-0.5/A)\n' );
fprintf ( 1, ' UR = ATAN( 0.5/A)\n' );
fprintf ( 1, ' A = %f\n', alpha );
end
%
% Start out with just 4 subintervals.
%
n = 4;
%
% Initialize values that define the problem.
%
[ ibc, tol, ul, ur, xl, xn, xr ] = init ( n );
%
% Start the iteration counter off at 0.
%
kount = 0;
%
% Begin the next iteration.
%
while ( 1 )
kount = kount + 1;
fprintf ( 1, '\n' );
fprintf ( 1, 'Begin new iteration with %d nodes.\n', n );
fprintf ( 1, '\n' );
%
% Solve the regular problem.
%
[ u, h, nu ] = solvex ( ibc, kount, n, nl, nmax, ul, ur, xn );
%
% Solve N subproblems to get the error estimators.
%
eta = solvey ( u, h, n, nu, ul, ur, xn );
%
% Examine the error estimators, and see how many intervals should
% be subdivided.
%
[ n, xn, status ] = subdiv ( eta, kount, n, nmax, tol, xn );
if ( status ~= 0 )
break
end
end
%
% Terminate.
%
fprintf ( 1, '\n' );
fprintf ( 1, 'FEM1D_ADAPTIVE:\n' );
fprintf ( 1, ' Normal end of execution.\n' );
fprintf ( 1, '\n' );
timestamp ( );
return
end
function [ adiag, aleft, arite, f ] = assemble ( h, n, indx, node, nu, nl, ...
nquad, nmax, ul, ur, wquad, xn, xquad )
%*****************************************************************************80
%
%% ASSEMBLE assembles the global matrix.
%
% Licensing:
%
% This code is distributed under the GNU LGPL license.
%
% Modified:
%
% 05 November 2006
%
% Author:
%
% MATLAB version by John Burkardt
%
% Parameters:
%
% Input, real H(N)
% H(I) is the length of subinterval I. This code uses
% equal spacing for all the subintervals.
%
% Input, integer N
% The number of subintervals into which the interval
% [XL,XR] is broken.
%
% Input, integer INDX(0:N).
% For a node I, INDX(I) is the index of the unknown
% associated with node I.
% If INDX(I) is equal to -1, then no unknown is associated
% with the node, because a boundary condition fixing the
% value of U has been applied at the node instead.
% Unknowns are numbered beginning with 1.
% If IBC is 2 or 4, then there is an unknown value of U
% at node 0, which will be unknown number 1. Otherwise,
% unknown number 1 will be associated with node 1.
% If IBC is 1 or 4, then there is an unknown value of U
% at node N, which will be unknown N or N+1,
% depending on whether there was an unknown at node 0.
%
% Input, integer NODE(NL,N).
% For each subinterval I:
% NODE(1,I) is the number of the left node, and
% NODE(2,I) is the number of the right node.
%
% Input, integer NU.
% NU is the number of unknowns in the linear system.
% Depending on the value of IBC, there will be N-1,
% N, or N+1 unknown values, which are the coefficients
% of basis functions.
%
% Input, integer NL.
% The number of basis functions used in a single
% subinterval. (NL-1) is the degree of the polynomials
% used. For this code, NL is fixed at 2, meaning that
% piecewise linear functions are used as the basis.
%
% Input, integer NQUAD
% The number of quadrature points used in a subinterval.
%
% Input, integer NMAX, the maximum number of unknowns that can be handled.
%
% Input, real UL.
% If IBC is 1 or 3, UL is the value that U is required
% to have at X = XL.
% If IBC is 2 or 4, UL is the value that U' is required
% to have at X = XL.
%
% Input, real UR.
% If IBC is 2 or 3, UR is the value that U is required
% to have at X = XR.
% If IBC is 1 or 4, UR is the value that U' is required
% to have at X = XR.
%
% Input, real WQUAD(NQUAD).
% WQUAD(I) is the weight associated with the I-th point
% of an NQUAD point Gaussian quadrature rule.
%
% Input, real XN(0:N).
% XN(I) is the location of the I-th node. XN(0) is XL,
% and XN(N) is XR.
%
% Input, real XQUAD(NQUAD,NMAX), the I-th quadrature point
% in interval J.
%
% Output, real ADIAG(NU).
% ADIAG(I) is the "diagonal" coefficient of the I-th
% equation in the linear system. That is, ADIAG(I) is
% the coefficient of the I-th unknown in the I-th equation.
%
% Output, real ALEFT(NU).
% ALEFT(I) is the "left hand" coefficient of the I-th
% equation in the linear system. That is, ALEFT(I) is the
% coefficient of the (I-1)-th unknown in the I-th equation.
% There is no value in ALEFT(1), since the first equation
% does not refer to a "0-th" unknown.
%
% Output, real ARITE(NU).
% ARITE(I) is the "right hand" coefficient of the I-th
% equation in the linear system. ARITE(I) is the coefficient
% of the (I+1)-th unknown in the I-th equation. There is
% no value in ARITE(NU) because the NU-th equation does not
% refer to an "NU+1"-th unknown.
%
% Output, real F(NU).
% ASSEMBLE stores into F the right hand side of the linear
% equations.
% SOLVE replaces those values of F by the solution of the
% linear equations.
%
%
% Zero out the entries.
%
f(1:nu) = 0.0;
aleft(1:nu) = 0.0;
arite(1:nu) = 0.0;
adiag(1:nu) = 0.0;
%
% For each interval,
%
for ie = 1 : n
he = h(ie);
xleft = xn(node(1,ie)+1);
xrite = xn(node(2,ie)+1);
%
% For each quadrature point in the interval,
%
for iq = 1 : nquad
xquade = xquad(iq,ie);
wquade = wquad(iq);
%
% Pick a basis function which defines the equation,
%
for il = 1 : nl
ig = node(il,ie);
iu = indx(ig+1);
if ( 0 < iu )
[ phii, phiix ] = phi ( il, xquade, xleft, xrite );
f(iu) = f(iu) + he * wquade * ff ( xquade ) * phii;
%
% Take care of boundary conditions specifying the value of U'.
%
if ( ig == 0 )
x = 0.0;
f(iu) = f(iu) - pp ( x ) * ul;
elseif ( ig == n )
x = 1.0;
f(iu) = f(iu) + pp ( x ) * ur;
end
%
% Pick a basis function which defines the coefficient
% being computed.
%
for jl = 1 : nl
jg = node(jl,ie);
ju = indx(jg+1);
[ phij, phijx ] = phi ( jl, xquade, xleft, xrite );
aij = he * wquade * ...
( pp ( xquade ) * phiix * phijx ...
+ qq ( xquade ) * phii * phij );
%
% Decide where the coefficient is to be added.
%
if ( ju <= 0 )
if ( jg == 0 )
f(iu) = f(iu) - aij * ul;
elseif ( jg == n )
f(iu) = f(iu) - aij * ur;
end
elseif ( iu == ju )
adiag(iu) = adiag(iu) + aij;
elseif ( ju < iu )
aleft(iu) = aleft(iu) + aij;
else
arite(iu) = arite(iu) + aij;
end
end
end
end
end
end
return
end
function value = ff ( x )
%*****************************************************************************80
%
%% FF evaluates the function F in the differential equation.
%
% Discussion:
%
% This is the function F(X) that appears on the right hand
% side of the equation:
%
% -d/dx ( P(x) * dU(x)/dx ) + Q(x) * U(x) = F(x)
%
% Licensing:
%
% This code is distributed under the GNU LGPL license.
%
% Modified:
%
% 05 November 2006
%
% Author:
%
% MATLAB version by John Burkardt
%
% Parameters:
%
% Input, real X, the evaluation point.
%
% Output, real VALUE, the value of F(X).
%
%
% Find out which problem we're working on.
%
problem = get_problem ( );
if ( problem == 1 )
value = 0.0;
elseif ( problem == 2 )
value = -2.0 * x;
elseif ( problem == 3 )
value = 0.25 * pi^2 * sin ( 0.5 * pi * x );
elseif ( problem == 4 )
value = 0.25 * pi^2 * cos ( 0.5 * pi * x );
elseif ( problem == 5 )
beta = get_beta ( );
value = - ( x^beta ) + ( x^( beta + 2.0 ) ) ...
/ ( ( beta + 2.0 ) * ( beta + 1.0 ) );
elseif ( problem == 6 )
alpha = get_alpha ( );
value = 2.0 * alpha * ( x - 0.5 ) ...
/ ( alpha^2 + ( x - 0.5 )^2 )^2;
end
return
end
function [ h, indx, node, nu, wquad, xquad ] = geometry ( ibc, n, nl, nmax, ...
nquad, xn )
%*****************************************************************************80
%
%% GEOMETRY sets up some of the geometric information for the problem.
%
% Discussion:
%
% Note, however, that the location of the nodes
% is done outside of this routine, and, in fact, before this
% routine is called.
%
% Licensing:
%
% This code is distributed under the GNU LGPL license.
%
% Modified:
%
% 04 November 2006
%
% Author:
%
% MATLAB version by John Burkardt
%
% Parameters:
%
% Input, integer IBC.
% IBC declares what the boundary conditions are.
% 1, at the left endpoint, U has the value UL,
% at the right endpoint, U' has the value UR.
% 2, at the left endpoint, U' has the value UL,
% at the right endpoint, U has the value UR.
% 3, at the left endpoint, U has the value UL,
% and at the right endpoint, U has the value UR.
% 4, at the left endpoint, U' has the value UL,
% at the right endpoint U' has the value UR.
%
% Input, integer N
% The number of subintervals into which the interval
% [XL,XR] is broken.
%
% Input, integer NL.
% The number of basis functions used in a single
% subinterval. (NL-1) is the degree of the polynomials
% used. For this code, NL is fixed at 2, meaning that
% piecewise linear functions are used as the basis.
%
% Input, integer NMAX, the maximum number of unknowns that can be handled.
%
% Input, integer NQUAD
% The number of quadrature points used in a subinterval.
%
% Input, real XN(0:N).
% XN(I) is the location of the I-th node. XN(0) is XL,
% and XN(N) is XR.
%
% Output, real H(N)
% H(I) is the length of subinterval I. This code uses
% equal spacing for all the subintervals.
%
% Output, integer INDX(0:N).
% For a node I, INDX(I) is the index of the unknown
% associated with node I.
% If INDX(I) is equal to -1, then no unknown is associated
% with the node, because a boundary condition fixing the
% value of U has been applied at the node instead.
% Unknowns are numbered beginning with 1.
% If IBC is 2 or 4, then there is an unknown value of U
% at node 0, which will be unknown number 1. Otherwise,
% unknown number 1 will be associated with node 1.
% If IBC is 1 or 4, then there is an unknown value of U
% at node N, which will be unknown N or N+1,
% depending on whether there was an unknown at node 0.
%
% Output, integer NODE(NL,N).
% For each subinterval I:
% NODE(1,I) is the number of the left node, and
% NODE(2,I) is the number of the right node.
%
% Output, integer NU.
% NU is the number of unknowns in the linear system.
% Depending on the value of IBC, there will be N-1,
% N, or N+1 unknown values, which are the coefficients
% of basis functions.
%
% Output, real WQUAD(NQUAD).
% WQUAD(I) is the weight associated with the I-th point
% of an NQUAD point Gaussian quadrature rule.
%
% Output, real XQUAD(NQUAD,NMAX), the I-th quadrature point
% in interval J.
%
%
% Store in NODE the fact that interval I has node I-1
% as its left endpoint, and node I as its right endpoint.
%
for i = 1 : n
node(1,i) = i-1;
node(2,i) = i;
end
%
% For every node that is associated with an unknown, we
% record the number of the unknown in INDX.
%
nu = 0;
for i = 0 : n
if ( i == 0 & ( ibc == 1 | ibc == 3 ) )
indx(i+1) = -1;
elseif ( i == n & ( ibc == 2 | ibc == 3 ) )
indx(i+1) = -1;
else
nu = nu + 1;
indx(i+1) = nu;
end
end
%
% We compute the width of each interval.
%
for i = 1 : n
igl = node(1,i);
igr = node(2,i);
h(i) = xn(igr+1) - xn(igl+1);
end
%
% We compute the location of the quadrature points in each
% interval.
%
for i = 1 : n
xl = xn(node(1,i)+1);
xr = xn(node(2,i)+1);
if ( nquad == 1 )
xquad(1,i) = 0.5 * ( xl + xr );
elseif ( nquad == 2 )
alfa = -0.577350;
xquad(1,i) = ( ( 1.0 - alfa ) * xl ...
+ ( 1.0 + alfa ) * xr ) ...
/ 2.0;
alfa = +0.577350;
xquad(2,i) = ( ( 1.0 - alfa ) * xl ...
+ ( 1.0 + alfa ) * xr ) ...
/ 2.0;
elseif ( nquad == 3 )
alfa = -0.774597;
xquad(1,i) = ( ( 1.0 - alfa ) * xl ...
+ ( 1.0 + alfa ) * xr ) ...
/ 2.0;
xquad(2,i) = 0.5 * ( xl + xr );
alfa = +0.774597;
xquad(3,i) = ( ( 1.0 - alfa ) * xl ...
+ ( 1.0 + alfa ) * xr ) ...
/ 2.0;
end
end
%
% Store the weights for the quadrature rule.
%
if ( nquad == 1 )
wquad(1) = 1.0;
elseif ( nquad == 2 )
wquad(1) = 0.5;
wquad(2) = 0.5;
elseif ( nquad == 3 )
wquad(1) = 4.0 / 9.0;
wquad(2) = 5.0 / 18.0;
wquad(3) = 4.0 / 9.0;
end
return
end
function alpha = get_alpha ( )
%*****************************************************************************80
%
%% GET_ALPHA returns the value of ALPHA, for use by problem 6.
%
% Licensing:
%
% This code is distributed under the GNU LGPL license.
%
% Modified:
%
% 05 November 2006
%
% Author:
%
% MATLAB version by John Burkardt
%
% Parameters:
%
% Output, real ALPHA, the value of ALPHA.
%
alpha = 0.01;
return
end
function beta = get_beta ( )
%*****************************************************************************80
%
%% GET_BETA returns the value of BETA, for use by problem 5.
%
% Licensing:
%
% This code is distributed under the GNU LGPL license.
%
% Modified:
%
% 22 February 2014
%
% Author:
%
% MATLAB version by John Burkardt
%
% Parameters:
%
% Output, real BETA, the value of BETA.
%
beta = -0.9;
return
end
function problem = get_problem ( )
%*****************************************************************************80
%
%% GETPRB returns the value of the current problem number.
%
% Licensing:
%
% This code is distributed under the GNU LGPL license.
%
% Modified:
%
% 05 November 2006
%
% Author:
%
% MATLAB version by John Burkardt
%
% Parameters:
%
% Output, integer PROBLEM, the index of the problem.
%
problem = 6;
return
end
function [ ibc, tol, ul, ur, xl, xn, xr ] = init ( n )
%*****************************************************************************80
%
%% INIT initializes some parameters that define the problem.
%
% Licensing:
%
% This code is distributed under the GNU LGPL license.
%
% Modified:
%
% 06 November 2006
%
% Author:
%
% MATLAB version by John Burkardt
%
% Parameters:
%
% Input, integer N
% The number of subintervals into which the interval
% [XL,XR] is broken.
%
% Output, integer IBC.
% IBC declares what the boundary conditions are.
% 1, at the left endpoint, U has the value UL,
% at the right endpoint, U' has the value UR.
% 2, at the left endpoint, U' has the value UL,
% at the right endpoint, U has the value UR.
% 3, at the left endpoint, U has the value UL,
% and at the right endpoint, U has the value UR.
% 4, at the left endpoint, U' has the value UL,
% at the right endpoint U' has the value UR.
%
% Output, real TOL.
% A tolerance that is used to determine whether the estimated
% error in an interval is so large that it should be subdivided
% and the problem solved again.
%
% Output, real UL.
% If IBC is 1 or 3, UL is the value that U is required
% to have at X = XL.
% If IBC is 2 or 4, UL is the value that U' is required
% to have at X = XL.
%
% Output, real UR.
% If IBC is 2 or 3, UR is the value that U is required
% to have at X = XR.
% If IBC is 1 or 4, UR is the value that U' is required
% to have at X = XR.
%
% Output, real XL.
% XL is the left endpoint of the interval over which the
% differential equation is being solved.
%
% Output, real XN(0:N).
% XN(I) is the location of the I-th node. XN(0) is XL,
% and XN(N) is XR.
%
% Output, real XR.
% XR is the right endpoint of the interval over which the
% differential equation is being solved.
%
tol = 0.01;
%
% Find out which problem we're working on.
%
problem = get_problem ( );
%
% Set the boundary conditions for the problem, and
% print out its title.
%
if ( problem == 1 )
ibc = 3;
ul = 0.0;
ur = 1.0;
xl = 0.0;
xr = 1.0;
fprintf ( 1, '\n' );
fprintf ( 1, 'Exact solution is U = X\n' );
elseif ( problem == 2 )
ibc = 3;
ul = 0.0;
ur = 1.0;
xl = 0.0;
xr = 1.0;
fprintf ( 1, '\n' );
fprintf ( 1, 'Exact solution is U = X*X\n' );
elseif ( problem == 3 )
ibc = 3;
ul = 0.0;
ur = 1.0;
xl = 0.0;
xr = 1.0;
fprintf ( 1, '\n' );
fprintf ( 1, 'Exact solution is U = SIN(PI*X/2)\n' );
elseif ( problem == 4 )
ibc = 3;
ul = 1.0;
ur = 0.0;
xl = 0.0;
xr = 1.0;
fprintf ( 1, '\n' );
fprintf ( 1, 'Exact solution is U = COS(PI*X/2)\n' );
elseif ( problem == 5 )
ibc = 3;
beta = get_beta ( 'DUMMY' );
ul = 0.0;
ur = 1.0 / ( ( beta + 2.0 ) * ( beta + 1.0 ) );
xl = 0.0;
xr = 1.0;
fprintf ( 1, '\n' );
fprintf ( 1, 'Rheinboldt problem\n' );
elseif ( problem == 6 )
ibc = 3;
alpha = get_alpha ( );
xl = 0.0;
xr = 1.0;
ul = uexact ( xl );
ur = uexact ( xr );
fprintf ( 1, '\n' );
fprintf ( 1, 'Arctangent problem\n' );
end
%
% The nodes are defined here, and not in the geometry routine.
% This is because each new iteration chooses the location
% of the new nodes in a special way.
%
for i = 0 : n
xn(i+1) = ( ( n - i ) * xl ...
+ ( i ) * xr ) ...
/ ( n );
end
fprintf ( 1, 'The equation is to be solved for\n' );
fprintf ( 1, 'X greater than %f\n', xl );
fprintf ( 1, ' and less than %f\n', xr );
fprintf ( 1, '\n' );
fprintf ( 1, 'The boundary conditions are:\n' );
fprintf ( 1, '\n' );
if ( ibc == 1 | ibc == 3 )
fprintf ( 1, ' At X = XL, U = %f\n', ul );
else
fprintf ( 1, ' At X = XL, U'' = %f\n', ul );
end
if ( ibc == 2 | ibc == 3 )
fprintf ( 1, ' At X = XR, U = %f\n', ur );
else
fprintf ( 1, ' At X = XR, U'' %f\n= ', ur );
end
return
end
function output ( f, ibc, indx, n, nu, ul, ur, xn )
%*****************************************************************************80
%
%% OUTPUT prints out the computed solution.
%
% Licensing:
%
% This code is distributed under the GNU LGPL license.
%
% Modified:
%
% 06 November 2006
%
% Author:
%
% MATLAB version by John Burkardt
%
% Parameters:
%
% Input, real F(NU).
% ASSEMBLE stores into F the right hand side of the linear
% equations.
% SOLVE replaces those values of F by the solution of the
% linear equations.
%
% Input, integer IBC.
% IBC declares what the boundary conditions are.
% 1, at the left endpoint, U has the value UL,
% at the right endpoint, U' has the value UR.
% 2, at the left endpoint, U' has the value UL,
% at the right endpoint, U has the value UR.
% 3, at the left endpoint, U has the value UL,
% and at the right endpoint, U has the value UR.
% 4, at the left endpoint, U' has the value UL,
% at the right endpoint U' has the value UR.
%
% Input, integer INDX(0:N).
% For a node I, INDX(I) is the index of the unknown
% associated with node I.
% If INDX(I) is equal to -1, then no unknown is associated
% with the node, because a boundary condition fixing the
% value of U has been applied at the node instead.
% Unknowns are numbered beginning with 1.
% If IBC is 2 or 4, then there is an unknown value of U
% at node 0, which will be unknown number 1. Otherwise,
% unknown number 1 will be associated with node 1.
% If IBC is 1 or 4, then there is an unknown value of U
% at node N, which will be unknown N or N+1,
% depending on whether there was an unknown at node 0.
%
% Input, integer N
% The number of subintervals into which the interval
% [XL,XR] is broken.
%
% Input, integer NU.
% NU is the number of unknowns in the linear system.
% Depending on the value of IBC, there will be N-1,
% N, or N+1 unknown values, which are the coefficients
% of basis functions.
%
% Input, real UL.
% If IBC is 1 or 3, UL is the value that U is required
% to have at X = XL.
% If IBC is 2 or 4, UL is the value that U' is required
% to have at X = XL.
%
% Input, real UR.
% If IBC is 2 or 3, UR is the value that U is required
% to have at X = XR.
% If IBC is 1 or 4, UR is the value that U' is required
% to have at X = XR.
%
% Input, real XN(0:N).
% XN(I) is the location of the I-th node. XN(0) is XL,
% and XN(N) is XR.
%
fprintf ( 1, '\n' );
fprintf ( 1, 'Node X(I) U(X(I)) Uexact Error\n' );
fprintf ( 1, '\n' );
for i = 0 : n
if ( i == 0 )
if ( ibc == 1 | ibc == 3 )
u = ul;
else
u = f(indx(i+1));
end
elseif ( i == n )
if ( ibc == 2 | ibc == 3 )
u = ur;
else
u = f(indx(i+1));
end
else
u = f(indx(i+1));
end
uex = uexact ( xn(i+1) );
error = u - uex;
fprintf ( 1, ' %4d %12f %12f %12f %12f\n', i, xn(i+1), u, uex, error );
end
return
end
function [ phii, phiix ] = phi ( il, x, xleft, xrite )
%*****************************************************************************80
%
%% PHI evaluates a linear basis function and its derivative.
%
% Discussion:
%
% The functions are evaluated at a point X in an interval. In any
% interval, there are just two basis functions. The first
% basis function is a line which is 1 at the left endpoint
% and 0 at the right. The second basis function is 0 at
% the left endpoint and 1 at the right.
%
% Licensing:
%
% This code is distributed under the GNU LGPL license.
%
% Modified:
%
% 05 November 2006
%
% Author:
%
% MATLAB version by John Burkardt
%
% Parameters:
%
% Input, integer IL, the local index of the basis function.
%
% Input, real X, the evaluation point.
%
% Input, real XLEFT, XRITE, the endpoints of the interval.
%
% Output, real PHII, PHIIX, the value of the basis function
% and its derivative.
%
if ( xleft <= x & x <= xrite )
if ( il == 1 )
phii = ( xrite - x ) / ( xrite - xleft );
phiix = -1.0 / ( xrite - xleft );
else
phii = ( x - xleft ) / ( xrite - xleft );
phiix = 1.0 / ( xrite - xleft );
end
%
% If X is outside of the interval, then the basis function
% is always zero.
%
else
phii = 0.0;
phiix = 0.0;
end
return
end
function value = pp ( x )
%*****************************************************************************80
%
%% PP evaluates the function P in the differential equation.
%
% Discussion:
%
% The function P(X) occurs in the differential equation:
%
% -d/dx ( P(x) * dU(x)/dx ) + Q(x) * U(x) = F(x)
%
% Licensing:
%
% This code is distributed under the GNU LGPL license.
%
% Modified:
%
% 05 November 2006
%
% Author:
%
% MATLAB version by John Burkardt
%
% Parameters:
%
% Input, real X, the evaluation point.
%
% Output, real VALUE, the value of P(X).
%
%
% Find out which problem we're working on.
%
problem = get_problem ( );
if ( problem == 1 )
value = 1.0;
elseif ( problem == 2 )
value = 1.0;
elseif ( problem == 3 )
value = 1.0;
elseif ( problem == 4 )
value = 1.0;
elseif ( problem == 5 )
value = 1.0;
elseif ( problem == 6 )
value = 1.0;
end
return
end
function prsys ( adiag, aleft, arite, f, nu )
%*****************************************************************************80
%
%% PRSYS prints out the tridiagonal linear system.
%
% Licensing:
%
% This code is distributed under the GNU LGPL license.
%
% Modified:
%
% 05 November 2006
%
% Author:
%
% MATLAB version by John Burkardt
%
% Parameters:
%
% Input, real ADIAG(NU).
% ADIAG(I) is the "diagonal" coefficient of the I-th
% equation in the linear system. That is, ADIAG(I) is
% the coefficient of the I-th unknown in the I-th equation.
%
% Input, real ALEFT(NU).
% ALEFT(I) is the "left hand" coefficient of the I-th
% equation in the linear system. That is, ALEFT(I) is the
% coefficient of the (I-1)-th unknown in the I-th equation.
% There is no value in ALEFT(1), since the first equation
% does not refer to a "0-th" unknown.
%
% Input, real ARITE(NU).
% ARITE(I) is the "right hand" coefficient of the I-th
% equation in the linear system. ARITE(I) is the coefficient
% of the (I+1)-th unknown in the I-th equation. There is
% no value in ARITE(NU) because the NU-th equation does not
% refer to an "NU+1"-th unknown.
%
% Input, real F(NU).
% ASSEMBLE stores into F the right hand side of the linear
% equations.
% SOLVE replaces those values of F by the solution of the
% linear equations.
%
% integer NU.
% NU is the number of unknowns in the linear system.
% Depending on the value of IBC, there will be N-1,
% N, or N+1 unknown values, which are the coefficients
% of basis functions.
%
fprintf ( 1, '\n' );
fprintf ( 1, 'Printout of tridiagonal linear system:\n' );
fprintf ( 1, '\n' );
fprintf ( 1, 'Equation A-Left A-Diag A-Rite RHS\n' );
fprintf ( 1, '\n' );
for i = 1 : nu
if ( i == 1 )
fprintf ( 1, '%3d %12f %12f %12f\n', ...
i, adiag(i), arite(i), f(i) );
elseif ( i < nu )
fprintf ( 1, '%3d %12f %12f %12f %12f\n', ...
i, aleft(i), adiag(i), arite(i), f(i) );
else
fprintf ( 1, '%3d %12f %12f %12f\n', ...
i, aleft(i), adiag(i), f(i) );
end
end
return
end
function value = qq ( x )
%*****************************************************************************80
%
%% QQ evaluates the function Q in the differential equation.
%
% Discussion:
%
% The function Q(X) occurs in the differential equation:
%
% -d/dx ( P(x) * dU(x)/dx ) + Q(x) * U(x) = F(x)
%
% Licensing:
%
% This code is distributed under the GNU LGPL license.
%
% Modified:
%
% 05 November 2006
%
% Author:
%
% MATLAB version by John Burkardt
%
% Parameters:
%
% Input, real X, the evaluation point.
%
% Output, real VALUE, the value of Q(X).
%
%
% Find out which problem we're working on.
%
problem = get_problem ( );
if ( problem == 1 )
value = 0.0;
elseif ( problem == 2 )
value = 0.0;
elseif ( problem == 3 )
value = 0.0;
elseif ( problem == 4 )
value = 0.0;
elseif ( problem == 5 )
value = 1.0;
elseif ( problem == 6 )
value = 0.0;
end
return
end
function u = solve ( adiag, aleft, arite, f, nu )
%*****************************************************************************80
%
%% SOLVE solves a tridiagonal matrix system of the form A*x = b.
%
% Licensing:
%
% This code is distributed under the GNU LGPL license.
%
% Modified:
%
% 06 November 2006
%
% Author:
%
% MATLAB version by John Burkardt
%
% Parameters:
%
% Input, real ADIAG(NU), ALEFT(NU), ARITE(NU).
% the diagonal, left and right entries of the equations.
%
% Input, real F(NU), the right hand side of the linear system.
%
% Input, INTEGER NU, the number of equations to be solved.
%
% Output, real U(NU), the solution of the linear system.
%
%
% Handle the special case of a single equation.
%
if ( nu == 1 )
u(1) = f(1) / adiag(1);
%
% The general case, when NU is greater than 1.
%
else
arite(1) = arite(1) / adiag(1);
for i = 2 : nu-1
adiag(i) = adiag(i) - aleft(i) * arite(i-1);
arite(i) = arite(i) / adiag(i);
end
adiag(nu) = adiag(nu) - aleft(nu) * arite(nu-1);
u(1) = f(1) / adiag(1);
for i = 2 : nu
u(i) = ( f(i) - aleft(i) * u(i-1) ) / adiag(i);
end
for i = nu-1 : -1 : 1
u(i) = u(i) - arite(i) * u(i+1);
end
end
return
end
function [ u, h, nu ] = solvex ( ibc, kount, n, nl, nmax, ul, ur, xn )
%*****************************************************************************80
%
%% SOLVEX discretizes and solves a differential equation given the nodes.
%
% Licensing:
%
% This code is distributed under the GNU LGPL license.
%
% Modified:
%
% 05 November 2006
%
% Author:
%
% MATLAB version by John Burkardt
%
% Parameters:
%
% Input, integer IBC.
% IBC declares what the boundary conditions are.
% 1, at the left endpoint, U has the value UL,
% at the right endpoint, U' has the value UR.
% 2, at the left endpoint, U' has the value UL,
% at the right endpoint, U has the value UR.
% 3, at the left endpoint, U has the value UL,
% and at the right endpoint, U has the value UR.
% 4, at the left endpoint, U' has the value UL,
% at the right endpoint U' has the value UR.
%
% Input, integer KOUNT, the number of adaptive steps that have been taken.
%
% Input, integer N
% The number of subintervals into which the interval
% [XL,XR] is broken.
%
% Input, integer NL.
% The number of basis functions used in a single
% subinterval. (NL-1) is the degree of the polynomials
% used. For this code, NL is fixed at 2, meaning that
% piecewise linear functions are used as the basis.
%
% Input, integer NMAX, the maximum number of unknowns that can be handled.
%
% Input, real UL.
% If IBC is 1 or 3, UL is the value that U is required
% to have at X = XL.
% If IBC is 2 or 4, UL is the value that U' is required
% to have at X = XL.
%
% Input, real UR.
% If IBC is 2 or 3, UR is the value that U is required
% to have at X = XR.
% If IBC is 1 or 4, UR is the value that U' is required
% to have at X = XR.
%
% Input, real XN(0:N).
% XN(I) is the location of the I-th node. XN(0) is XL,
% and XN(N) is XR.
%
% Output, real U(NU), the finite element coefficients of the
% solution of the differential equation.
%
% Output, real H(N)
% H(I) is the length of subinterval I. This code uses
% equal spacing for all the subintervals.
%
% Output, integer NU.
% NU is the number of unknowns in the linear system.
% Depending on the value of IBC, there will be N-1,
% N, or N+1 unknown values, which are the coefficients
% of basis functions.
%
% Local parameters:
%
% Local, real ADIAG(NU).
% ADIAG(I) is the "diagonal" coefficient of the I-th
% equation in the linear system. That is, ADIAG(I) is
% the coefficient of the I-th unknown in the I-th equation.
%
% Local, real ALEFT(NU).
% ALEFT(I) is the "left hand" coefficient of the I-th
% equation in the linear system. That is, ALEFT(I) is the
% coefficient of the (I-1)-th unknown in the I-th equation.
% There is no value in ALEFT(1), since the first equation
% does not refer to a "0-th" unknown.
%
% Local, real ARITE(NU).
% ARITE(I) is the "right hand" coefficient of the I-th
% equation in the linear system. ARITE(I) is the coefficient
% of the (I+1)-th unknown in the I-th equation. There is
% no value in ARITE(NU) because the NU-th equation does not
% refer to an "NU+1"-th unknown.
%
% Local, integer INDX(0:N).
% For a node I, INDX(I) is the index of the unknown
% associated with node I.
% If INDX(I) is equal to -1, then no unknown is associated
% with the node, because a boundary condition fixing the
% value of U has been applied at the node instead.
% Unknowns are numbered beginning with 1.
% If IBC is 2 or 4, then there is an unknown value of U
% at node 0, which will be unknown number 1. Otherwise,
% unknown number 1 will be associated with node 1.
% If IBC is 1 or 4, then there is an unknown value of U
% at node N, which will be unknown N or N+1,
% depending on whether there was an unknown at node 0.
%
% Local, integer NODE(NL,N).
% For each subinterval I:
% NODE(1,I) is the number of the left node, and
% NODE(2,I) is the number of the right node.
%
% Local, integer NQUAD
% The number of quadrature points used in a subinterval.
%
% Local, real WQUAD(NQUAD).
% WQUAD(I) is the weight associated with the I-th point
% of an NQUAD point Gaussian quadrature rule.
%
% Local, real XQUAD(NQUAD,NMAX), the I-th quadrature point
% in interval J.
%
nquad = 2;
%
% Given a set of N nodes (where N increases on each iteration),
% compute the other geometric information.
%
[ h, indx, node, nu, wquad, xquad ] = geometry ( ibc, n, nl, nmax, nquad, xn );
%
% Assemble the linear system.
%
[ adiag, aleft, arite, f ] = assemble ( h, n, indx, node, nu, nl, ...
nquad, nmax, ul, ur, wquad, xn, xquad );
%
% Print out the linear system, just once.
%
if ( kount == 1 )
prsys ( adiag, aleft, arite, f, nu );
end
%
% Solve the linear system.
%
u = solve ( adiag, aleft, arite, f, nu );
%
% Print out the solution.
%
fprintf ( 1, '\n' );
fprintf ( 1, 'Basic solution\n' );
output ( u, ibc, indx, n, nu, ul, ur, xn );
return
end
function eta = solvey ( u, h, n, nu, ul, ur, xn )
%*****************************************************************************80
%
%% SOLVEY computes error estimators for a finite element solution.
%
% Discussion:
%
% SOLVEY accepts information about the solution of a finite element
% problem on a grid of nodes with coordinates XN. It then starts
% at node 0, and for each node, computes two "error estimators",
% one for the left, and one for the right interval associated with the
% node. These estimators are found by solving a finite element problem
% over the two intervals, using the known values of the original
% solution as boundary data, and using a mesh that is "slightly"
% refined over the original one.
%
% Note that the computations at the 0-th and N-th nodes only involve
% a single interval.
%
% Licensing:
%
% This code is distributed under the GNU LGPL license.
%
% Modified:
%
% 06 November 2006
%
% Author:
%
% MATLAB version by John Burkardt
%
% Parameters:
%
% Input, real U(NU), the finite element representation of the
% solution of the current problem.
%
% Input, real H(N)
% H(I) is the length of subinterval I. This code uses
% equal spacing for all the subintervals.
%
% Input, integer N
% The number of subintervals into which the interval
% [XL,XR] is broken.
%
% Input, integer NU.
% NU is the number of unknowns in the linear system.
% Depending on the value of IBC, there will be N-1,
% N, or N+1 unknown values, which are the coefficients
% of basis functions.
%
% Input, real UL.
% If IBC is 1 or 3, UL is the value that U is required
% to have at X = XL.
% If IBC is 2 or 4, UL is the value that U' is required
% to have at X = XL.
%
% Input, real UR.
% If IBC is 2 or 3, UR is the value that U is required
% to have at X = XR.
% If IBC is 1 or 4, UR is the value that U' is required
% to have at X = XR.
%
% Input, real XN(0:N).
% XN(I) is the location of the I-th node. XN(0) is XL,
% and XN(N) is XR.
%
% Output, real ETA(N).
% ETA(I) is the error estimate for interval I. It is computed
% as the sum of two quantities, one associated with the left
% and one with the right node of the interval.
%
nl = 2;
ny = 2;
nquad = 2;
nmay = 2 * ny;
%
% Initialize the error estimators to zero.
%
eta(1:n) = 0.0;
%
% Set the boundary conditions for each subproblem to be
% known values of U at the left and right.
%
%
% For each node, subdivide its left and right hand intervals
% into NY subintervals.
%
% Set up and solve the differential equation again on this
% smaller region.
%
% The 0-th and N-th nodes are special cases.
%
ibcy = 3;
for j = 0 : n
if ( j == 0 )
m = ny;
jlo = j;
jmid = j + 1;
jhi = j + 1;
elseif ( j == n )
m = ny;
jlo = j - 1;
jmid = j;
jhi = j;
else
m = 2 * ny;
jlo = j - 1;
jmid = j;
jhi = j + 1;
end
%
% Set the location of the nodes in the subintervals.
%
yl = xn(jlo+1);
ym = xn(jmid+1);
yr = xn(jhi+1);
for i = 0 : ny
yn(i+1) = ( ( ny - i ) * yl ...
+ ( i ) * ym ) ...
/ ( ny );
end
for i = ny+1 : m
yn(i+1) = ( ( m - i ) * ym ...
+ ( i - ny ) * yr ) ...
/ ( m - ny );
end
%
% Set up the geometry of the sub-problem.
%
[ hy, indy, nodey, nuy, wquad, xquady ] = geometry ( ibcy, m, nl, ...
nmay, nquad, yn );
%
% Set the boundary values for the sub-problem.
%
if ( j <= 1 )
uly = ul;
else
uly = u(j-1);
end
if ( n - 1 <= j )
ury = ur;
else
ury = u(j+1);
end
%
% Assemble the matrix for the sub-problem.
%
[ adiag, aleft, arite, fy ] = assemble ( hy, m, indy, nodey, nuy, nl, ...
nquad, nmay, uly, ury, wquad, yn, xquady );
%
% Solve the system.
%
uy = solve ( adiag, aleft, arite, fy, nuy );
%
% Compute the weighted sum of the squares of the differences
% of the original computed slope and the refined computed slopes.
%
% Calculation for left interval.
%
if ( 1 <= j )
if ( j <= 1 )
uleft = ul;
urite = u(1);
elseif ( j == n )
uleft = u(j-1);
urite = ur;
else
uleft = u(j-1);
urite = u(j);
end
uprime = ( urite - uleft ) / h(j);
total = 0.0;
for i = 1 : ny
yl = yn(i-1+1);
yr = yn(i+1);
if ( i == 1 )
vlval = uly;
vrval = uy(i);
elseif ( i == m )
vlval = uy(i-1);
vrval = ury;
else
vlval = uy(i-1);
vrval = uy(i);
end
vprime = ( vrval - vlval ) / hy(i);
ulval = ( real ( ny - i + 1 ) * uleft ...
+ real ( i - 1 ) * urite ) ...
/ real ( ny );
urval = ( real ( ny - i ) * uleft ...
+ real ( i ) * urite ) ...
/ real ( ny );
%
% Compute the integral of
%
% p(x)*(u'(x)-v'(x))^2 + q(x)*(u(x)-v(x))^2
%
for k = 1 : nquad
y = xquady(k,i);
uval = ( ( yl - y ) * urval ...
+ ( y - yr ) * ulval ) ...
/ ( yl - yr );
vval = ( ( yl - y ) * vrval ...
+ ( y - yr ) * vlval ) ...
/ ( yl - yr );
total = total + 0.5 * wquad(k) * hy(i) * ...
( pp ( y ) * ( uprime - vprime )^2 ...
+ qq ( y ) * ( uval - vval )^2 );
end
end
eta(j) = eta(j) + 0.5 * sqrt ( total );
end
%
% Calculation for right interval.
%
if ( j <= n - 1 )
if ( j == 0 )
uleft = ul;
urite = u(j+1);
elseif ( n - 1 <= j )
uleft = u(j);
urite = ur;
else
uleft = u(j);
urite = u(j+1);
end
uprime = ( urite - uleft ) / h(j+1);
total = 0.0;
for i = m+1-ny : m
yl = yn(i);
yr = yn(i+1);
if ( i == 1 )
vlval = uly;
vrval = uy(i);
elseif ( i == m )
vlval = uy(i-1);
vrval = ury;
else
vlval = uy(i-1);
vrval = uy(i);
end
vprime = ( vrval - vlval ) / hy(i);
ulval = ( ( m - i + 1 ) * uleft ...
+ ( ny - m + i - 1 ) * urite ) ...
/ ( ny );
urval = ( ( m - i ) * uleft ...
+ ( ny - m + i ) * urite ) ...
/ ( ny );
%
% Compute the integral of
%
% p(x)*(u'(x)-v'(x))^2 + q(x)*(u(x)-v(x))^2
%
for k = 1 : nquad
y = xquady(k,i);
uval = ( ( yl - y ) * urval ...
+ ( y - yr ) * ulval ) ...
/ ( yl - yr );
vval = ( ( yl - y ) * vrval ...
+ ( y - yr ) * vlval ) ...
/ ( yl - yr );
total = total + 0.5 * wquad(k) * hy(i) * ...
( pp ( y ) * ( uprime - vprime )^2 ...
+ qq ( y ) * ( uval - vval )^2 );
end
end
eta(j+1) = eta(j+1) + 0.5 * sqrt ( total );
end
end
%
% Print out the error estimators.
%
fprintf ( 1, '\n' );
fprintf ( 1, 'ETA\n' );
fprintf ( 1, '\n' );
for j = 1 : n
fprintf ( 1, ' %12f\n', eta(j) );
end
return
end
function [ n, xn, status ] = subdiv ( eta, kount, n, nmax, tol, xn )
%*****************************************************************************80
%
%% SUBDIV decides which intervals should be subdivided.
%
% Licensing:
%
% This code is distributed under the GNU LGPL license.
%
% Modified:
%
% 06 November 2006
%
% Author:
%
% MATLAB version by John Burkardt
%
% Parameters:
%
% Input, real ETA(N).
% ETA(I) is the error estimate for interval I. It is computed
% as the sum of two quantities, one associated with the left
% and one with the right node of the interval.
%
% Input, integer KOUNT, the number of adaptive steps that have been taken.
%
% Input, integer N
% The number of subintervals into which the interval
% [XL,XR] is broken.
%
% Input, integer NMAX, the maximum number of unknowns that can be handled.
%
% Input, real TOL.
% A tolerance that is used to determine whether the estimated
% error in an interval is so large that it should be subdivided
% and the problem solved again.
%
% Input, real XN(0:N).
% XN(I) is the location of the I-th node. XN(0) is XL,
% and XN(N) is XR.
%
% Output, integer N, the updated number of nodes.
%
% Output, real XN(0:N), the updated nodes.
%
% Output, integer STATUS, reports status of subdivision.
% 0, a new subdivision was carried out.
% 1, no more subdivisions are needed.
% -1, no more subdivisions can be carried out.
%
% Local Parameters:
%
% Local, integer JADD(N).
% JADD(I) is 1 if the error estimates show that interval I
% should be subdivided.
%
status = 0;
%
% Add up the ETA's, and get their average.
%
ave = sum ( eta(1:n) ) / n;
%
% Look for intervals whose ETA value is relatively large,
% and note in JADD that these intervals should be subdivided.
%
k = 0;
temp = max ( 1.2 * ave + 0.00001, tol^2 / n );
fprintf ( 1, '\n' );
fprintf ( 1, 'Tolerance = %f\n', temp );
fprintf ( 1, '\n' );
for j = 1 : n
if ( temp < eta(j) )
k = k + 1;
jadd(j) = 1;
fprintf ( 1, 'Subdivide interval %d\n', j );
else
jadd(j) = 0;
end
end
%
% If no subdivisions needed, we're done.
%
if ( k <= 0 )
fprintf ( 1, 'Success on step %d\n', kount );
status = 1;
return
end
%
% See if we're about to go over our limit.
%
if ( nmax < n + k )
fprintf ( 1, '\n' );
fprintf ( 1, 'The iterations did not reach their goal.\n' );
fprintf ( 1, 'The next value of N is %d\n', n + k );
fprintf ( 1, 'which exceeds NMAX = %d\n', nmax );
status = -1;
return
end
%
% Insert new nodes where needed.
%
k = 0;
xt(1) = xn(1);
for j = 1 : n
if ( 0 < jadd(j) )
xt(j+k+1) = 0.5 * ( xn(j+1) + xn(j) );
k = k + 1;
end
xt(j+k+1) = xn(j+1);
end
%
% Update the value of N, and copy the new nodes into XN.
%
n = n + k;
xn(1:n+1) = xt(1:n+1);
return
end
function timestamp ( )
%*****************************************************************************80
%
%% TIMESTAMP prints the current YMDHMS date as a timestamp.
%
% Licensing:
%
% This code is distributed under the GNU LGPL license.
%
% Modified:
%
% 14 February 2003
%
% Author:
%
% John Burkardt
%
t = now;
c = datevec ( t );
s = datestr ( c, 0 );
fprintf ( 1, '%s\n', s );
return
end
function value = uexact ( x )
%*****************************************************************************80
%
%% UEXACT returns the value of the exact solution at any point X.
%
% Licensing:
%
% This code is distributed under the GNU LGPL license.
%
% Modified:
%
% 05 November 2006
%
% Parameters:
%
% Input, real X, the evaluation point.
%
% Output, real VALUE, the value of the exact solution at X.
%
%
% Find out which problem we're working on.
%
problem = get_problem ( );
if ( problem == 1 )
value = x;
elseif ( problem == 2 )
value = x^2;
elseif ( problem == 3 )
value = sin ( pi * x / 2.0 );
elseif ( problem == 4 )
value = cos ( pi * x / 2.0 );
elseif ( problem == 5 )
beta = get_beta ( );
value = ( x^( beta + 2.0 ) ) ...
/ ( ( beta + 2.0 ) * ( beta + 1.0 ) );
elseif ( problem == 6 )
alpha = get_alpha ( );
value = atan ( ( x - 0.5 ) / alpha );
else
value = 0.0;
end
return
end自适应1D FEM
最新推荐文章于 2018-03-29 11:30:30 发布
1万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
