onedimension.h

偏微分方程数值解 有限元法 面向对象 变分问题 剖分问题 边值处理 误差分析 椭圆方程

#ifndef Onedimension_h
#define Onedimension_h
#include "FinityElement.h"
#include "mymath.h"
#define filename1 "oneinitial.txt"
#define E 1.0e-8
double h;
double* tout;
class Onedimension:public FinityElement<double>{
private:
 	double height,low,high,steplen,left,right;//区间长度，端点，步长，左右边值
	int N;//剖分点个数
	double *t,*sol,*exactsol,*noncoef;//自变量、解、精确解、非齐次系数
	double stiffness[100][100];//刚度矩阵
    double exactf(double x);//精确解函数
protected:
public:
	Onedimension();
	~Onedimension(){}
	void Equation();//方程
	void Function();//变分问题
	void Dispart();//剖分
	void Radix();//基函数
	void Finity();//有限元方程
	void Calculate();//计算
    void Exact();//精确解;
	void Error();//误差
	friend double F(double x);
    //double f3(double x);
    friend double f4(double x);
    friend double f5(double x); 
};
Onedimension::Onedimension(){
    ifstream fin(filename1,ios::in|ios::nocreate);
	if(!fin){
		cout<<filename1<<" cannot be openned!"<<endl;exit(1);
	}
	cout<<"Input low:";fin>>low;	
	cout<<"Input high:";fin>>high;
	cout<<"Input leftedgevalue:";fin>>left;
	cout<<"Input rightedgevalue:";fin>>right;
	cout<<endl<<"输入剖分点数 N（如5,10,20,40）:";cin>>N;
	fin.close();
}
void Onedimension::Dispart(){
	int i;
	height=high-low;
	h=steplen=height/N;
	t=new double[N+1];
	tout=new double[N+1];
	sol=new double[N+1];
	exactsol=new double[N+1];
	noncoef=new double[N+1];
	for(i=0;i<N+1;i++)
	t[i]=low+i*steplen;
	for(i=0;i<N+1;i++)
	tout[i]=t[i];
}
double F(double x){
	return 2*sin(x*pi/2);
	//return 0;
}
double f1(double x){
	return -1.0/h+h*pi*pi*(1-x)*x/4.0;
}
double f2(double x){
	return 1.0/h+h*pi*pi*x*x/4.0;
}
double f3(double x){
	return 1.0/h+h*pi*pi*(1-x)*(1-x)/4.0;
}
double f4(double x,int j){
	return F(tout[j-1]+h*x)*x;
}
double f5(double x,int j){
	return F(tout[j]+h*x)*(1-x);
}
void Onedimension::Radix(){//初始化刚度矩阵为0
	int i,j;
    for(i=0;i<N+1;i++){
		for(j=0;j<N+1;j++)
			stiffness[i][j]=0;
	}
	for(i=0;i<N+1;i++)
		sol[i]=0;
}
double Onedimension::exactf(double x){
	return (4.0*sin((pi*x)/2))/(pi*pi)+
		((left*exp(-pi/2.0)+2.0*right/pi)/(exp(pi/2)+exp(-pi/2)))*exp(pi*x/2.0)+
		((left*exp(pi/2.0)-2.0*right/pi)/(exp(pi/2)+exp(-pi/2)))*exp(-pi*x/2.0);
	//return 0;
}
void Onedimension::Exact(){
	int i;
	for(i=0;i<N+1;i++)
		exactsol[i]=exactf(t[i]);
}
void Onedimension::Calculate(){
	int j,i;
 	//以下是高斯勒让德积分求刚度矩阵;
	//stiffness[1][1]=GL(f2,0,1,4)+GL(f3,0,1,4);因为有u[0],所以可以统一到下面的循环
	for(j=1;j<N;j++){
		stiffness[j-1][j]=GL(f1,0,1,4);
		stiffness[j][j]=GL(f2,0,1,4)+GL(f3,0,1,4);
		stiffness[j+1][j]=GL(f1,0,1,4);
	}//因为没有u[n+1]所以要单独考虑以下两个元素
	stiffness[N-1][N]=GL(f1,0,1,4);
	stiffness[N][N]=GL(f2,0,1,4);
	//以下用带两个参数的龙贝格积分公式计算有限元方程的非齐次项系数
	for(j=1;j<N;j++)
	{noncoef[j]=h*(LBG2(f4,j,0,1)+LBG2(f5,j,0,1));}
	noncoef[N]=h*LBG2(f4,N,0,1);
	noncoef[1]=noncoef[1]-left*stiffness[0][1];//对左边界作处理
    noncoef[N]=noncoef[N]+right*(1+(right-t[N])/h);//对右边界作处理，这里有可能扩充p(x)
	delete []tout;
    cout<<endl<<"以下输出刚度矩阵及其非齐次项系数："<<endl;
    for(i=1;i<N+1;i++){
		for(j=1;j<N+1;j++)
			cout<<stiffness[i][j]<<" ";cout<<setw(14)<<noncoef[i]<<endl;//
	}
	//用LU分解法解三对角系数方程组理论上精度不高,假如有a[0][0]=0则无法进行
    /*for(i=1;i<N+1;i++){
		for(j=1;j<N+1;j++)
			stiffness[i-1][j-1]=stiffness[i][j];
	}
	for(i=1;i<N+1;i++)
		noncoef[i-1]=noncoef[i]; 
    sol=LU(stiffness,noncoef,N);
    //GS(stiffness,noncoef,sol,N,15);//很难保证高斯勒让德迭代法收敛
	for(i=N;i>0;i--)
	{sol[i]=sol[i-1];
	}
	sol[0]=left;*/
    double *a,*b,*c;
	a=new double[N+1];
	b=new double[N+1];
	c=new double[N+1];
	for(i=2;i<N+1;i++)a[i]=stiffness[i][i-1];
	for(i=1;i<N+1;i++)b[i]=stiffness[i][i];
	for(i=1;i<N;i++)c[i]=stiffness[i][i+1];
    sol[0]=left;
    Chase(sol, a, b, c,noncoef,N);//用追赶法解三对角系数方程组	
	delete []a;
	delete []b;
	delete []c;/**/
}
void Onedimension::Error(){
	double x=0;
	cout<<"以下输出误差分析："<<endl;
    cout<<setw(6)<<"自变量"<<setw(14)<<"有限元解"<<setw(14)<<"精确解"
		<<setw(16)<<"误差"<<endl;
	for(int i=0;i<N+1;i++)
	{
		x=fabs(sol[i]-exactsol[i]);
		cout<<setw(6)<<t[i]<<setw(14)<<sol[i]<<setw(14)<<exactsol[i]<<
			setw(16)<<x<<endl;
	}

}
#endif