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>
#pragma hdrstop

//---------------------------------------------------------------------------
//global variables
int n;      //节点数
double a,b;
double *x;  //节点坐标
double *h;  //单元长度
double *u;  //节点处函数值
double **A; //总刚度矩阵
double *B;
//---------------------------
//问题定义
double p(double y)
{ return 1; }
double q(double y)
{ return 1; }
double f(double y)
{ return 1; }
//-------------------------------------
//开辟内存与释放
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 s,p0,ep,h0,aa,bb,w,y,g;
 double a0=0,b0=1;
 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;
 h0=b0-a0; s=fabs(0.001*h0);
 p0=1.0e+35; ep=eps+1.0;
 while ((ep>=eps)&&(fabs(h0)>s))
 { g=0.0;
  for (i=1;i<=m;i++)
  { aa=a0+(i-1.0)*h0; bb=a0+i*h0;
  w=0.0;
   for (j=0;j<=4;j++)
   { y=((bb-aa)*t[j]+(bb+aa))/2.0;
     w=w+func(k,y)*c[j];
   }
   g=g+w;
  }
  g=g*h0/2.0;
  ep=fabs(g-p0)/(1.0+fabs(g));
  p0=g; m=m+1; h0=(b0-a0)/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)
{
  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)
{}
#pragma argsused
int main(int argc, char* argv[])
{
        return 0;
}
//---------------------------------------------------------------------------