comput.cpp

自己写的高精度紧致差分格式

/*////////////////////////////////////////////////////////////////

       各种格式的比较!!!!! 正确!!!

////////////////////////////////////////////////////////////////*/

#include <iostream>
#include <cmath>
#include <fstream>
#include <iomanip>

using namespace std;

class ComputScheme 
{
public:
    ComputScheme();
    void init();

    void PadeScheme();       //Pade格式(FST)
    void SCScheme();         //二阶精度中心格式(SLW)
    void SUWScheme();        //二阶精度迎风格式(MXD)
    void WSCScheme();        //宽网格基二阶精度中心格式(SLW)
    void FUWCScheme();       //五阶精度迎风紧致格式(MXD)
	
    void Output(char* filename);
private:
	int     n;
	int     nt;
	double  dx;
	double* x;
	double* u;
	double* y;

	double  *a;
	double  *b;
	double  *c;
	double  *F;

       double  *v1;
	double  *v2;
	double  *v3;
	double  *v4;
	double  *v5;

       double  *k1;
	double  *k2;
	double  *k3;
	double  *k4;
	double  *k5;
	double  *k6;
	
	double  tinit;
	double  tend;
	double  dt;
	double  cfl;
	
	void    tridiag(double *a, double *b, double *c, double *F, double *y, int n);
};

ComputScheme::ComputScheme()
{
	int i;
	int a0=100;
//	int a0=157;
//	int a0=194;
	int a1=25;
	
	n=1600;
	dx=8.0/n;
	x=new double [n+2];
	y=new double [n+2];
	u=new double [n+2];
	for(i=0;i<=n+1;i++){
		x[i]=i*dx-4;
		y[i]=x[i];
	}
	for(i=0;i<=n+1;i++)
		u[i]=exp(-16*x[i]*x[i])* sin(a0*x[i]) + ((exp(-16*(x[i]-1.5)*(x[i]-1.5))) +(exp(-16*(x[i]+1.5)*(x[i]+1.5))))*sin(a1*x[i]);
	Output("init.plt");
}

void ComputScheme::tridiag(double *a, double *b, double *c, double *F, double *y, int n)
{
	int i;
	double *g=new double [n+1];
	double *w=new double [n+1];

	g[1]=F[1]/b[1];
	w[1]=c[1]/b[1];

	for (i=2; i<n-1; i++)
	{
		g[i]=(F[i]+a[i]*g[i-1])/(b[i]-a[i]*w[i-1]);
		w[i]=c[i]/(b[i]-a[i]*w[i-1]);
	}
	y[n-1]=(F[n-1]+a[n-1]*g[n-2])/(b[n-1]-a[n-1]*w[n-2]);
	for (i=n-2; i>0; i--)
		y[i]=g[i]+w[i]*y[i+1];
}

void ComputScheme::init()
{
	tinit=0;
	tend=1.0;
	dt=0.1*dx;
	
	a=new double [n+2];
	b=new double [n+2];
	c=new double [n+2];
	F=new double [n+2];

	v1=new double [n+2];
	v2=new double [n+2];
	v3=new double [n+2];
	v4=new double [n+2];
	v5=new double [n+2];

	k1=new double [n+2];
	k2=new double [n+2];
	k3=new double [n+2];
	k4=new double [n+2];
	k5=new double [n+2];
	k6=new double [n+2];

	for (int i=0; i<=n+1; i++)
	{
		a[i]=b[i]=c[i]=F[i]=0;
		v1[i]=v2[i]=v3[i]=v4[i]=v5[i]=0;
		k1[i]=k2[i]=k3[i]=k4[i]=k5[i]=k6[i]=0;
	}	
}



void ComputScheme::SCScheme()
{
	double t=0;
	int i;
	int j=0;
 

    while (t<tend)
    {
		
	    for (i=1; i<n+1; i++)
	        y[i]=(u[i+1]-u[i-1])/2;
		for (i=1; i<n+1; i++)
			v1[i]=u[i]-dt*y[i]/dx/2;

		for (i=1; i<n+1; i++)
	        y[i]=(v1[i+1]-v1[i-1])/2;
		for (i=1; i<n+1; i++)
			u[i]=u[i]-dt*y[i]/dx;

        u[0]=0;
		u[n+1]=0;
		t+=dt;

		cout<<setw(8)<<t<<endl;		
    }
	
}

void ComputScheme::SUWScheme()
{
	double t=0;
	int i;

    while (t<tend)
    {
		
	    for (i=2; i<=n+1; i++)
	        y[i]=(3*u[i]-4*u[i-1]+u[i-2])/2;
		for (i=2; i<=n+1; i++)
			v1[i]=u[i]-dt*y[i]/dx/2;

		for (i=2; i<=n+1; i++)
	        y[i]=(3*v1[i]-4*v1[i-1]+v1[i-2])/2;
		for (i=2; i<=n+1; i++)
			u[i]=u[i]-dt*y[i]/dx;

        u[0]=0;
//		u[1]=0;
		u[n+1]=0;
		t+=dt;
		cout<<setw(8)<<t<<endl;
    }	
}

void ComputScheme::WSCScheme()
{
    double t=0;
	int i;
    while (t<tend)
    {
	    for (i=2; i<n; i++)
	        y[i]=(u[i+2]-u[i-2])/4;
		y[1]=y[n]=0;
		for (i=1; i<n+1; i++)
			v1[i]=u[i]-dt*y[i]/dx/2;
		

		for (i=2; i<n; i++)
	        y[i]=(v1[i+2]-v1[i-2])/4;
		y[1]=y[n]=0;
		for (i=1; i<n+1; i++)
			u[i]=u[i]-dt*y[i]/dx;

        u[0]=0;
//		u[1]=0;
//		u[n]=0;
		u[n+1]=0;
		t+=dt;
		cout<<setw(8)<<t<<endl;
    }	
}

void ComputScheme::PadeScheme()
{
	double t=0;
	int i;
	int j=0;
	
	for (i=1; i<n+1; i++)
	{
		a[i]=-1.0/2;
		b[i]=1.0/2;
		c[i]=0;
//		F[i]=u[i]-u[i-1];
	} 
	a[0]=0;
	c[n]=0;

    while (t<tend)
    {
		
		for (i=1; i<=n; i++)
	        F[i]=u[i]-u[i-1];
		y[0]=0;
		y[n+1]=0;

		tridiag(a, b, c, F, y, n+1);
		for (i=1; i<=n; i++)
			v1[i]=u[i]-dt*y[i]/dx/2;
		v1[0]=0.0;
		v1[n]=0.0;

		for (i=1; i<=n; i++)
	        F[i]=v1[i]-v1[i-1];
		tridiag(a, b, c, F, y, n+1);
		for (i=1; i<=n; i++)
			u[i]=u[i]-dt*y[i]/dx;

        u[0]=0;
		u[n]=0;
		t+=dt;
		cout<<setw(8)<<t<<endl;		
    }
}

void ComputScheme::FUWCScheme()    
{
	 double t=0;
	 int i;
	
	for (i=1; i<n+1; i++)
	{
		a[i]=-2.0/5;
		b[i]=3.0/5;
		c[i]=0;
	}
	a[0]=0;
	c[n]=0;

	while (t<tend)
	{
		for (i=2; i<n; i++)
			F[i]=(-u[i+2]+12*u[i+1]+36*u[i]-44*u[i-1]-3*u[i-2])/60;
		F[1]=F[n]=0;
		tridiag(a, b, c, F, y, n+1);
		for (i=1; i<n+1; i++)
			k1[i]=-y[i]/dx;

		for (i=1; i<n+1; i++)
			v1[i]=u[i]+dt*k1[i]/2;
		for (i=2; i<n; i++)
			F[i]=(-v1[i+2]+12*v1[i+1]+36*v1[i]-44*v1[i-1]-3*v1[i-2])/60;
		F[1]=F[n]=0;
		tridiag(a, b, c, F, y, n+1);
		for (i=1; i<n+1; i++)
			k2[i]=-y[i]/dx;

		for (i=1; i<n+1; i++)
			v2[i]=u[i]+dt*(3*k1[i]+k2[i])/16;
		for (i=2; i<n; i++)
			F[i]=(-v2[i+2]+12*v2[i+1]+36*v2[i]-44*v2[i-1]-3*v2[i-2])/60;
		F[1]=F[n]=0;
		tridiag(a, b, c, F, y, n+1);
		for (i=1; i<n+1; i++)
			k3[i]=-y[i]/dx;

		for (i=1; i<n+1; i++)
			v3[i]=u[i]+dt*k3[i]/2;
		for (i=2; i<n; i++)
			F[i]=(-v3[i+2]+12*v3[i+1]+36*v3[i]-44*v3[i-1]-3*v3[i-2])/60;
		F[1]=F[n]=0;
		tridiag(a, b, c, F, y, n+1);
		for (i=1; i<n+1; i++)
			k4[i]=-y[i]/dx;

		for (i=1; i<n+1; i++)
			v4[i]=u[i]+dt*(-3*k2[i]+6*k3[i]+9*k4[i])/16;
		for (i=2; i<n; i++)
			F[i]=(-v4[i+2]+12*v4[i+1]+36*v4[i]-44*v4[i-1]-3*v4[i-2])/60;
		F[1]=F[n]=0;
		tridiag(a, b, c, F, y, n+1);
		for (i=1; i<n+1; i++)
			k5[i]=-y[i]/dx;

	    for (i=1; i<n+1; i++)
			v5[i]=u[i]+dt*(k1[i]+4*k2[i]+6*k3[i]-12*k4[i]+8*k5[i])/7;
		for (i=2; i<n; i++)
			F[i]=(-v5[i+2]+12*v5[i+1]+36*v5[i]-44*v5[i-1]-3*v5[i-2])/60;
		F[1]=F[n]=0;
		tridiag(a, b, c, F, y, n+1);
		for (i=1; i<n+1; i++)
			k6[i]=-y[i]/dx;
		
		for (i=1; i<n+1; i++)
			u[i] += (7*k1[i] + 32*k3[i] + 12*k4[i] + 32*k5[i] + 7*k6[i])*dt/90;
		
		t += dt;
		u[0]=0;
		u[n+1]=0;
	}	 
}

void ComputScheme::Output(char* filename)
{
	ofstream fout;
	fout.open(filename,ios::out);
	for(int i=1;i<=n;i++)
		fout<<x[i]<<'\t'<<u[i]<<endl;
	fout.close();
}

int main()
{
	ComputScheme cp;
	cp.init();
//	cp.SCScheme();
//	cp.Output("SCS.plt");

//	cp.SUWScheme();
//	cp.Output("SUWS.plt");

//	cp.PadeScheme();
//	cp.Output("PadeS.plt");

//	cp.WSCScheme();
//	cp.Output("WSCS.plt");

	cp.FUWCScheme();
	cp.Output("FUWCS.plt");

    return 1;
}