Finite Element Method with Adaptive Refinement
2024-05-25 19:38:00
FEM1D_ADAPTIVE is a MATLAB program which applies the finite element method to a linear two point boundary value problem in one spatial dimension, using adaptive refinement to estimate the error, refine the mesh, and produce an improved solution.
这里是原文链接
1.1D Time-Dependent Reaction-Diffusion Equation Neumann Boundary Conditions Finite Element Method
function fem1d ( )%*****************************************************************************80
%
%% MAIN is the main program for FEM1D.
%
% Discussion:
%
% FEM1D solves a one dimensional ODE using the finite element method.
%
% The differential equation has the form:
%
% -d/dx ( p(x) du/dx ) + q(x) * u = f(x)
%
% The finite-element method uses piecewise linear basis functions.
%
% 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 or U' at XL and XR are also specified.
%
% The interval [XL,XR] is "meshed" with NSUB+1 points,
%
% XN(0) = XL, XN(1)=XL+H, XN(2)=XL+2*H, ..., XN(NSUB)=XR.
%
% This creates NSUB subintervals, with interval number 1
% having endpoints XN(0) and XN(1), and so on up to interval
% NSUB, which has endpoints XN(NSUB-1) and XN(NSUB).
%
% Licensing:
%
% This code is distributed under the GNU LGPL license.
%
% Modified:
%
% 05 October 2008
%
% 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 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 H(N), the length of the subintervals.
%
% integer 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 INDX(1:N+1).
% 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 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 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 NQUAD.
% The number of quadrature points used in a subinterval.
% This code uses NQUAD = 1.
%
% integer NSUB, the number of subintervals into which the interval
% [XL,XR] is broken.
%
% integer NU, 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.
%
% 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 XL, the left endpoint of the interval over which the
% differential equation is being solved.
%
% real XN(1:N+1).
% XN(I) is the location of the I-th node. XN(1) is XL,
% and XN(N+1) is XR.
%
% real XQUAD(N)
% XQUAD(I) is the location of the single quadrature point
% in interval I.
%
% real XR, the right endpoint of the interval over which the
% differential equation is being solved.
%n = 100;nl = 2;timestamp ( );fprintf ( 1, '\n' );fprintf ( 1, 'FEM1D\n' );fprintf ( 1, ' MATLAB version\n' );fprintf ( 1, '\n' );fprintf ( 1, ' Solve the two-point boundary value problem:\n' );fprintf ( 1, ' -d/dx (p(x) du/dx) + q(x)*u = f(x)\n' );fprintf ( 1, ' on an interval [xl,xr], with the values of\n' );fprintf ( 1, ' u or u'' specified at xl and xr.\n' );fprintf ( 1, '\n' );fprintf ( 1, ' The interval is broken into %d subintervals.\n', n );fprintf ( 1, ' The number of basis functions per element is %d\n', nl );
%
% Initialize variables that define the problem.
%[ ibc, nquad, ul, ur, xl, xr ] = init ( );
%
% Compute the quantities that describe the geometry of the problem.
%[ h, indx, node, nu, xn, xquad ] = geometry ( ibc, nl, n, xl, xr );
%
% Assemble the matrix.
%[ adiag, aleft, arite, f ] = assemble ( h, indx, nl, node, ...nu, nquad, n, ul, ur, xn, xquad );
%
% Print out the linear system.
%system_print ( adiag, aleft, arite, f, nu );
%
% Solve the linear system.
%u = solve ( adiag, aleft, arite, f, nu );
%
% Print the current solution.
%output ( u, ibc, indx, n, nu, ul, ur, xn );
%
% Terminate.
%fprintf ( 1, '\n' );fprintf ( 1, 'FEM1D:\n' );fprintf ( 1, ' Normal end of execution.\n' );fprintf ( 1, '\n' );timestamp ( );
return
end
function [ adiag, aleft, arite, f ] = assemble ( h, indx, nl, ...node, nu, nquad, n, ul, ur, xn, xquad )%*****************************************************************************80
%
%% ASSEMBLE assembles the matrix and right hand side of the linear system.
%
% Discussion:
%
% Note that a 1 point quadrature rule, which is sometimes used to
% assemble the matrix and right hand side, is just barely accurate
% enough for simple problems. If you want better results, you
% should use a quadrature rule that is more accurate.
%
% Licensing:
%
% This code is distributed under the GNU LGPL license.
%
% Modified:
%
% 29 April 2007
%
% Author:
%
% MATLAB version by John Burkardt
%
% Parameters:
%
% Input, real H(N), the length of the subintervals.
%
% Input, integer INDX(1:N+1).
% 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 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 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 NQUAD, the number of quadrature points in a subinterval.
% This code uses NQUAD = 1.
%
% Input, integer N, the number of subintervals into which the interval
% [XL,XR] is broken.
%
% 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(1:N+1), the location of the I-th node. XN(1) is XL,
% and XN(N+1) is XR.
%
% Input, real XQUAD(N), the location of the single quadrature point
% in interval I.
%
% 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), the right hand side of the linear
% equations.
%f(1:nu) = 0.0;adiag(1:nu) = 0.0;aleft(1:nu) = 0.0;arite(1:nu) = 0.0;
%
% For element IE...
%for ie = 1 : nhe = h(ie);xleft = xn(node(1,ie)+1);xrite = xn(node(2,ie)+1);
%
% For quadrature point IQ...
%for iq = 1 : nquadxqe = xquad(ie);
%
% For basis function IL...
%for il = 1 : nlig = node(il,ie);iu = indx(ig+1);if ( 0 < iu )[ phii, phiix ] = phi ( il, xqe, xleft, xrite );f(iu) = f(iu) + he * ff ( xqe ) * phii;
%
% Handle boundary conditions.
%if ( ig == 0 )x = xn(1);f(iu) = f(iu) - pp ( x ) * ul;elseif ( ig == n )x = xn(n+1);f(iu) = f(iu) + pp ( x ) * ur;end
%
% For basis function JL...
%for jl = 1 : nljg = node(jl,ie);ju = indx(jg+1);[ phij, phijx ] = phi ( jl, xqe, xleft, xrite );aij = he * ( pp ( xqe ) * phiix * phijx ...+ qq ( xqe ) * phii * phij );if ( ju <= 0 )if ( jg == 0 )f(iu) = f(iu) - aij * ul;elseif ( jg == n )f(iu) = f(iu) - aij * ur;endelseif ( iu == ju )adiag(iu) = adiag(iu) + aij;elseif ( ju < iu )aleft(iu) = aleft(iu) + aij;elsearite(iu) = arite(iu) + aij;endendendendendend
return
end
function value = ff ( x )%*****************************************************************************80
%
%% FF returns the right hand side of the differential equation.
%
% Licensing:
%
% This code is distributed under the GNU LGPL license.
%
% Modified:
%
% 09 November 2006
%
% Author:
%
% MATLAB version by John Burkardt
%
% Parameters:
%
% Input, real X, the evaluation point.
%
% Output, real VALUE, the value of F(X).
%value = 0.0;
return
end
function [ h, indx, node, nu, xn, xquad ] = geometry ( ibc, nl, nsub, ...xl, xr )%*****************************************************************************80
%
%% GEOMETRY sets up the geometry for the interval [XL,XR].
%
% Licensing:
%
% This code is distributed under the GNU LGPL license.
%
% Modified:
%
% 09 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 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 NSUB.
% The number of subintervals into which the interval
% [XL,XR] is broken.
%
% Input, real XL.
% XL is the left endpoint of the interval over which the
% differential equation is being solved.
%
% Input, real XR.
% XR is the right endpoint of the interval over which the
% differential equation is being solved.
%
% Output, real H(NSUB), the length of the subintervals.
%
% Output, integer INDX(1:NSUB+1).
% 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,NSUB).
% 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 NSUB-1,
% NSUB, or NSUB+1 unknown values, which are the coefficients
% of basis functions.
%
% Output, real XN(1:NSUB+1).
% XN(I) is the location of the I-th node. XN(1) is XL,
% and XN(N+1) is XR.
%
% Output, real XQUAD(NSUB)
% XQUAD(I) is the location of the single quadrature point
% in interval I.
%%
% Set the value of XN, the locations of the nodes.
%xn = zeros ( nsub + 1, 1 );fprintf ( 1, '\n' );fprintf ( 1, ' Node Location\n' );fprintf ( 1, '\n' );for i = 0 : nsubxn(i+1) = ( ( nsub - i ) * xl ...+ ( i ) * xr ) .../ ( nsub );endr8vec_print_some ( nsub + 1, xn, 1, 10, ' First 10 nodes:' );
%
% Set the lengths of each subinterval.
%fprintf ( 1, '\n' );fprintf ( 1, 'Subint Length\n' );fprintf ( 1, '\n' );h = zeros ( nsub, 1 );for i = 1 : nsubh(i) = xn(i+1) - xn(i);endr8vec_print_some ( nsub, h, 1, 10, ' First 10 interval widths:' );
%
% Set the quadrature points, each of which is the midpoint of its subinterval.
%fprintf ( 1, '\n' );fprintf ( 1, 'Subint Quadrature point\n' );fprintf ( 1, '\n' );xquad = zeros ( nsub );for i = 1 : nsubxquad(i) = 0.5 * ( xn(i) + xn(i+1) );endr8vec_print_some ( nsub, xquad, 1, 10, ' First 10 quadrature points:' );
%
% Set the value of NODE, which records, for each interval,
% the node numbers at the left and right.
%node = zeros ( 2, nsub );for i = 1 : nsubnode(1,i) = i - 1;node(2,i) = i;endfprintf ( 1, '\n' );fprintf ( 1, ' First 10 pairs of nodes defining intervals:\n' );fprintf ( 1, 'Subint Left Node Right Node\n' );fprintf ( 1, '\n' );for i = 1 : min ( nsub, 10 )fprintf ( 1, ' %6d %6d %6d\n', i, node(1,i), node(2,i) );end
%
% Starting with node 0, see if an unknown is associated with
% the node. If so, give it an index.
%nu = 0;
%
% Handle first node.
%i = 0;if ( ibc == 1 || ibc == 3 )indx(i+1) = -1;elsenu = nu + 1;indx(i+1) = nu;end
%
% Handle nodes 1 through nsub-1
%for i = 1 : nsub-1nu = nu + 1;indx(i+1) = nu;end
%
% Handle the last node.
%i = nsub;if ( ibc == 2 || ibc == 3 )indx(i+1) = -1;elsenu = nu + 1;indx(i+1) = nu;endfprintf ( 1, '\n' );fprintf ( 1, ' First 10 unknown indices\n' );fprintf ( 1, ' Node Unknown\n' );fprintf ( 1, '\n' );for i = 0 : min ( nsub, 9 )fprintf ( 1, ' %6d %6d\n', i, indx(i+1) );end
return
end
function [ ibc, nquad, ul, ur, xl, xr ] = init ( )%*****************************************************************************80
%
%% INIT initializes variables that define the problem.
%
% Licensing:
%
% This code is distributed under the GNU LGPL license.
%
% Modified:
%
% 09 November 2006
%
% Author:
%
% MATLAB version by John Burkardt
%
% Parameters:
%
% 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, integer NQUAD.
% The number of quadrature points used in a subinterval.
% This code uses NQUAD = 1.
%
% 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 XR.
% XR is the right endpoint of the interval over which the
% differential equation is being solved.
%ibc = 1;nquad = 1;ul = 0.0;ur = 1.0;xl = 0.0;xr = 1.0;
%
% Print out the values that have been set.
%fprintf ( 1, '\n' );fprintf ( 1, 'The equation is to be solved for\n' );fprintf ( 1, 'X greater than XL = %f\n', xl );fprintf ( 1, ' and less than XR = %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 );elsefprintf ( 1, ' At X = XL, U'' = %f\n', ul );endif ( ibc == 2 || ibc == 3 )fprintf ( 1, ' At X = XR, U = %f\n', ur );elsefprintf ( 1, ' At X = XR, U'' = %f\n', ur );endfprintf ( 1, '\n' );fprintf ( 1, 'Number of quadrature points per element is %d\n', nquad );
return
end
function output ( f, ibc, indx, nsub, nu, ul, ur, xn )%*****************************************************************************80
%
%% OUTPUT prints out the computed solution at the nodes.
%
% Licensing:
%
% This code is distributed under the GNU LGPL license.
%
% Modified:
%
% 09 November 2006
%
% Author:
%
% MATLAB version by John Burkardt
%
% Parameters:
%
% Input, real F(NU), 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(1:N+1).
% 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 NSUB.
% 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(1:N+1).
% XN(I) is the location of the I-th node. XN(1) is XL,
% and XN(N+1) is XR.
%fprintf ( 1, '\n' );fprintf ( 1, 'First 10 entries of computed solution:\n' );fprintf ( 1, '\n' );fprintf ( 1, ' Node X(I) U(I)\n' );fprintf ( 1, '\n' );for i = 0 : min ( nsub, 9 )if ( i == 0 )if ( ibc == 1 || ibc == 3 )u = ul;elseu = f(indx(i+1));endelseif ( i == nsub )if ( ibc == 2 || ibc == 3 )u = ur;elseu = f(indx(i+1));endelseu = f(indx(i+1));endfprintf ( 1, ' %6d %12f %12f\n', i, xn(i+1), u );end
return
end
function [ phii, phiix ] = phi ( il, x, xleft, xrite )%*****************************************************************************80
%
%% PHI evaluates a linear basis function and its derivative.
%
% Discussion:
%
% 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:
%
% 09 November 2006
%
% Author:
%
% MATLAB version by John Burkardt
%
% Parameters:
%
% Input, integer IL, the index of the basis function.
% 1, the function which is 1 at XLEFT and 0 at XRITE.
% 2, the function which is 0 at XLEFT and 1 at XRITE.
%
% Input, real X, the evaluation point.
%
% Input, real XLEFT, XRITE, the left and right
% endpoints of the interval.
%
% Output, real PHII, PHIIX, the value of the
% basis function and its derivative at X.
%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.
%elsephii = 0.0;phiix = 0.0;end
return
end
function value = pp ( x )%*****************************************************************************80
%
%% PP returns the value of the coefficient function P(X).
%
% Licensing:
%
% This code is distributed under the GNU LGPL license.
%
% Modified:
%
% 09 November 2006
%
% Author:
%
% MATLAB version by John Burkardt
%
% Parameters:
%
% Input, real X, the evaluation point.
%
% Output, real VALUE, the value of P(X).
%value = 1.0;
return
end
function value = qq ( x )%*****************************************************************************80
%
%% QQ returns the value of the coefficient function Q(X).
%
% Licensing:
%
% This code is distributed under the GNU LGPL license.
%
% Modified:
%
% 09 November 2006
%
% Author:
%
% MATLAB version by John Burkardt
%
% Parameters:
%
% Input, real X, the evaluation point.
%
% Output, real VALUE, the value of Q(X).
%value = 0.0;
return
end
function r8vec_print_some ( n, a, i_lo, i_hi, title )%*****************************************************************************80
%
%% R8VEC_PRINT_SOME prints "some" of an R8VEC.
%
% Discussion:
%
% An R8VEC is a vector of R8 values.
%
% Licensing:
%
% This code is distributed under the GNU LGPL license.
%
% Modified:
%
% 10 September 2009
%
% Author:
%
% John Burkardt
%
% Parameters:
%
% Input, integer N, the dimension of the vector.
%
% Input, real A(N), the vector to be printed.
%
% Input, integer MAX_PRINT, the maximum number of lines to print.
%
% Input, string TITLE, a title.
%fprintf ( 1, '\n' );fprintf ( 1, '%s\n', title );fprintf ( 1, '\n' );for i = max ( 1, i_lo ) : min ( n, i_hi )fprintf ( 1, ' %8d: %12f\n', i, a(i) );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:
%
% 09 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.
% Note that for the first equation, there is no ALEFT
% coefficient, and for the last, there is no ARITE.
% So there is no need to store a value in ALEFT(1), nor
% in ARITE(NU).
%
% Input, real F(NU), the right hand side of the linear
% system to be solved.
%
% Input, integer NU.
% NU is the number of equations to be solved.
%
% Output, real U(NU), the solution of the linear system.
%arite(1) = arite(1) / adiag(1);for i = 2 : nu-1adiag(i) = adiag(i) - aleft(i) * arite(i-1);arite(i) = arite(i) / adiag(i);endadiag(nu) = adiag(nu) - aleft(nu) * arite(nu-1);u = zeros ( nu, 1 );u(1) = f(1) / adiag(1);for i = 2 : nuu(i) = ( f(i) - aleft(i) * u(i-1) ) / adiag(i);endfor i = nu-1 : -1 : 1u(i) = u(i) - arite(i) * u(i+1);end
return
end
function system_print ( adiag, aleft, arite, f, nu )%*****************************************************************************80
%
%% SYSTEM_PRINT prints out the tridiagonal linear system to be solved.
%
% Licensing:
%
% This code is distributed under the GNU LGPL license.
%
% Modified:
%
% 05 October 2008
%
% Author:
%
% 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.
%fprintf ( 1, '\n' );fprintf ( 1, 'First 10 rows of tridiagonal linear system:\n' );fprintf ( 1, '\n' );fprintf ( 1, 'Equation ALEFT ADIAG ARITE RHS\n' );fprintf ( 1, '\n' );for i = 1 : min ( nu, 10 )fprintf ( 1, '%3d', i );if ( i == 1 )fprintf ( 1, ' ' );elsefprintf ( 1, ' %12f', aleft(i) );endfprintf ( 1, ' %12f', adiag(i) );if ( i < nu )fprintf ( 1, ' %12f', arite(i) );elsefprintf ( 1, ' ' );endfprintf ( 1, ' %12f\n', f(i) );end
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
endfunction 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 )breakendend
%
% 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 : nhe = h(ie);xleft = xn(node(1,ie)+1);xrite = xn(node(2,ie)+1);
%
% For each quadrature point in the interval,
%for iq = 1 : nquadxquade = xquad(iq,ie);wquade = wquad(iq);
%
% Pick a basis function which defines the equation,
%for il = 1 : nlig = 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 : nljg = 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;endelseif ( iu == ju )adiag(iu) = adiag(iu) + aij;elseif ( ju < iu )aleft(iu) = aleft(iu) + aij;elsearite(iu) = arite(iu) + aij;endendendendendend
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 : nnode(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 : nif ( i == 0 & ( ibc == 1 | ibc == 3 ) )indx(i+1) = -1;elseif ( i == n & ( ibc == 2 | ibc == 3 ) )indx(i+1) = -1;elsenu = nu + 1;indx(i+1) = nu;endend
%
% We compute the width of each interval.
%for i = 1 : nigl = 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 : nxl = 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;endend
%
% 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 : nxn(i+1) = ( ( n - i ) * xl ...+ ( i ) * xr ) .../ ( n );endfprintf ( 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 );elsefprintf ( 1, ' At X = XL, U'' = %f\n', ul );endif ( ibc == 2 | ibc == 3 )fprintf ( 1, ' At X = XR, U = %f\n', ur );elsefprintf ( 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 : nif ( i == 0 )if ( ibc == 1 | ibc == 3 )u = ul;elseu = f(indx(i+1));endelseif ( i == n )if ( ibc == 2 | ibc == 3 )u = ur;elseu = f(indx(i+1));endelseu = f(indx(i+1));enduex = 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 );elsephii = ( x - xleft ) / ( xrite - xleft );phiix = 1.0 / ( xrite - xleft );end
%
% If X is outside of the interval, then the basis function
% is always zero.
%elsephii = 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 : nuif ( 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) );elsefprintf ( 1, '%3d %12f %12f %12f\n', ...i, aleft(i), adiag(i), f(i) );endend
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.
%elsearite(1) = arite(1) / adiag(1);for i = 2 : nu-1adiag(i) = adiag(i) - aleft(i) * arite(i-1);arite(i) = arite(i) / adiag(i);endadiag(nu) = adiag(nu) - aleft(nu) * arite(nu-1);u(1) = f(1) / adiag(1);for i = 2 : nuu(i) = ( f(i) - aleft(i) * u(i-1) ) / adiag(i);endfor i = nu-1 : -1 : 1u(i) = u(i) - arite(i) * u(i+1);endend
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 : nif ( j == 0 )m = ny;jlo = j;jmid = j + 1;jhi = j + 1;elseif ( j == n )m = ny;jlo = j - 1;jmid = j;jhi = j;elsem = 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 : nyyn(i+1) = ( ( ny - i ) * yl ...+ ( i ) * ym ) .../ ( ny );endfor i = ny+1 : myn(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;elseuly = u(j-1);endif ( n - 1 <= j )ury = ur;elseury = 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;elseuleft = u(j-1);urite = u(j);enduprime = ( urite - uleft ) / h(j);total = 0.0;for i = 1 : nyyl = 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;elsevlval = uy(i-1);vrval = uy(i);endvprime = ( 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 : nquady = 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 );endendeta(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;elseuleft = u(j);urite = u(j+1);enduprime = ( urite - uleft ) / h(j+1);total = 0.0;for i = m+1-ny : myl = yn(i);yr = yn(i+1);if ( i == 1 )vlval = uly;vrval = uy(i);elseif ( i == m )vlval = uy(i-1);vrval = ury;elsevlval = uy(i-1);vrval = uy(i);endvprime = ( 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 : nquady = 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 );endendeta(j+1) = eta(j+1) + 0.5 * sqrt ( total );endend
%
% Print out the error estimators.
%fprintf ( 1, '\n' );fprintf ( 1, 'ETA\n' );fprintf ( 1, '\n' );for j = 1 : nfprintf ( 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 : nif ( temp < eta(j) )k = k + 1;jadd(j) = 1;fprintf ( 1, 'Subdivide interval %d\n', j );elsejadd(j) = 0;endend
%
% If no subdivisions needed, we're done.
%if ( k <= 0 )fprintf ( 1, 'Success on step %d\n', kount );status = 1;
returnend
%
% 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;
returnend
%
% Insert new nodes where needed.
%k = 0;xt(1) = xn(1);for j = 1 : nif ( 0 < jadd(j) )xt(j+k+1) = 0.5 * ( xn(j+1) + xn(j) );k = k + 1;endxt(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 );elsevalue = 0.0;end
return
endFinite Element Method with Adaptive Refinement相关推荐
- 二阶偏微分方程组 龙格库塔法_有限单元法(Finite Element Method)实现声波方程模拟（Part 2）...
2.1 前言 承接上一篇文章,前面我们已经介绍了一维声波方程有限元求解: 蓝不是蓝:有限单元法(Finite Element Method)实现声波方程模拟(Part 1)zhuanlan.zhih ...
- 文献学习(part87)--Linearized Alternating Direction Method with Adaptive Penalty for Low-Rank ...
学习笔记,仅供参考,有错必纠 文章目录 Linearized Alternating Direction Method with Adaptive Penalty for Low-Rank Repre ...
- 关于元素定位使用class_name定位报错的部分问题Unable to locate element: {“method“:“css selector“
在使用selenium的定位web元素的时候,有时候会提示找不到定位元素,报错:selenium.common.exceptions.NoSuchElementException: Message: ...
- 有限差分法（Finite Difference Method）解方程：边界和内部结点的控制方程
本文转载自本人发布在简书平台上的文章 1. FDM解常微分方程 1.1. 问题描述 d 2 ϕ d x 2 = S ϕ (1) \frac{d^2\phi}{dx^2}=S_{\phi} \tag{1 ...
- 常微分方程中混合边界及狄利克雷边界条件下的Shooting method以及Finite difference method
声明 本部分是一个学习笔记,主要内容来自于华冬英老师编写的<微分方程的数值解法与程序实践>.如果觉得内容不错,可自行购买价格良心的官方正版教材.http://www.hxedu.com.c ...
- 【论文合集】2022年10月医学影像期刊论文合集
★ 本月IEEE Transactions on Medical Imaging(1区 top if 11.037) 共32篇, Medical Image Analysis(1区 top if 13 ...
- 《DeepXDE:a deep learning library for solving differential equations》梳理
本论文向我们介绍了一个求解微分方程的神经网络PINNs,和PINNs的Python库DeepXDE 摘要 1.PINNs(Physics-informed neural networks物理信息神经网 ...
- matlab 工具箱下载地址
1.平面操作工具箱 http://cathy.ijs.si/~leon/planman.html 2.SimMechanics 工具箱 (这个好像不是免费的) http://www.mathwork ...
- Computer Vision_2_Active Shape Models：Active Shape Models-Their Training and Application——1995
此为计算机视觉部分,主要侧重在底层特征提取,视频分析,跟踪,目标检测和识别方面等方面. 1. Active Appearance Models 活动表观模型和活动轮廓模型基本思想来源 Snake,现在 ...
最新文章
- GAN版马里奥创作家来了：一个样本即可训练，生成关卡要素丰富 | 开源
- 小程序 a标签_微慕WordPress小程序增强版V2.0新版上线
- 我所理解的Java NIO
- 重写系统中的UINavigationController返回按钮的事件
- python做什么项目好_推荐两个牛逼的Python项目
- 为什么说嵌入式开发比单片机要难很多？
- bom .dom_MicroProfile 2.2 BOM导入支持
- Java课程设计和sql数据库_数据库SQLserver+java课程设计
- go设置后端启动_Vue 之前后端分离的跨域
- 怎么设置php 中图片的大小写,php中强制字母转换大小写的方法有哪些
- 美容院管理系统高效管理门店店务？
- 第十二课，assimp模型加载（数据加载篇）
- 光棍.com市场推广策划书(爆笑)
- transformer中的位置嵌入
- 电脑打字习惯让人提笔忘字
- 查询最近12个月的数据SQL语句
- UART、I2C、USB、SPI、CAN、Jtag、PCI/PCIE协议汇总
- 流体动力学控制方程（详细推导）
- 6号团队-团队任务5:项目总结会
- Android应用全屏显示