burgers.cpp

burgers方程的数值解

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

#define UR 1
#define UL 0

using namespace std;

class Burgers{
private:
	int     n;
	double  dx;
	double* x;
	double* u;
	double* fn;
	double  time;
	double  cfl;
	double  dt;
	long    step;
	double  fun(double u);
public:
	Burgers();
	void init();
	void output(char* filename);
	void Lax_Friderichs();
	void Lax_Wendroff();
	void UpWind();
	void Godunov();
};

double Burgers::fun(double u)
{
	return 0.5*u*u;
}

Burgers::Burgers()
{
	int i;
	double pi=4.0*atan(1.0);
	n=400;
	dx=4.0/n;
	x=new double [n+2];
	u=new double [n+2];
	for(i=0;i<=n+1;i++)
		x[i]=-1.0+i*dx;
	for(i=0;i<=n+1;i++)
		if(x[i]>=0.0)
			u[i]=UR;
		else
			u[i]=UL;
}

void Burgers::init()
{
	fn=new double [n+2];
	step=0;
	cfl=0.9;
	time=2.0;
	dt=cfl*dx;
}

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

void Burgers::Lax_Friderichs()
{
	int i;
	double* v=new double [n+1];
	double t=0.0,r;
	while(t<time){
		if(t+dt>time) dt=time-t;
		t+=dt;
		step++;
		if(step%1000==0) 
			cout<<"t="<<t<<endl;
		r=dt/dx;
		for(i=1;i<=n;i++)			
			fn[i]=0.5*(fun(u[i])+fun(u[i-1]))-(u[i]-u[i-1])/r/2.0;
		for(i=1;i<n;i++)
			v[i]=u[i]-(fn[i+1]-fn[i])*r;
		for(i=1;i<n;i++)
			u[i]=v[i];
		u[0]=UL;
		u[n]=UR;
	}
}

void Burgers::Lax_Wendroff()
{
	int i;
	double* v=new double [n+1];
	double t=0.0,r,mu;
	while(t<time){
		if(t+dt>time) dt=time-t;
		t+=dt;
		r=dt/dx;
		step++;
		if(step%1000==0) 
			cout<<"t="<<t<<endl;
		for(i=1;i<=n;i++){
			if(fabs(u[i]-u[i-1])>0.00001)
				mu=r*(fun(u[i])-fun(u[i-1]))/(u[i]-u[i-1]);
			else
				mu=u[i-1];
			fn[i]=0.5*(fun(u[i])+fun(u[i-1]))-mu*mu*(u[i]-u[i-1])/r/2.0;
		}
		for(i=1;i<n;i++)
			v[i]=u[i]-(fn[i+1]-fn[i])*r;
		for(i=1;i<n;i++)
			u[i]=v[i];
		u[0]=UL;
		u[n]=UR;
	}
}

void Burgers::UpWind()
{
	int i;
	double* v=new double [n+1];
	double t=0.0,r,mu;
	while(t<time){
		if(t+dt>time) dt=time-t;
		t+=dt;
		r=dt/dx;
		step++;
		if(step%1000==0) 
			cout<<"t="<<t<<endl;
		for(i=1;i<=n;i++){
			if(fabs(u[i]-u[i-1])>0.00001)
				mu=r*(fun(u[i])-fun(u[i-1]))/(u[i]-u[i-1]);
			else
				mu=u[i-1];
			fn[i]=0.5*(fun(u[i])+fun(u[i-1]))-fabs(mu)*(u[i]-u[i-1])/r/2.0;
		}
		fn[0]=fn[n];
		fn[n+1]=fn[1];
		for(i=1;i<n;i++)
			v[i]=u[i]-(fn[i+1]-fn[i])*r;
		for(i=1;i<n;i++)
			u[i]=v[i];
		u[0]=UL;
		u[n]=UR;
	}
}

void Burgers::Godunov()
{
	int i;
	double* v=new double [n+1];
	double t=0.0,r;
	double a,b;
	while(t<time){
		if(t+dt>time) dt=time-t;
		t+=dt;
		r=dt/dx;
		step++;
		if(step%1000==0) 
			cout<<"t="<<t<<endl;
		for(i=1;i<=n;i++){
			a=u[i-1];
			b=u[i];
			if(a<=b){
				if(a>0)
					fn[i]=fun(a);
				else if(b<0)
					fn[i]=fun(b);
				else
					fn[i]=0.0;
			}
			if(a>b){
				if(a<0)
					fn[i]=fun(b);
				else if(b>0)
					fn[i]=fun(a);
				else{
					a=fun(a);
					b=fun(b);					
					fn[i]=(a>b)?a:b;
				}
			}
		}
		fn[0]=fn[n];
		fn[n+1]=fn[1];
		for(i=1;i<n;i++)
			v[i]=u[i]-(fn[i+1]-fn[i])*r;
		for(i=1;i<n;i++)
			u[i]=v[i];
		u[0]=UL;
		u[n]=UR;
	}
}

void main()
{
	Burgers p;
	p.init();
	p.output("init.plt");
	p.Godunov();
	p.output("GD.plt");
}