advection2.cpp

求解对流微分方程：u_t+a(x,t)u_x=0

/*****************************************************************************
    求解微分方程：
	    u_t+a(x,t)u_x=0
	其中
	                   1+x^2
		a(x,t)=---------------------
		         1+2*x*t+2*x^2+x^4
	初始条件
                / 1,   0.2<=x<=0.4
        u(x,0)=<
		        \ 0,   else
	边界条件
	    u(x,t)=u(x*,0)
	其中
	               t   
		x*=x - ---------
		         1+x^2
*****************************************************************************/
#include <iostream>
#include <fstream>
#include <cmath>

using namespace std;

class LinAdv{
private:
	double* u;
	double* x;
	int     n;
	double  cfl;	
	double  dx;
	double  dt;
	double  tend;
	double  a(double x,double t){return (1+x*x)/(1+2*x*t+2*x*x+x*x*x*x);};
	double  da(double x,double t){return -2*x*(1+x*x)/(1+2*x*t+2*x*x+x*x*x*x)/(1+2*x*t+2*x*x+x*x*x*x);}
public:
	LinAdv();
	void Lax_Friedrichs();
	void UpWind();
	void Lax_Wendroff();
	double getError();
	void output();
};

LinAdv::LinAdv()
{
	int i;
	n=399;
	dx=1.0/(n+1);
	x=new double [n+2];
	for(i=0;i<n+2;i++)
		x[i]=i*dx;
	u=new double [n+2];
	for(i=0;i<n+2;i++){		
		u[i]=0.0;
		if(x[i]>=0.2&&x[i]<=0.4)
			u[i]=1.0;
//		u[i]=exp(-10*(4*x[i]-1)*(4*x[i]-1));
	}
	cfl=1.0;
	dt=cfl*dx;
	tend=1.0;
}

void LinAdv::output()
{
	ofstream fout;
	fout.open("result.plt",ios::out);
	for(int i=0;i<n+2;i++)
		fout<<x[i]<<'\t'<<u[i]<<endl;
	fout.close();
}

void LinAdv::UpWind()
{
	int i;
	double* v=new double [n+1];
	for(i=0;i<n+1;i++)
		v[i]=0.0;
	double t=0.0;
	while(t<tend){
		if(t+dt>tend) dt=tend-t;
		t+=dt;
		cout<<t<<endl;
		for(i=1;i<=n;i++)
			v[i]=u[i]-a(x[i],t)*dt*(u[i]-u[i-1])/dx;
		for(i=1;i<=n;i++)
			u[i]=v[i];
		u[0]=0.0;
		u[n+1]=0.0;
	}
}

void LinAdv::Lax_Friedrichs()
{
	int i;
	double* v=new double [n+1];
	for(i=0;i<n+1;i++)
		v[i]=0.0;
	double t=0.0;
	while(t<tend){
		if(t+dt>tend) dt=tend-t;
		t+=dt;
		cout<<t<<endl;
		for(i=1;i<=n;i++)
			v[i]=0.5*(u[i+1]+u[i-1])-0.5*a(x[i],t)*dt*(u[i+1]-u[i-1])/dx;
		for(i=1;i<=n;i++)
			u[i]=v[i];
		u[0]=0.0;
		u[n+1]=0.0;
	}
}

void LinAdv::Lax_Wendroff()
{
	int i;
	double* v=new double [n+1];
	for(i=0;i<n+1;i++)
		v[i]=0.0;
	double t=0.0;
	while(t<tend){
		if(t+dt>tend) dt=tend-t;
		t+=dt;
		cout<<t<<endl;
		for(i=1;i<=n;i++)
			v[i]=u[i]-0.5*a(x[i],t)*dt*(u[i+1]-u[i-1])/dx
			         -0.25*dt*dt*da(x[i],t)*(u[i+1]-u[i-1])/dx
					 +0.5*dt*dt*a(x[i],t)*(a(x[i]+0.5*dx,t)*(u[i+1]-u[i])/dx-a(x[i]-0.5*dx,t)*(u[i]-u[i-1])/dx)/dx;
		for(i=1;i<=n;i++)
			u[i]=v[i];
		u[0]=0.0;
		u[n+1]=0.0;
	}
}

double LinAdv::getError()
{
	double* v=new double [n+2];	
	double  xx;
	for(int i=0;i<n+2;i++){
		xx=x[i]-tend/(1+x[i]*x[i]);
		if(xx>=0.2&&xx<=0.4)
			v[i]=1.0;
		else
			v[i]=0.0;
	}
	ofstream fout;
	fout.open("true.plt",ios::out);
	for(i=0;i<n+2;i++)
		fout<<x[i]<<'\t'<<v[i]<<endl;
	fout.close();
	double error=0.0;
	for(i=1;i<n;i++)
		error+=(u[i]-v[i])*(u[i]-v[i]);
	return sqrt(error*dx);
}

void main()
{
	LinAdv P;
	P.Lax_Friedrichs();
	P.output();	
	cout<<P.getError()<<endl;
}