# include <math.h>
# include <stdio.h>
# include <stdlib.h>
# include <time.h>
# include "fem1d_bvp_linear.h"
/******************************************************************************/
double *fem1d_bvp_linear ( int n, double a ( double x ), double c ( double x ),
double f ( double x ), double x[] )
/******************************************************************************/
/*
Purpose:
FEM1D_BVP_LINEAR solves a two point boundary value problem.
Discussion:
The program uses the finite element method, with piecewise linear basis
functions to solve a boundary value problem in one dimension.
The problem is defined on the region 0 <= x <= 1.
The following differential equation is imposed between 0 and 1:
- d/dx a(x) du/dx + c(x) * u(x) = f(x)
where a(x), c(x), and f(x) are given functions.
At the boundaries, the following conditions are applied:
u(0.0) = 0.0
u(1.0) = 0.0
A set of N equally spaced nodes is defined on this
interval, with 0 = X(1) < X(2) < ... < X(N) = 1.0.
At each node I, we associate a piecewise linear basis function V(I,X),
which is 0 at all nodes except node I. This implies that V(I,X) is
everywhere 0 except that
for X(I-1) <= X <= X(I):
V(I,X) = ( X - X(I-1) ) / ( X(I) - X(I-1) )
for X(I) <= X <= X(I+1):
V(I,X) = ( X(I+1) - X ) / ( X(I+1) - X(I) )
We now assume that the solution U(X) can be written as a linear
sum of these basis functions:
U(X) = sum ( 1 <= J <= N ) U(J) * V(J,X)
where U(X) on the left is the function of X, but on the right,
is meant to indicate the coefficients of the basis functions.
To determine the coefficient U(J), we multiply the original
differential equation by the basis function V(J,X), and use
integration by parts, to arrive at the I-th finite element equation:
Integral A(X) * U'(X) * V'(I,X) + C(X) * U(X) * V(I,X) dx
= Integral F(X) * V(I,X) dx
We note that the functions U(X) and U'(X) can be replaced by
the finite element form involving the linear sum of basis functions,
but we also note that the resulting integrand will only be nonzero
for terms where J = I - 1, I, or I + 1.
By writing this equation for basis functions I = 2 through N - 1,
and using the boundary conditions, we have N linear equations
for the N unknown coefficients U(1) through U(N), which can
be easily solved.
Licensing:
This code is distributed under the GNU LGPL license.
Modified:
18 June 2014
Author:
John Burkardt
Parameters:
Input, int N, the number of nodes.
Input, double A ( double X ), evaluates a(x);
Input, double C ( double X ), evaluates c(x);
Input, double F ( double X ), evaluates f(x);
Input, double X[N], the mesh points.
Output, double FEM1D_BVP_LINEAR[N], the finite element coefficients,
which are also the value of the computed solution at the mesh points.
*/
{
# define QUAD_NUM 2
double abscissa[QUAD_NUM] = {
-0.577350269189625764509148780502,
+0.577350269189625764509148780502 };
double *amat;
double axq;
double *b;
double cxq;
int e;
int e_num;
double fxq;
int i;
int ierror;
int j;
int l;
int q;
int quad_num = QUAD_NUM;
int r;
double *u;
double weight[QUAD_NUM] = { 1.0, 1.0 };
double wq;
double vl;
double vlp;
double vr;
double vrp;
double xl;
double xq;
double xr;
/*
Zero out the matrix and right hand side.
*/
amat = r8mat_zero_new ( n, n );
b = r8vec_zero_new ( n );
e_num = n - 1;
for ( e = 0; e < e_num; e++ )
{
l = e;
r = e + 1;
xl = x[l];
xr = x[r];
for ( q = 0; q < quad_num; q++ )
{
xq = ( ( 1.0 - abscissa[q] ) * xl
+ ( 1.0 + abscissa[q] ) * xr )
/ 2.0;
wq = weight[q] * ( xr - xl ) / 2.0;
vl = ( xr - xq ) / ( xr - xl );
vlp = - 1.0 / ( xr - xl );
vr = ( xq - xl ) / ( xr - xl );
vrp = + 1.0 / ( xr - xl );
axq = a ( xq );
cxq = c ( xq );
fxq = f ( xq );
amat[l+l*n] = amat[l+l*n] + wq * ( vlp * axq * vlp + vl * cxq * vl );
amat[l+r*n] = amat[l+r*n] + wq * ( vlp * axq * vrp + vl * cxq * vr );
b[l] = b[l] + wq * ( vl * fxq );
amat[r+l*n] = amat[r+l*n] + wq * ( vrp * axq * vlp + vr * cxq * vl );
amat[r+r*n] = amat[r+r*n] + wq * ( vrp * axq * vrp + vr * cxq * vr );
b[r] = b[r] + wq * ( vr * fxq );
}
}
/*
Equation 1 is the left boundary condition, U(0.0) = 0.0;
*/
for ( j = 0; j < n; j++ )
{
amat[0+j*n] = 0.0;
}
b[0] = 0.0;
for ( i = 1; i < n; i++ )
{
b[i] = b[i] - amat[i+0*n] * b[0];
}
for ( i = 0; i < n; i++ )
{
amat[i+0*n] = 0.0;
}
amat[0+0*n] = 1.0;
/*
Equation N is the right boundary condition, U(1.0) = 0.0;
*/
for ( j = 0; j < n; j++ )
{
amat[n-1+j*n] = 0.0;
}
b[n-1] = 0.0;
for ( i = 0; i < n - 1; i++ )
{
b[i] = b[i] - amat[i+(n-1)*n] * b[n-1];
}
for ( i = 0; i < n; i++ )
{
amat[i+(n-1)*n] = 0.0;
}
amat[n-1+(n-1)*n] = 1.0;
/*
Solve the linear system.
*/
u = r8mat_solve2 ( n, amat, b, &ierror );
free ( amat );
free ( b );
return u;
# undef QUAD_NUM
}
/******************************************************************************/
double h1s_error_linear ( int n, double x[], double u[],
double exact_ux ( double x ) )
/******************************************************************************/
/*
Purpose:
H1S_ERROR_LINEAR estimates the seminorm error of a finite element solution.
Discussion:
We assume the finite element method has been used, over an interval [A,B]
involving N nodes, with piecewise linear elements used for the basis.
The coefficients U(1:N) have been computed, and a formula for the
exact derivative is known.
This function estimates the seminorm of the error:
SEMINORM = Integral ( A <= X <= B ) ( dU(X)/dx - EXACT_UX(X) )^2 dX
Licensing:
This code is distributed under the GNU LGPL license.
Modified:
18 June 2014
Author:
John Burkardt
Parameters:
Input, int N, the number of nodes.
Input, double X(N), the mesh points.
Input, double U(N), the finite element coefficients.
Input, function EQ = EXACT_UX ( X ), returns the value of the exact
derivative at the point X.
Output, double H1S_ERROR_LINEAR, the estimated seminorm of
the error.
*/
{
# define QUAD_NUM 2
double abscissa[QUAD_NUM] = {
-0.577350269189625764509148780502,
+0.577350269189625764509148780502 };
double exq;
double h1s;
int i;
int q;
int quad_num = QUAD_NUM;
double ul;
double ur;
double uxq;
double weight[QUAD_NUM] = { 1.0, 1.0 };
double wq;
double xl;
double xq;
double xr;
h1s = 0.0;
/*
Integrate over each interval.
*/
for ( i = 0; i < n - 1; i++ )
{
xl = x[i];
xr = x[i+1];
ul = u[i];
ur = u[i+1];
for ( q = 0; q < quad_num; q++ )
{
xq = ( ( 1.0 - abscissa[q] ) * xl
+ ( 1.0 + abscissa[q] ) * xr )
/ 2.0;
wq = weight[q] * ( xr - xl ) / 2.0;
/*
The piecewise linear derivative is a constant in the interval.
*/
uxq = ( ur - ul ) / ( xr - xl );
exq = exact_ux ( xq );
h1s = h1s + wq * pow ( uxq - exq, 2);
}
}
h1s = sqrt ( h1s );
return h1s;
# undef QUAD_NUM
}
/******************************************************************************/
int i4_power ( int i, int j )
/******************************************************************************/
/*
Purpose:
I4_POWER returns the value of I^J.
Licensing:
This code is distributed under the GNU LGPL license.
Modified:
23 October 2007
Author:
John Burkardt
Parameters:
Input, int I, J, the base and the power. J should be nonnegative.
Output, int I4_POWER, the value of I^J.
*/
{
int k;
int valu