# include
# include
# include
# include
int main ( );
void assemble ( double adiag[], double aleft[], double arite[], double f[],
double h[], int n, int indx[], int node[], int nu, int nl, int nquad,
int nmax, double ul, double ur, double wquad[], double xn[], double xquad[] );
double ff ( double x );
void geometry ( double h[], int ibc, int indx[], int n, int nl, int nmax,
int node[], int nquad, int *nu, double wquad[], double xn[], double xquad[] );
double get_alpha ( );
double get_beta ( );
int get_problem ( );
void init ( int *ibc, int n, double *tol, double *ul, double *ur, double *xl,
double xn[], double *xr );
void output ( double f[], int ibc, int indx[], int n, int nu, double ul,
double ur, double xn[] );
void phi ( int il, double x, double *phii, double *phiix, double xleft,
double xrite );
double pp ( double x );
void prsys ( double adiag[], double aleft[], double arite[], double f[],
int nu );
double qq ( double x );
void solve ( double adiag[], double aleft[], double arite[], double f[],
int nu );
void solvex ( double adiag[], double aleft[], double arite[], double f[],
double h[], int ibc, int indx[], int kount, int n, int nl, int nmax,
int node[], int nquad, int *nu, double ul, double ur, double wquad[],
double xn[], double xquad[] );
void solvey ( double eta[], double f[], double h[], int n, int nu, double ul,
double ur, double xn[] );
int subdiv ( double eta[], int kount, int *n, int nmax, double tol,
double xn[] );
void timestamp ( );
double uexact ( double x );
/******************************************************************************/
int main ( )
/******************************************************************************/
/*
Purpose:
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:
04 July 2013
Author:
C version by John Burkardt
Parameters:
double 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.
double 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.
double 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.
double 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.
double 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.
double FY(M).
FY is the right hand side of the linear system of the refined
problem.
double H(N)
H(I) is the length of subinterval I. This code uses
equal spacing for all the subintervals.
double HY(M).
HY(I) is the length of subinterval I in the refined problem.
int 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.
int 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.
int 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.
int INDY(0:M).
INDY(I) records the index of the unknown associated with
node I for the refined problem.
int JADD(N).
JADD(I) is 1 if the error estimates show that interval I
should be subdivided.
int KOUNT, the number of adaptive steps that have been taken.
int 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.
int N
The number of subintervals into which the interval
[XL,XR] is broken.
int 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.
int NMAX, the maximum number of unknowns that can be handled.
int 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.
int 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.
int NQUAD
The number of quadrature points used in a subinterval.
int 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.
int NUY.
The number of unknowns in the refined problem.
int NY.
NY is the number of subintervals into which a given interval
will be subdivided, before solving the refined probelm.
int 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)
int 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.
double 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.
double 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.
double 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.
double WQUAD(NQUAD).
WQUAD(I) is the weight associated with the I-th point
of an NQUAD point Gaussian quadrature rule.
double XL.
XL is the left endpoint of the interval over which the
differential equation is being solved.
double XN(0:N).
XN(I) is the location of the I-th node. XN(0) is XL,
and XN(N) is XR.
double XQUAD(NQUAD,NMAX), the I-th quadrature point
in interval J.
double XQUADY(NQUAD,NMAY ), the I-th quadrature point
in subinterval J of the refined problem.
double 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.
double YN(0:M).
YN(I) is the location of the I-th node in the refined
problem.
*/
{
# define NL 2
# define NMAX 30
# define NQUAD 2
double alpha;
double adiag[NMAX+1];
double aleft[NMAX+1];
double arite[NMAX+1];
double beta;
double eta[NMAX];
double f[NMAX+1];
double h[NMAX];
int ibc;
int indx[NMAX+1];
int kount;
int n;
int node[NL*NMAX];
int nu;
int problem;
int status;
double tol;
double ul;
double ur;
double wquad[NQUAD];
double xl;
double xn[NMAX+1];
double xquad[NQUAD*NMAX];
double xr;
timestamp ( );
printf ( "\n" );
printf ( "FEM1D_ADAPTIVE\n" );
printf ( " C version\n" );
printf ( "\n" );
printf ( "Solve the two-point boundary value problem:\n" );
printf ( "\n" );
printf ( " -d/dx ( P(x) * dU(x)/dx ) + Q(x) * U(x) = F(x)\n" );
printf ( "\n" );
printf ( "on the interval [0,1], specifying the value\n" );
printf ( "of U at each endpoint.\n" );
printf ( "\n" );
printf ( " The number of basis functions per element is %d\n", NL );
printf ( "\n" );
printf ( " The number of quadrature points per element is %d\n", NQUAD );
problem = get_problem ( );
printf ( "\n" );
printf ( " Problem index = %d\n", problem );
printf ( "\n" );
if ( problem == 1 )
{
printf ( "\n" );
printf ( " \"Linear\" problem:\n" );
printf ( " (No refinement needed)\n" );
printf ( "\n" );
printf ( " U(X) = X\n" );
printf ( " P(X) = 1.0\n" );
printf ( " Q(X) = 0.0\n" );
printf ( " F(X) = 0.0\n" );
printf ( " IBC = 3\n" );
printf ( " UL = 0.0\n" );
printf ( " UR = 1.0\n" );
}
else if ( problem == 2 )
{
printf ( "\n" );
printf ( " \"Quadratic\" problem:\n" );
printf ( " (No refinement needed)\n" );
printf ( "\n" );
printf ( " U(X) = X*X\n" );
printf ( " P(X) = 1.0\n" );
printf ( " Q(X) = 0.0\n" );
printf ( " F(X) = -2.0\n" );
printf ( " IBC = 3\n" );
printf ( " UL = 0.0\n" );
printf ( " UR = 1.\n" );
}
else if ( problem == 3 )
{
printf ( "\n" );
printf ( " \"SINE\" problem:\n" );
printf ( "\n" );
printf ( " U(X) = SIN(PI*X/2)\n" );
printf ( " P(X) = 1.0\n" );
printf ( " Q(X) = 0.0\n" );
printf ( " F(X) = PI*PI*SIN(PI*X/2)/4\n" );
printf ( " IBC = 3\n" );
printf ( " UL = 0.0\n" );
printf ( " UR = 1.0\n" );
}
else if ( problem == 4 )
{
printf ( "\n" );
printf ( " \"COSINE\" problem:\n" );
printf ( "\n" );
printf ( " U(X) = COS(PI*X/2)\n" );
printf ( " P(X) = 1.0\n" );
printf ( " Q(X) = 0.0\n" );
printf ( " F(X) = PI*PI*COS(PI*X/2)/4\n" );
printf ( " IBC = 3\n" );
printf ( " UL = 0.0\n" );
printf ( " UR = 1.0\n" );
}
else if ( problem == 5 )
{
beta = get_beta ( );
printf ( "\n" );
printf ( " \"RHEINBOLDT\" problem:\n" );
printf ( "\n" );
printf ( " U(X) = X**(B+2)/((B+2)*(B+1))\n" );
printf ( " P(X) = 1.0\n" );
printf ( " Q(X) = 1.0\n" );
printf ( " F(X) = -X**B+(X**B+2))/((B+2)*(B+1))\n" );
printf ( " IBC = 3\n" );
printf ( " UL = 0.0\n" );
printf ( " UR = 1/((B+2)*(B+1))\n" );
printf ( " B = %g\n", beta );
}
else if ( problem == 6 )
{
alpha = get_alpha ( );
printf ( "\n" );
printf ( " \"ARCTAN\" problem:\n" );
printf ( "\n" );
printf ( " U(X) = ATAN((X-0.5)/A)\n" );
printf ( " P(X) = 1.0\n" );
printf ( " Q(X) = 0.0\n" );
printf ( " F(X) = 2*A*(X-0.5)/(A**2+(X-0.5)**2)**2\n" );
printf ( " IBC = 3\n" );
printf ( " UL = ATAN(-0.5/A)\n" );
printf ( " UR = ATAN( 0.5/A)\n" );
printf ( " A = %g\n", alpha );
}
/*
Start out with just 4 subintervals.
*/
n = 4;
/*
Initialize values that define the problem.
*/
init ( &ibc, n, &tol, &ul, &ur, &xl, xn, &xr );
/*
Start the iteration counter off at 0.
*/
kount = 0;
/*
Begin the next iteration.
*/
for ( ; ; )
{
kount = kount + 1;
printf ( "\n" );
printf ( " Begin new iteration with %d nodes.\n", n );
printf ( "\n" );
/*
Solve the regular problem.
*/
solvex ( adiag, aleft, arite, f, h, ibc, indx, kount, n, NL, NMAX,
node, NQUAD, &nu, ul, ur, wquad, xn, xquad );
/*
Solve N subproblems to get the error estimators.
*/
solvey ( eta, f, h, n, nu, ul, ur, xn );
/*
Examine the error estimators, and see how many intervals should
be subdivided.
*/
status = subdiv ( eta, kount, &n, NMAX, tol, xn );
if ( status != 0 )
{
break;
}
}
/*
Terminate.
*/
printf ( "\n" );
printf ( "FEM1D_ADAPTIVE:\n" );
printf ( " Normal end of execution.\n" );
printf ( "\n" );
timestamp ( );
return 0;
# undef NL
# undef NMAX
# undef NQUAD
}
/******************************************************************************/
void assemble ( double adiag[], double aleft[], double arite[], double f[],
double h[], int n, int indx[], int node[], int nu, int nl, int nquad,
int nmax, double ul, double ur, double wquad[], double xn[], double xquad[] )
/******************************************************************************/
/*
Purpose:
ASSEMBLE assembles the global matrix.
Licensing:
This code is distributed under the GNU LGPL license.
Modified:
04 July 2013
Author:
C version by John Burkardt
Parameters:
Output, double 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, double 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, double 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, double 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, double H(N)
H(I) is the length of subinterval I. This code uses
equal spacing for all the subintervals.
Input, int N
The number of subintervals into which the interval
[XL,XR] is broken.
Input, int 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, int 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, int 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, int 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, int NQUAD
The number of quadrature points used in a subinterval.
Input, int NMAX, the maximum number of unknowns that can be handled.
Input, double 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, double 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, double WQUAD(NQUAD).
WQUAD(I) is the weight associated with the I-th point
of an NQUAD point Gaussian quadrature rule.
Input, double XN(0:N).
XN(I) is the location of the I-th node. XN(0) is XL,
and XN(N) is XR.
Input, double XQUAD(NQUAD,NMAX), the I-th quadrature point
in interval J.
*/
{
double aij;
double he;
int i;
int ie;
int ig;
int il;
int iq;
int iu;
int jg;
int jl;
int ju;
double phii;
double phiix;
double phij;
double phijx;
double wquade;
double x;
double xleft;
double xquade;
double xrite;
/*
Zero out the entries.
*/
for ( i = 0; i < nu; i++ )
{
f[i] = 0.0;
}
for ( i = 0; i < nu; i++ )
{
aleft[i] = 0.0;
}
for ( i = 0; i < nu; i++ )
{
arite[i] = 0.0;
}
for ( i = 0; i < nu; i++ )
{
adiag[i] = 0.0;
}
/*
For each interval,
*/
for ( ie = 0; ie < n; ie++ )
{
he = h[ie];
xleft = xn[node[0+ie*2]];
xrite = xn[node[1+ie*2]];
/*
For each quadrature point in the interval,
*/
for ( iq = 0; iq < nquad; iq++ )
{
xquade = xquad[iq+ie*nquad];
wquade = wquad[iq];
/*
Pick a basis function which defines the equation,
*/
for ( il = 1; il <= nl; il++ )
{
ig = node[il-1+ie*nl];
iu = indx[ig];
if ( 0 < iu )
{
phi ( il, xquade, &phii, &phiix, xleft, xrite );
f[iu-1] = f[iu-1] + he * wquade * ff ( xquade ) * phii;
/*
Take care of boundary conditions specifying the value of U'.
*/
if ( ig == 0 )
{
x = 0.0;
f[iu-1] = f[iu-1] - pp ( x ) * ul;
}
else if ( ig == n )
{
x = 1.0;
f[iu-1] = f[iu-1] + pp ( x ) * ur;
}
/*
Pick a basis function which defines the coefficient
being computed.
*/
for ( jl = 1; jl <= nl; jl++ )
{
jg = node[jl-1+ie*nl];
ju = indx[jg];
phi ( jl, xquade, &phij, &phijx, 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-1] = f[iu-1] - aij * ul;
}
else if ( jg == n )
{
f[iu-1] = f[iu-1] - aij * ur;
}
}
else if ( iu == ju )
{
adiag[iu-1] = adiag[iu-1] + aij;
}
else if ( ju < iu )
{
aleft[iu-1] = aleft[iu-1] + aij;
}
else
{
arite[iu-1] = arite[iu-1] + aij;
}
}
}
}
}
}
return;
}
/******************************************************************************/
double ff ( double x )
/******************************************************************************/
/*
Purpose:
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:
04 July 2013
Author:
C version by John Burkardt
Parameters:
Input, double X, the evaluation point.
Output, double FF, the value of F(X).
*/
{
double alpha;
double beta;
int problem;
double pi = 3.141592653589793;
double value;
/*
Find out which problem we're working on.
*/
problem = get_problem ( );
if ( problem == 1 )
{
value = 0.0;
}
else if ( problem == 2 )
{
value = -2.0 * x;
}
else if ( problem == 3 )
{
value = 0.25 * pi * pi * sin ( 0.5 * pi * x );
}
else if ( problem == 4 )
{
value = 0.25 * pi * pi * cos ( 0.5 * pi * x );
}
else if ( problem == 5 )
{
beta = get_beta ( );
value = - pow ( x, beta ) + ( pow ( x, beta + 2.0 ) )
/ ( ( beta + 2.0 ) * ( beta + 1.0 ) );
}
else if ( problem == 6 )
{
alpha = get_alpha ( );
value = 2.0 * alpha * ( x - 0.5 )
/ pow ( alpha * alpha + ( x - 0.5 ) * ( x - 0.5 ), 2 );
}
return value;
}
/******************************************************************************/
void geometry ( double h[], int ibc, int indx[], int n, int nl, int nmax,
int node[], int nquad, int *nu, double wquad[], double xn[], double xquad[] )
/******************************************************************************/
/*
Purpose:
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 July 2013
Author:
C version by John Burkardt
Parameters:
Output, double H(N)
H(I) is the length of subinterval I. This code uses
equal spacing for all the subintervals.
Input, int 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, int 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, int N
The number of subintervals into which the interval
[XL,XR] is broken.
Input, int 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, int NMAX, the maximum number of unknowns that can be handled.
Output, int 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, int NQUAD
The number of quadrature points used in a subinterval.
Output, int 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, double WQUAD(NQUAD).
WQUAD(I) is the weight associated with the I-th point
of an NQUAD point Gaussian quadrature rule.
Input, double XN(0:N).
XN(I) is the location of the I-th node. XN(0) is XL,
and XN(N) is XR.
Output, double XQUAD(NQUAD,NMAX), the I-th quadrature point
in interval J.
*/
{
double alfa;
int i;
int igl;
int igr;
double xl;
double xr;
/*
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 = 0; i < n; i++ )
{
node[0+i*2] = i;
node[1+i*2] = i + 1;
}
/*
For every node that is associated with an unknown, we
record the number of the unknown in INDX.
*/
*nu = 0;
for ( i = 0; i <= n; i++ )
{
if ( i == 0 && ( ibc == 1 || ibc == 3 ) )
{
indx[i] = -1;
}
else if ( i == n && ( ibc == 2 || ibc == 3 ) )
{
indx[i] = -1;
}
else
{
*nu = *nu + 1;
indx[i] = *nu;
}
}
/*
We compute the width of each interval.
*/
for ( i = 0; i < n; i++ )
{
igl = node[0+i*2];
igr = node[1+i*2];
h[i] = xn[igr] - xn[igl];
}
/*
We compute the location of the quadrature points in each
interval.
*/
for ( i = 0; i < n; i++ )
{
xl = xn[node[0+i*2]];
xr = xn[node[1+i*2]];
if ( nquad == 1 )
{
xquad[0+i*nquad] = 0.5 * ( xl + xr );
}
else if ( nquad == 2 )
{
alfa = -0.577350;
xquad[0+i*nquad] = ( ( 1.0 - alfa ) * xl
+ ( 1.0 + alfa ) * xr )
/ 2.0;
alfa = +0.577350;
xquad[1+i*nquad] = ( ( 1.0 - alfa ) * xl
+ ( 1.0 + alfa ) * xr )
/ 2.0;
}
else if ( nquad == 3 )
{
alfa = -0.774597;
xquad[0+i*nquad] = ( ( 1.0 - alfa ) * xl
+ ( 1.0 + alfa ) * xr )
/ 2.0;
xquad[1+i*nquad] = 0.5 * ( xl + xr );
alfa = +0.774597;
xquad[2+i*nquad] = ( ( 1.0 - alfa ) * xl
+ ( 1.0 + alfa ) * xr )
/ 2.0;
}
}
/*
Store the weights for the quadrature rule.
*/
if ( nquad == 1 )
{
wquad[0] = 1.0;
}
else if ( nquad == 2 )
{
wquad[0] = 0.5;
wquad[1] = 0.5;
}
else if ( nquad == 3 )
{
wquad[0] = 4.0 / 9.0;
wquad[1] = 5.0 / 18.0;
wquad[2] = 4.0 / 9.0;
}
return;
}
/******************************************************************************/
double get_alpha ( void )
/******************************************************************************/
/*
Purpose:
GET_ALPHA returns the value of ALPHA, for use by problem 6.
Licensing:
This code is distributed under the GNU LGPL license.
Modified:
04 July 2013
Author:
C version by John Burkardt
Parameters:
Output, double GET_ALPHA, the value of ALPHA.
*/
{
double value;
value = 0.01;
return value;
}
/******************************************************************************/
double get_beta ( void )
/******************************************************************************/
/*
Purpose:
GET_BETA returns the value of BETA, for use by problem 5.
Licensing:
This code is distributed under the GNU LGPL license.
Modified:
04 July 2013
Author:
C version by John Burkardt
Parameters:
Output, double VALUE, the value of BETA.
*/
{
double value;
value = -0.9;
return value;
}
/******************************************************************************/
int get_problem ( void )
/******************************************************************************/
/*
Purpose:
GETPRB returns the value of the current problem number.
Licensing:
This code is distributed under the GNU LGPL license.
Modified:
04 July 2013
Author:
C version by John Burkardt
Parameters:
Output, int GET_PROBLEM, the index of the problem.
*/
{
int value;
value = 6;
return value;
}
/******************************************************************************/
void init ( int *ibc, int n, double *tol, double *ul, double *ur, double *xl,
double xn[], double *xr )
/******************************************************************************/
/*
Purpose:
INIT initializes some parameters that define the problem.
Licensing:
This code is distributed under the GNU LGPL license.
Modified:
04 July 2013
Author:
C version by John Burkardt
Parameters:
Output, int *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, int N
The number of subintervals into which the interval
[XL,XR] is broken.
Output, double *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, double *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, double *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, double *XL.
XL is the left endpoint of the interval over which the
differential equation is being solved.
Output, double XN(0:N).
XN(I) is the location of the I-th node. XN(0) is XL,
and XN(N) is XR.
Output, double *XR.
XR is the right endpoint of the interval over which the
differential equation is being solved.
*/
{
double beta;
int i;
int problem;
*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;
printf ( "\n" );
printf ( "Exact solution is U = X\n" );
}
else if ( problem == 2 )
{
*ibc = 3;
*ul = 0.0;
*ur = 1.0;
*xl = 0.0;
*xr = 1.0;
printf ( "\n" );
printf ( "Exact solution is U = X*X\n" );
}
else if ( problem == 3 )
{
*ibc = 3;
*ul = 0.0;
*ur = 1.0;
*xl = 0.0;
*xr = 1.0;
printf ( "\n" );
printf ( "Exact solution is U = SIN(PI*X/2)\n" );
}
else if ( problem == 4 )
{
*ibc = 3;
*ul = 1.0;
*ur = 0.0;
*xl = 0.0;
*xr = 1.0;
printf ( "\n" );
printf ( "Exact solution is U = COS(PI*X/2)\n" );
}
else if ( problem == 5 )
{
*ibc = 3;
beta = get_beta ( );
*ul = 0.0;
*ur = 1.0 / ( ( beta + 2.0 ) * ( beta + 1.0 ) );
*xl = 0.0;
*xr = 1.0;
printf ( "\n" );
printf ( "Rheinboldt problem\n" );
}
else if ( problem == 6 )
{
*ibc = 3;
*xl = 0.0;
*xr = 1.0;
*ul = uexact ( *xl );
*ur = uexact ( *xr );
printf ( "\n" );
printf ( "Arctangent problem\n" );
}
/*
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; i <= n; i++ )
{
xn[i] = ( ( double ) ( n - i ) * ( *xl )
+ ( double ) ( i ) * ( *xr ) )
/ ( double ) ( n );
}
printf ( "The equation is to be solved for \n" );
printf ( "X greater than %g\n", *xl );
printf ( " and less than %g\n", *xr );
printf ( "\n" );
printf ( "The boundary conditions are:\n" );
printf ( "\n" );
if ( *ibc == 1 || *ibc == 3 )
{
printf ( " At X = XL, U = %g\n", *ul );
}
else
{
printf ( " At X = XL, U' = %g\n", *ul );
}
if ( *ibc == 2 || *ibc == 3 )
{
printf ( " At X = XR, U= %g\n", *ur );
}
else
{
printf ( " At X = XR, U' = %g\n", *ur );
}
return;
}
/******************************************************************************/
void output ( double f[], int ibc, int indx[], int n, int nu, double ul,
double ur, double xn[] )
/******************************************************************************/
/*
Purpose:
OUTPUT prints out the computed solution.
Licensing:
This code is distributed under the GNU LGPL license.
Modified:
04 July 2013
Author:
C version by John Burkardt
Parameters:
Input, double 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, int 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, int 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, int N
The number of subintervals into which the interval
[XL,XR] is broken.
Input, int 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, double 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, double 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, double XN(0:N).
XN(I) is the location of the I-th node. XN(0) is XL,
and XN(N) is XR.
*/
{
double error;
int i;
double u;
double uex;
printf ( "\n" );
printf ( "Node X(I) U(X(I)) Uexact Error\n" );
printf ( "\n" );
for ( i = 0; i <= n; i++ )
{
if ( i == 0 )
{
if ( ibc == 1 || ibc == 3 )
{
u = ul;
}
else
{
u = f[indx[i]-1];
}
}
else if ( i == n )
{
if ( ibc == 2 || ibc == 3 )
{
u = ur;
}
else
{
u = f[indx[i]-1];
}
}
else
{
u = f[indx[i]-1];
}
uex = uexact ( xn[i] );
error = u - uex;
printf ( " %4d %13g %12g %12g %12g\n", i, xn[i], u, uex, error );
}
return;
}
/******************************************************************************/
void phi ( int il, double x, double *phii, double *phiix, double xleft,
double xrite )
/******************************************************************************/
/*
Purpose:
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:
04 July 2013
Author:
C version by John Burkardt
Parameters:
Input, int IL, the local index of the basis function.
Input, double X, the evaluation point.
Output, double *PHII, *PHIIX, the value of the basis function
and its derivative.
Input, double XLEFT, XRITE, the endpoints of the interval.
*/
{
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 );
}
}
/*
If X is outside of the interval, then the basis function
is always zero.
*/
else
{
*phii = 0.0;
*phiix = 0.0;
}
return;
}
/******************************************************************************/
double pp ( double x )
/******************************************************************************/
/*
Purpose:
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:
04 July 2013
Author:
C version by John Burkardt
Parameters:
Input, double X, the evaluation point.
Output, double PP, the value of P(X).
*/
{
int problem;
double value;
/*
Find out which problem we're working on.
*/
problem = get_problem ( );
if ( problem == 1 )
{
value = 1.0;
}
else if ( problem == 2 )
{
value = 1.0;
}
else if ( problem == 3 )
{
value = 1.0;
}
else if ( problem == 4 )
{
value = 1.0;
}
else if ( problem == 5 )
{
value = 1.0;
}
else if ( problem == 6 )
{
value = 1.0;
}
return value;
}
/******************************************************************************/
void prsys ( double adiag[], double aleft[], double arite[], double f[],
int nu )
/******************************************************************************/
/*
Purpose:
PRSYS prints out the tridiagonal linear system.
Licensing:
This code is distributed under the GNU LGPL license.
Modified:
04 July 2013
Author:
C version by John Burkardt
Parameters:
Input, double 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, double 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, double 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, double 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.
int 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.
*/
{
int i;
printf ( "\n" );
printf ( "Printout of tridiagonal linear system:\n" );
printf ( "\n" );
printf ( "Equation A-Left A-Diag A-Rite RHS\n" );
printf ( "\n" );
for ( i = 0; i < nu; i++ )
{
printf ( "%4d", i + 1 );
if ( i == 0 )
{
printf ( " " );
}
else
{
printf ( " %12g", aleft[i] );
}
printf ( " %12g", adiag[i] );
if ( i < nu - 1 )
{
printf ( " %12g", arite[i] );
}
else
{
printf ( " " );
}
printf ( " %12g\n", f[i] );
}
return;
}
/******************************************************************************/
double qq ( double x )
/******************************************************************************/
/*
Purpose:
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:
04 July 2013
Author:
C version by John Burkardt
Parameters:
Input, double X, the evaluation point.
Output, double QQ, the value of Q(X).
*/
{
int problem;
double value;
/*
Find out which problem we're working on.
*/
problem = get_problem ( );
if ( problem == 1 )
{
value = 0.0;
}
else if ( problem == 2 )
{
value = 0.0;
}
else if ( problem == 3 )
{
value = 0.0;
}
else if ( problem == 4 )
{
value = 0.0;
}
else if ( problem == 5 )
{
value = 1.0;
}
else if ( problem == 6 )
{
value = 0.0;
}
return value;
}
/******************************************************************************/
void solve ( double adiag[], double aleft[], double arite[], double f[],
int nu )
/******************************************************************************/
/*
Purpose:
SOLVE solves a tridiagonal matrix system of the form A*x = b.
Licensing:
This code is distributed under the GNU LGPL license.
Modified:
04 July 2013
Author:
C version by John Burkardt
Parameters:
Input/output, double ADIAG(NU), ALEFT(NU), ARITE(NU).
On input, ADIAG, ALEFT, and ARITE contain the diagonal,
left and right entries of the equations.
On output, ADIAG and ARITE have been changed in order
to compute the solution.
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/output, double F(NU).
On input, F contains the right hand side of the linear
system to be solve
On output, F contains the solution of the linear system.
Input, int NU, the number of equations to be solved.
*/
{
int i;
/*
Handle the special case of a single equation.
*/
if ( nu == 1 )
{
f[0] = f[0] / adiag[0];
}
/*
The general case, when NU is greater than 1.
*/
else
{
arite[0] = arite[0] / adiag[0];
for ( i = 2; i <= nu - 1; i++ )
{
adiag[i-1] = adiag[i-1] - aleft[i-1] * arite[i-2];
arite[i-1] = arite[i-1] / adiag[i-1];
}
adiag[nu-1] = adiag[nu-1] - aleft[nu-1] * arite[nu-2];
f[0] = f[0] / adiag[0];
for ( i = 2; i <= nu; i++ )
{
f[i-1] = ( f[i-1] - aleft[i-1] * f[i-2] ) / adiag[i-1];
}
for ( i = nu-1; 1 <= i; i-- )
{
f[i-1] = f[i-1] - arite[i-1] * f[i];
}
}
return;
}
/******************************************************************************/
void solvex ( double adiag[], double aleft[], double arite[], double f[],
double h[], int ibc, int indx[], int kount, int n, int nl, int nmax,
int node[], int nquad, int *nu, double ul, double ur, double wquad[],
double xn[], double xquad[] )
/******************************************************************************/
/*
Purpose:
SOLVEX discretizes and solves a differential equation given the nodes.
Licensing:
This code is distributed under the GNU LGPL license.
Modified:
04 July 2013
Author:
C version by John Burkardt
Parameters:
Workspace, double 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.
Workspace, double 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.
Workspace, double 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, double 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.
Output, double H(N)
H(I) is the length of subinterval I. This code uses
equal spacing for all the subintervals.
Input, int 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.
Workspace, int 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, int KOUNT, the number of adaptive steps that have been taken.
Input, int N
The number of subintervals into which the interval
[XL,XR] is broken.
Input, int 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, int NMAX, the maximum number of unknowns that can be handled.
Workspace, int 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.
Workspace, int NQUAD
The number of quadrature points used in a subinterval.
Output, int *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, double 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, double 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.
Workspace, double WQUAD(NQUAD).
WQUAD(I) is the weight associated with the I-th point
of an NQUAD point Gaussian quadrature rule.
Input, double XN(0:N).
XN(I) is the location of the I-th node. XN(0) is XL,
and XN(N) is XR.
Workspace, double XQUAD(NQUAD,NMAX), the I-th quadrature point
in interval J.
*/
{
/*
Given a set of N nodes (where N increases on each iteration),
compute the other geometric information.
*/
geometry ( h, ibc, indx, n, nl, nmax, node, nquad, nu, wquad, xn, xquad );
/*
Assemble the linear system.
*/
assemble ( adiag, aleft, arite, f, 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 );
}
/*
Solve the linear system.
*/
solve ( adiag, aleft, arite, f, *nu );
/*
Print out the solution.
*/
printf ( "\n" );
printf ( "Basic solution\n" );
output ( f, ibc, indx, n, *nu, ul, ur, xn );
return;
}
/******************************************************************************/
void solvey ( double eta[], double f[], double h[], int n, int nu, double ul,
double ur, double xn[] )
/******************************************************************************/
/*
Purpose:
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:
04 July 2013
Author:
C version by John Burkardt
Parameters:
Output, double 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, double 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, double H(N)
H(I) is the length of subinterval I. This code uses
equal spacing for all the subintervals.
Input, int N
The number of subintervals into which the interval
[XL,XR] is broken.
Input, int 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, double 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, double 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, double XN(0:N).
XN(I) is the location of the I-th node. XN(0) is XL,
and XN(N) is XR.
*/
{
# define NL 2
# define NY 2
# define NQUAD 2
# define NMAY ( 2 * NY )
double adiag[NMAY];
double aleft[NMAY];
double arite[NMAY];
double fy[NMAY];
double hy[NMAY];
int i;
int ibcy;
int indy[NMAY+1];
int j;
int jhi;
int jlo;
int jmid;
int k;
int m;
int nodey[NL*NMAY];
int nuy;
double total;
double uleft;
double ulval;
double uly;
double uprime;
double urite;
double urval;
double ury;
double uval;
double vlval;
double vprime;
double vrval;
double vval;
double wquad[NQUAD];
double xquady[NQUAD*NMAY];
double y;
double yl;
double ym;
double yn[NMAY+1];
double yr;
/*
Initialize the error estimators to zero.
*/
for ( j = 0; j < n; j++ )
{
eta[j] = 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; j <= n; j++ )
{
if ( j == 0 )
{
m = NY;
jlo = j;
jmid = j + 1;
jhi = j + 1;
}
else if ( j == n )
{
m = NY;
jlo = j - 1;
jmid = j;
jhi = j;
}
else
{
m = 2 * NY;
jlo = j - 1;
jmid = j;
jhi = j + 1;
}
/*
Set the location of the nodes in the subintervals.
*/
yl = xn[jlo];
ym = xn[jmid];
yr = xn[jhi];
for ( i = 0; i <= NY; i++ )
{
yn[i] = ( ( double ) ( NY - i ) * yl
+ ( double ) ( i ) * ym )
/ ( double ) ( NY );
}
for ( i = NY+1; i <= m; i++ )
{
yn[i] = ( ( double ) ( m - i ) * ym
+ ( double ) ( i - NY ) * yr )
/ ( double ) ( m - NY );
}
/*
Set up the geometry of the sub-problem.
*/
geometry ( hy, ibcy, indy, m, NL, NMAY, nodey, NQUAD, &nuy,
wquad, yn, xquady );
/*
Set the boundary values for the sub-problem.
*/
if ( j <= 1 )
{
uly = ul;
}
else
{
uly = f[j-2];
}
if ( n - 1 <= j )
{
ury = ur;
}
else
{
ury = f[j];
}
/*
Assemble the matrix for the sub-problem.
*/
assemble ( adiag, aleft, arite, fy, hy, m, indy, nodey, nuy, NL,
NQUAD, NMAY, uly, ury, wquad, yn, xquady );
/*
Solve the system.
*/
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 = f[0];
}
else if ( j == n )
{
uleft = f[j-2];
urite = ur;
}
else
{
uleft = f[j-2];
urite = f[j-1];
}
uprime = ( urite - uleft ) / h[j-1];
total = 0.0;
for ( i = 1; i <= NY; i++ )
{
yl = yn[i-1];
yr = yn[i];
if ( i == 1 )
{
vlval = uly;
vrval = fy[i-1];
}
else if ( i == m )
{
vlval = fy[i-2];
vrval = ury;
}
else
{
vlval = fy[i-2];
vrval = fy[i-1];
}
vprime = ( vrval - vlval ) / hy[i-1];
ulval = ( ( double ) ( NY - i + 1 ) * uleft
+ ( double ) ( i - 1 ) * urite )
/ ( double ) ( NY );
urval = ( ( double ) ( NY - i ) * uleft
+ ( double ) ( i ) * urite )
/ ( double ) ( NY );
/*
Compute the integral of
p(x)*(u'(x)-v'(x))**2 + q(x)*(u(x)-v(x))**2
*/
for ( k = 0; k < NQUAD; k++ )
{
y = xquady[k+(i-1)*NQUAD];
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-1] *
( pp ( y ) * pow ( uprime - vprime, 2 )
+ qq ( y ) * pow ( uval - vval, 2 ) );
}
}
eta[j-1] = eta[j-1] + 0.5 * sqrt ( total );
}
/*
Calculation for right interval.
*/
if ( j <= n - 1 )
{
if ( j == 0 )
{
uleft = ul;
urite = f[j];
}
else if ( n - 1 <= j )
{
uleft = f[j-1];
urite = ur;
}
else
{
uleft = f[j-1];
urite = f[j];
}
uprime = ( urite - uleft ) / h[j];
total = 0.0;
for ( i = m+1-NY; i <= m; i++ )
{
yl = yn[i-1];
yr = yn[i];
if ( i == 1 )
{
vlval = uly;
vrval = fy[i-1];
}
else if ( i == m )
{
vlval = fy[i-2];
vrval = ury;
}
else
{
vlval = fy[i-2];
vrval = fy[i-1];
}
vprime = ( vrval - vlval ) / hy[i-1];
ulval = ( ( double ) ( m - i + 1 ) * uleft
+ ( double ) ( NY - m + i - 1 ) * urite )
/ ( double ) ( NY );
urval = ( ( double ) ( m - i ) * uleft
+ ( double ) ( NY - m + i ) * urite )
/ ( double ) ( NY );
/*
Compute the integral of
p(x)*(u'(x)-v'(x))^2 + q(x)*(u(x)-v(x))^2
*/
for ( k = 0; k < NQUAD; k++ )
{
y = xquady[k+(i-1)*NQUAD];
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-1] *
( pp ( y ) * pow ( uprime - vprime, 2 )
+ qq ( y ) * pow ( uval - vval, 2 ) );
}
}
eta[j] = eta[j] + 0.5 * sqrt ( total );
}
}
/*
Print out the error estimators.
*/
printf ( "\n" );
printf ( "ETA\n" );
printf ( "\n" );
for ( j = 0; j < n; j++ )
{
printf ( "%12g\n", eta[j] );
}
return;
# undef NL
# undef NMAY
# undef NQUAD
# undef NY
}
/******************************************************************************/
int subdiv ( double eta[], int kount, int *n, int nmax, double tol,
double xn[] )
/******************************************************************************/
/*
Purpose:
SUBDIV decides which intervals should be subdivided.
Licensing:
This code is distributed under the GNU LGPL license.
Modified:
04 July 2013
Author:
C version by John Burkardt
Parameters:
Input, double 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, int KOUNT, the number of adaptive steps that have been taken.
Input/output, int N
The number of subintervals into which the interval
[XL,XR] is broken.
Input, int NMAX, the maximum number of unknowns that can be handled.
Input, double 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/output, double XN(0:N).
XN(I) is the location of the I-th node. XN(0) is XL,
and XN(N) is XR.
Output, int SUBDIV, 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, int JADD(N).
JADD(I) is 1 if the error estimates show that interval I
should be subdivided.
*/
{
double ave;
int j;
int *jadd;
int k;
int status;
double temp;
double *xt;
status = 0;
/*
Add up the ETA's, and get their average.
*/
ave = 0.0;
for ( j = 0; j < *n; j++ )
{
ave = ave + eta[j];
}
ave = ave / ( double ) ( *n );
/*
Look for intervals whose ETA value is relatively large,
and note in JADD that these intervals should be subdivided.
*/
k = 0;
temp = fmax ( 1.2 * ave + 0.00001, tol * tol / ( double ) ( *n ) );
printf ( "\n" );
printf ( "Tolerance = %g\n", temp );
printf ( "\n" );
jadd = ( int * ) malloc ( *n * sizeof ( int ) );
for ( j = 0; j < *n; j++ )
{
if ( temp < eta[j] )
{
k = k + 1;
jadd[j] = 1;
printf ( "Subdivide interval %d\n", j + 1 );
}
else
{
jadd[j] = 0;
}
}
/*
If no subdivisions needed, we're done.
*/
if ( k <= 0 )
{
printf ( "Success on step %d\n", kount );
status = 1;
return status;
}
/*
See if we're about to go over our limit.
*/
if ( nmax < *n + k )
{
printf ( "\n" );
printf ( "The iterations did not reach their goal.\n" );
printf ( "The next value of N is %d\n", *n + k );
printf ( "which exceeds NMAX = %d\n", nmax );
status = -1;
return status;
}
/*
Insert new nodes where needed.
*/
xt = ( double * ) malloc ( ( nmax + 1 ) * sizeof ( double ) );
k = 0;
xt[0] = xn[0];
for ( j = 0; j < *n; j++ )
{
if ( 0 < jadd[j] )
{
xt[j+1+k] = 0.5 * ( xn[j+1] + xn[j] );
k = k + 1;
}
xt[j+1+k] = xn[j+1];
}
/*
Update the value of N, and copy the new nodes into XN.
*/
*n = *n + k;
for ( j = 0; j <= *n; j++ )
{
xn[j] = xt[j];
}
free ( jadd );
free ( xt );
return status;
}
/******************************************************************************/
void timestamp ( void )
/******************************************************************************/
/*
Purpose:
TIMESTAMP prints the current YMDHMS date as a time stamp.
Example:
31 May 2001 09:45:54 AM
Licensing:
This code is distributed under the GNU LGPL license.
Modified:
24 September 2003
Author:
John Burkardt
Parameters:
None
*/
{
# define TIME_SIZE 40
static char time_buffer[TIME_SIZE];
const struct tm *tm;
time_t now;
now = time ( NULL );
tm = localtime ( &now );
strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm );
printf ( "%s\n", time_buffer );
return;
# undef TIME_SIZE
}
/******************************************************************************/
double uexact ( double x )
/******************************************************************************/
/*
Purpose:
UEXACT returns the value of the exact solution at any point X.
Licensing:
This code is distributed under the GNU LGPL license.
Modified:
04 July 2013
Author:
C version by John Burkardt
Parameters:
Input, double X, the evaluation point.
Output, double UEXACT, the value of the exact solution at X.
*/
{
double alpha;
double beta;
int problem;
double pi = 3.141592653589793;
double value;
/*
Find out which problem we're working on.
*/
problem = get_problem ( );
if ( problem == 1 )
{
value = x;
}
else if ( problem == 2 )
{
value = x * x;
}
else if ( problem == 3 )
{
value = sin ( pi * x / 2.0 );
}
else if ( problem == 4 )
{
value = cos ( pi * x / 2.0 );
}
else if ( problem == 5 )
{
beta = get_beta ( );
value = ( pow ( x, beta + 2.0 ) ) / ( ( beta + 2.0 ) * ( beta + 1.0 ) );
}
else if ( problem == 6 )
{
alpha = get_alpha ( );
value = atan ( ( x - 0.5 ) / alpha );
}
else
{
value = 0.0;
}
return value;
}