C 代码 应用 的 p 方法版本 有限元法 （FEM） 到线性 一维中的两点边界值问题 （BVP）.rar

# include <stdlib.h> # include <stdio.h> # include <math.h> # include <time.h> int main ( ); void alpbet ( double a[], double alpha[], double beta[], int np, int problem, int quad_num, double quad_w[], double quad_x[] ); void exact ( double alpha[], double beta[], double f[], int np, int nprint, int problem, int quad_num, double quad_w[], double quad_x[] ); double ff ( double x, int problem ); void ortho ( double a[], double alpha[], double beta[], int np, int problem, int quad_num, double quad_w[], double quad_x[] ); void out ( double alpha[], double beta[], double f[], int np, int nprint ); void phi ( double alpha[], double beta[], int i, int np, double *phii, double *phiix, double x ); double pp ( double x, int problem ); void quad ( int quad_num, double quad_w[], double quad_x[] ); double qq ( double x, int problem ); void sol ( double a[], double alpha[], double beta[], double f[], int np, int problem, int quad_num, double quad_w[], double quad_x[] ); void timestamp ( void ); double uex ( double x, int problem ); /******************************************************************************/ int main ( void ) /******************************************************************************/ /* Purpose: MAIN is the main program for FEM1D_PMETHOD. Discussion: FEM1D_PMETHOD implements the P-version of the finite element method. Program to solve the one dimensional problem: - d/dX (P dU/dX) + Q U = F by the finite-element method using a sequence of polynomials which satisfy the boundary conditions and are orthogonal with respect to the inner product: (U,V) = Integral (-1 to 1) P U' V' + Q U V dx Here U is an unknown scalar function of X defined on the interval [-1,1], and P, Q and F are given functions of X. The boundary values are U(-1) = U(1)=0. Sample problem #1: U=1-x^4, P=1, Q=1, F=1.0+12.0*x^2-x^4 Sample problem #2: U=cos(0.5*pi*x), P=1, Q=0, F=0.25*pi*pi*cos(0.5*pi*x) The program should be able to get the exact solution for the first problem, using NP = 2. Licensing: This code is distributed under the GNU LGPL license. Modified: 04 July 2013 Author: Original FORTRAN77 version by Max Gunzburger, Teresa Hodge. C version by John Burkardt. Local Parameters: Local, double A[NP+1], the squares of the norms of the basis functions. Local, double ALPHA[NP], BETA[NP], the basis function recurrence coefficients. Local, double F[NP+1]. F contains the basis function coefficients that form the representation of the solution U. That is, U(X) = SUM (I=0 to NP) F(I) * BASIS(I)(X) where "BASIS(I)(X)" means the I-th basis function evaluated at the point X. Local, int NP. The highest degree polynomial to use. Local, int NPRINT. The number of points at which the computed solution should be printed out at the end of the computation. Local, int PROBLEM, indicates the problem being solved. 1, U=1-x^4, P=1, Q=1, F=1.0+12.0*x^2-x^4. 2, U=cos(0.5*pi*x), P=1, Q=0, F=0.25*pi*pi*cos(0.5*pi*x). Local, int QUAD_NUM, the order of the quadrature rule. Local, double QUAD_W[QUAD_NUM], the quadrature weights. Local, double QUAD_X[QUAD_NUM], the quadrature abscissas. */ { # define NP 2 # define QUAD_NUM 10 double a[NP+1]; double alpha[NP]; double beta[NP]; double f[NP+1]; int nprint = 10; int problem = 2; double quad_w[QUAD_NUM]; double quad_x[QUAD_NUM]; timestamp ( ); printf ( "\n" ); printf ( "FEM1D_PMETHOD\n" ); printf ( " C version\n" ); printf ( "\n" ); printf ( " Solve the two-point boundary value problem\n" ); printf ( "\n" ); printf ( " - d/dX (P dU/dX) + Q U = F\n" ); printf ( "\n" ); printf ( " on the interval [-1,1], with\n" ); printf ( " U(-1) = U(1) = 0.\n" ); printf ( "\n" ); printf ( " The P method is used, which represents U as\n" ); printf ( " a weighted sum of orthogonal polynomials.\n" ); printf ( "\n" ); printf ( " Highest degree polynomial to use is %d\n", NP ); printf ( " Number of points to be used for output = %d\n", nprint ); if ( problem == 1 ) { printf ( "\n" ); printf ( " Problem #1:\n" ); printf ( " U=1-x^4,\n" ); printf ( " P=1,\n" ); printf ( " Q=1,\n" ); printf ( " F=1 + 12 * x^2 - x^4\n" ); } else if ( problem == 2 ) { printf ( "\n" ); printf ( " Problem #2:\n" ); printf ( " U=cos(0.5*pi*x),\n" ); printf ( " P=1,\n" ); printf ( " Q=0,\n" ); printf ( " F=0.25*pi*pi*cos(0.5*pi*x)\n" ); } /* Get quadrature abscissas and weights for interval [-1,1]. */ quad ( QUAD_NUM, quad_w, quad_x ); /* Compute the constants for the recurrence relationship that defines the basis functions. */ alpbet ( a, alpha, beta, NP, problem, QUAD_NUM, quad_w, quad_x ); /* Test the orthogonality of the basis functions. */ ortho ( a, alpha, beta, NP, problem, QUAD_NUM, quad_w, quad_x ); /* Solve for the solution of the problem, in terms of coefficients of the basis functions. */ sol ( a, alpha, beta, f, NP, problem, QUAD_NUM, quad_w, quad_x ); /* Print out the solution, evaluated at each of the NPRINT points. */ out ( alpha, beta, f, NP, nprint ); /* Compare the computed and exact solutions. */ exact ( alpha, beta, f, NP, nprint, problem, QUAD_NUM, quad_w, quad_x ); /* Terminate. */ printf ( "\n" ); printf ( "FEM1D_PMETHOD\n" ); printf ( " Normal end of execution.\n" ); printf ( "\n" ); timestamp ( ); return 0; # undef NP # undef QUAD_NUM } /******************************************************************************/ void alpbet ( double a[], double alpha[], double beta[], int np, int problem, int quad_num, double quad_w[], double quad_x[] ) /******************************************************************************/ /* Purpose: ALPBET calculates the coefficients in the recurrence relationship. Discussion: ALPHA and BETA are the coefficients in the three term recurrence relation for the orthogonal basis functions on [-1,1]. The routine also calculates A, the square of the norm of each basis function. Licensing: This code is distributed under the GNU LGPL license. Modified: 04 July 2013 Author: Original FORTRAN77 version by Max Gunzburger, Teresa Hodge. C version by John Burkardt. Parameters: Output, double A[NP+1], the squares of the norms of the basis functions. Output, double ALPHA[NP], BETA[NP], the basis function recurrence coefficients. Input, int NP. The highest degree polynomial to use. Input, int PROBLEM, indicates the problem being solved. 1, U=1-x^4, P=1, Q=1, F=1.0+12.0*x^2-x^4. 2, U=cos(0.5*pi*x), P=1, Q=0, F=0.25*pi*pi*cos(0.5*pi*x). Input, int QUAD_NUM, the order of the quadrature rule. Input, double QUAD_W(QUAD_NUM), the quadrature weights. Input, double QUAD_X(QUAD_NUM), the quadrature abscissas. */ { int i; int iq; int k; double q; double qm1; double qm1x; double qm2; double qm2x; double qx; double s; double ss; double su; double sv; double t; double u; double v; double x; ss = 0.0; su = 0.0; for ( iq = 0; iq < quad_num; iq++ ) { x = quad_x[iq]; s = 4.0 * pp ( x, problem ) * x * x + qq ( x, problem ) * ( 1.0 - x * x ) * ( 1.0 - x * x ); u = 2.0 * pp ( x, problem ) * x * ( 3.0 * x * x - 1.0 ) + x * qq ( x, problem ) * ( 1.0 - x * x ) * ( 1.0 - x * x ); ss = ss + s * quad_w[iq]; su = su + u * quad_w[iq]; } a[0] = ss; alpha[0] = su / ss; beta[0] = 0.0; for ( i = 1; i <= np; i++ ) { ss = 0.0; su = 0.0; sv = 0.0; for ( iq = 0; iq < quad_num; iq++ ) { x = quad_x[iq]; /* Three term recurrence for Q. */ qm1 = 0.0; q = 1.0; for ( k = 0; k <= i-1; k++ ) { qm2 = qm1;