delaunay.txt

二维有限元方法的计算机实现

二维有限元方法
实验目的： 
           －(uxx+uyy)+u=f(x,y)
           u(x,0)=0,        0≤x≤1
           u(0,y)=0,      0≤y≤1
           u(x,1-x)=0,      0≤x≤1
的真解    u*=sinπxsinπysinπ(1-x-y)
1.	求出f(x,y).
2.	用平行于x轴，y轴，直线y＝1－x的三种直线对Ω作三角剖分，用线性元求解元问题的uh（x，y）.
实验原理：
先求得f（x，y）＝（1＋4π2）sinπxsinπysinπ(1-x-y)－π2sin2(πx+πy).
A(u,v)=∮(uxvx+uyvy+uv)dxdy
b(v)= ∮f(x)vdx.
A(u,v)=∑∮ei(uxvx+uyvy+uv)dxdy=uT∑∮ei(NTxNx+NTyNy+NTN)dxdyV.
b(v)= ∑∮ei f(x)vdx=∑∮ei f(x)NTdxV.
其中
      N =(Ni，Nj，Nm）.
      Nx=(Nix,Njx,Nmx). 
Ny=(Niy,Njy,Nmy).
误差：
	N=6	N=10	N=15
ε1	0.105637	0.0626752	0.0402863
ε2	0.218895	0.190746	0.1767
ε3	0.168505	0.0881016	0.0544093
源代码：
#include "math.h"
#define pi 3.141593
void Get_Ab(double **A,double *b,int n);
void GuassSeidel(double **A,double *b,double *x,int n);//用高斯-赛德尔迭代法解线性方程组;
double abs_double(double x);
void Get_AA_bb(double **A,double *b,int n,double **&AA,double *&bb,int &LofXiaBiao,double *&uu,int *&XiaBiao);
double f(double x,double y);
void main()
{
	int i,j,n,dim,*XiaBiao,LofXiaBiao=0;
	double **A,**AA,*b,*bb,*u,*uu,e1,e2,e3;
	cout<<"请输入[0,1]区间段数N：\n";
	cin>>n;
	dim=(n+1)*(n+2)/2;
	A=new double*[dim];
	b=new double[dim];
	u=new double[dim];
	for(i=0;i<dim;i++)
	{
		A[i]=new double[dim];
		for(j=0;j<dim;j++)
			A[i][j]=0;
		b[i]=0;
		u[i]=0;
	}
	Get_Ab(A,b,n);
	Get_AA_bb(A,b,n,AA,bb,LofXiaBiao,uu,XiaBiao);
	GuassSeidel(AA,bb,uu,LofXiaBiao);
	for(i=0;i<LofXiaBiao;i++)
		u[XiaBiao[i]]=uu[i];
	for(i=0;i<=n;i++)
	{
		for(j=0;j<=i;j++)
			cout<<u[i*(i+1)/2+j]<<"       ";
		cout<<endl;
	}
	//误差估计
	e1=abs_double(u[0]-f(0,1));
	for(i=1;i<=n;i++)
	{
		for(j=0;j<=i;j++)
		if(e1<abs_double(u[i*(i+1)/2+j]-f(j/double(n),1-i/double(n))))
			e1=abs_double(u[i*(i+1)/2+j]-f(j/double(n),1-i/double(n)));
	}
	cout<<"e1="<<e1<<endl;
	e1=e2=0;
    for(i=0;i<=n;i++)
		for(j=0;j<=i;j++)
	{
		e1+=f(j/double(n),1-i/double(n))*f(j/double(n),1-i/double(n));
		e2+=(u[i*(i+1)/2+j]-f(j/double(n),1-i/double(n)))*(u[i*(i+1)/2+j]-f(j/double(n),1-i/double(n)));
	}
	e1=sqrt(e1);
	e2=sqrt(e2);
	e3=e2/e1;
	cout<<"e2="<<e2<<endl;
	cout<<"e3="<<e3<<endl;
}
void Get_Ab(double **A,double *b,int n)//求得矩阵A和向量b
{
	int i,j,k,l,m;
	double x,y;
	for(k=0;k<n;k++)
		for(l=0;l<=k;l++)
		{
			x=(3*l+1)/double(3*n);
			y=1-(3*k+2)/double(3*n);
			i=k*(k+1)/2+l;
			j=(k+1)*(k+2)/2+l;
			m=(k+1)*(k+2)/2+l+1;
			A[i][i]=A[i][i]+0.5+1/double(12*n*n);
			A[j][i]=A[i][j]=A[i][j]-0.5+1/double(24*n*n);
			A[j][j]=A[j][j]+1+1/double(12*n*n);
			A[m][j]=A[j][m]=A[j][m]-0.5+1/double(24*n*n);
			A[m][m]=A[m][m]+0.5+1/double(12*n*n);
			A[i][m]=A[m][i]=A[m][i]+1/double(24*n*n);
			b[i]+=((1+4*pi*pi)*sin(pi*x)*sin(pi*y)*sin(pi*x+pi*y)-pi*pi*sin(2*pi*x+2*pi*y))/double(6*n*n);
			b[j]+=((1+4*pi*pi)*sin(pi*x)*sin(pi*y)*sin(pi*x+pi*y)-pi*pi*sin(2*pi*x+2*pi*y))/double(6*n*n);
			b[m]+=((1+4*pi*pi)*sin(pi*x)*sin(pi*y)*sin(pi*x+pi*y)-pi*pi*sin(2*pi*x+2*pi*y))/double(6*n*n);
			if(l<k)
			{
				i=(k+1)*(k+2)/2+l+1;
			    j=(k+1)*k/2+l+1;
			    m=(k+1)*k/2+l;
			    A[i][i]=A[i][i]+0.5+1/double(12*n*n);
			    A[j][i]=A[i][j]=A[i][j]-0.5+1/double(24*n*n);
			    A[j][j]=A[j][j]+1+1/double(12*n*n);
			    A[m][j]=A[j][m]=A[j][m]-0.5+1/double(24*n*n);
			    A[m][m]=A[m][m]+0.5+1/double(12*n*n);
			    A[i][m]=A[m][i]=A[m][i]+1/double(24*n*n);
			    b[i]+=((1+4*pi*pi)*sin(pi*x)*sin(pi*y)*sin(pi*x+pi*y)-pi*pi*sin(2*pi*x+2*pi*y))/double(6*n*n);
			    b[j]+=((1+4*pi*pi)*sin(pi*x)*sin(pi*y)*sin(pi*x+pi*y)-pi*pi*sin(2*pi*x+2*pi*y))/double(6*n*n);
			    b[m]+=((1+4*pi*pi)*sin(pi*x)*sin(pi*y)*sin(pi*x+pi*y)-pi*pi*sin(2*pi*x+2*pi*y))/double(6*n*n);
			}
		}
}
//求得矩阵AA和向量bb，u中非0点的方程组
void Get_AA_bb(double **A,double *b,int n,double **&AA,double *&bb,int &LofXiaBiao,double *&uu,int *&XiaBiao)
{
	int i=0,j,k,l;
	for(k=1;k<n;k++)
		for(l=1;l<k;l++)
			LofXiaBiao++;
	XiaBiao=new int [LofXiaBiao];
	AA=new double*[LofXiaBiao];
	for(j=0;j<LofXiaBiao;j++)
		AA[j]=new double[LofXiaBiao];
	bb=new double[LofXiaBiao];
	uu=new double[LofXiaBiao];
	for(k=1;k<n;k++)
		for(l=1;l<k;l++)
			XiaBiao[i++]=k*(k+1)/2+l;
	for(i=0;i<LofXiaBiao;i++)
	{
		for(j=0;j<LofXiaBiao;j++)
			AA[i][j]=A[XiaBiao[i]][XiaBiao[j]];
		bb[i]=b[XiaBiao[i]];
	}
}
void GuassSeidel(double **A,double *b,double *x,int n)//用高斯-赛德尔迭代法解线性方程组
{
	double p0,p,s;
	int i,j;
	while(1)
	{
		p0=0;
		for(i=0;i<n;i++)
		{
			s=0;
			for(j=0;j<n;j++)
				s+=A[i][j]*x[j];
			p=(b[i]-s)/A[i][i];
			if(abs_double(p)>abs_double(p0))
				p0=p;
			x[i]+=p;
			
		}
		if(abs_double(p0)<0.000001)
			break;
	}
}

double abs_double(double x)
{
	if(x<0)
		return -x;
	else
		return x;
}
double f(double x,double y)//精确解u*
{
	return sin(pi*x)*sin(pi*y)*sin(pi-pi*x-pi*y);
}