C 代码 应用有限元法， 对于分段线性元素， 到一个空间维度上的两点边界值问题 （BVP）.rar

# 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