femone1.cpp

本程序采用有限元法(一次区间元)解决一个经典的两点边值问题

//---------------------------------------------------------------------------
/*问题描述
 本程序采用有限元法(一次区间元)解决如下经典的两点边值问题
  -d(p*du/dx)/dx+q*u=f,x(a,b)
  u(a)=0;u(b)=0
  其中:p(x)>=pmin>0,q=q(x)>=0;
*/
//整个程序调试成功
#include<string.h>
#include<math.h>
#include<stdio.h>
#pragma hdrstop

//---------------------------------------------------------------------------
//global variables
int n;      //节点数
double a=0,b=1;
double *x;  //节点坐标
double *h;  //单元长度
double *u;  //节点处函数值
double **A; //总刚度矩阵
double *B;
//---------------------------
//作为测试程序 定义精确解极其一阶 二阶导数
double u0(double x)
{ return x*exp(x)-x*exp(1);}
double u1(double x)
{ return exp(x)*(x+1)-exp(1);}
double u2(double x)
{ return exp(x)*(x+2);}
//---------------------------
//问题定义
double p(double x)
{ return exp(x); }
double p1(double x)
{ return exp(x);}
double q(double x)
{ return x-1; }
double f(double x)
{ return -p1(x)*u1(x)-p(x)*u2(x)+q(x)*u0(x); }
//-------------------------------------
//开辟内存与释放
void allocate(int m)
{
 n=m;
 x=new double[m+1];
 h=new double[m+1];
 u=new double[m+1];
 memset(x,0,(m+1)*sizeof(double));
 memset(h,0,(m+1)*sizeof(double));
 memset(u,0,(m+1)*sizeof(double));
 A=new double*[m+1];
 for(int i=0;i<=m;i++)
  {
   A[i]=new double[m+1];
   memset(A[i],0,(m+1)*sizeof(double));
  }
 B=new double[m+1];
 memset(B,0,(m+1)*sizeof(double));
}
void release(void)
{
 delete[] x;
 delete[] h;
 delete[] u;
 delete[] B;
 for(int i=0;i<=n;i++)
  delete[] A[i];
 delete[] A;
}
//---------------------------
//单元划分(默认为等分)
void divide(int m=10)
{
 allocate(m);
 x[0]=a;
 x[m]=b;
 for(int i=1;i<=m;i++)
 {
  x[i]=x[0]+i*(b-a)/m;
  h[i]=x[i]-x[i-1];
 }
}
//---------------------------
//标准单元基函数
double N0(double y)
{ return 1-y; }
double N1(double y)
{ return y; }
//---------------------------
//几个被积函数
double fa1(int i,double y)
{
 return p(x[i-1]+h[i]*y)/h[i]+
        h[i]*q(x[i-1]+h[i]*y)*N0(y)*N0(y);
}
double fa2(int i,double y)
{
 return -p(x[i-1]+h[i]*y)/h[i]+
        h[i]*q(x[i-1]+h[i]*y)*N0(y)*N1(y);
}
double fa3(int i,double y)
{
 return p(x[i-1]+h[i]*y)/h[i]+
        h[i]*q(x[i-1]+h[i]*y)*N1(y)*N1(y);
}
double fb1(int i,double y)
{
 return f(x[i-1]+h[i]*y)*N0(y)*h[i];
}
double fb2(int i,double y)
{
 return f(x[i-1]+h[i]*y)*N1(y)*h[i];
}
//---------------------------
//高斯积分
double Gauss(double(*func)(int,double),int k,double eps=1e-6)
{
 int m,i,j;
 double a=0,b=1;
 double s,p,ep,h,aa,bb,w,x,g;
 static double t[5]={-0.9061798459,-0.5384693101,0.0,
                         0.5384693101,0.9061798459};
 static double c[5]={0.2369268851,0.4786286705,0.5688888889,
                        0.4786286705,0.2369268851};
 m=1;
 h=b-a; s=fabs(0.001*h);
 p=1.0e+35; ep=eps+1.0;
 while ((ep>=eps)&&(fabs(h)>s))
 { g=0.0;
 for (i=1;i<=m;i++)
 { aa=a+(i-1.0)*h; bb=a+i*h;
 w=0.0;
 for (j=0;j<=4;j++)
 { x=((bb-aa)*t[j]+(bb+aa))/2.0;
   w=w+func(k,x)*c[j];
 }
 g=g+w;
 }
 g=g*h/2.0;
 ep=fabs(g-p)/(1.0+fabs(g));
 p=g; m=m+1; h=(b-a)/m;
 }
 return(g);
}
//---------------------------
//计算有限元方程
void FEE(void)
{
 double I[5];
 for (int i=1;i<=n;i++)
 {
  I[0]=Gauss(fa1,i);
  I[1]=Gauss(fa2,i);
  I[2]=Gauss(fa3,i);
  I[3]=Gauss(fb1,i);
  I[4]=Gauss(fb2,i);
  A[i-1][i-1]+=I[0];
  A[i-1][i]+=I[1];
  A[i][i-1]+=I[1];
  A[i][i]+=I[2];
  B[i-1]+=I[3];
  B[i]+=I[4];
 }
}
//-------------------------
//追赶法解三对角方程组
/*
        |b0  c0         ...||x0  |=d0
        |a1  b1  c1     ...||x1  |=d1
        |... ... ...    ...||... |=...
        |... an-2 bn-2 cn-2||xn-2|=dn-2
        |... ...  an-1 bn-1||xn-1|=dn-1
*/
void TriDiagonal(double *a,double *b,double *c,double *d,int n,double *x)
//注意 a[0] c[n-1]未使用
{
  double *beta,*y;
  beta=new double[n];
  y=new double[n];

  beta[0]=c[0]/b[0];
  for(int i=1;i<=n-2;i++)
   beta[i]=c[i]/(b[i]-a[i]*beta[i-1]);

  y[0]=d[0]/b[0];
  for(int i=1;i<=n-1;i++)
   y[i]=(d[i]-a[i]*y[i-1])/(b[i]-a[i]*beta[i-1]);

  x[n-1]=y[n-1];
  for(int i=n-2;i>=0;i--)
   x[i]=y[i]-beta[i]*x[i+1];

  delete []beta;
  delete []y;
}
//-------------------------
//解有限元方程
void solve(void)
{
 divide(10);
 FEE();
 //调整A B(掐头去尾)
 double *a0,*b0,*c0,*d0,*x0;
 a0=new double[n-1];
 b0=new double[n-1];
 c0=new double[n-1];
 d0=new double[n-1];
 x0=new double[n-1];

 for(int i=0;i<=n-2;i++)
 {
   if(i>0)
   a0[i]=A[i+1][i];
   b0[i]=A[i+1][i+1];
   if(i<n-2)
   c0[i]=A[i+1][i+2];
   d0[i]=B[i+1];
 }
 TriDiagonal(a0,b0,c0,d0,n-1,x0);
 for(int i=0;i<=n-2;i++)
  u[i+1]=x0[i];
 delete[] a0;
 delete[] b0;
 delete[] c0;
 delete[] d0;
 delete[] x0;
}
//-------------------------
//测试追赶法
//Success
void testTri(void)
{
 double a[5]={0,11,8,5,2},
        b[5]={13,10,7,4,1},
        c[5]={12,9,6,3},
        d[5]={3,0,-2,6,8},
        x[5];
 TriDiagonal(a,b,c,d,5,x);
 for(int i=0;i<5;i++)
 printf("x(%d)=%e\n",i+1,x[i]);
}
//-----------------------------
//测试高斯积分
//Success
double g(int i,double x)
{ return exp(x)+i;}
void testGauss(void)
{
 printf("%e",Gauss(g,0));
}
#pragma argsused
int main(int argc, char* argv[])
{
 solve();
 printf("i      x(i)       u(i)(accurate)    ui(calculte)    u(i)-ui\n");
 for(int i=0;i<=n;i++)
 {
  printf("%d    %4f    %e       %e      %e\n",i,x[i],u0(x[i]),u[i],u0(x[i])-u[i]);
 }
 release();
 getchar();
 return 0;
}
//---------------------------------------------------------------------------