nls.cpp

三维非线性薛定鄂方程的有限差分解法的源程序

# include <math.h>
# include <stdio.h>
# include <fstream.h>
# include <iostream.h>
void chafen2D(double * data, int nt, int nhx, int nhy,
			  double tau, double hx, double hy,
			  double reb, double imb, double nlsre, double nlsim)
			  
{
	double tbhxre=tau*reb/hx/hx, tbhxim=tau*imb/hx/hx;
    double tbhyre=tau*reb/hy/hx, tbhyim=tau*imb/hy/hy;
	double tnlsre=tau*nlsre, tnlsim=tau*nlsim;
	int it, ihx, ihy, nytm, nxtm;
	double tmxim, tmxre, tmyre, tmyim, tmnre, tmnim, temp;
	double e20, e21, e22, e11, e12;
	for(it=1; it<nt;it++)
	{
		for(ihy=1;ihy<nhy;ihy++)
		{
			nytm=(ihy+ihy)*nhx;
			for(ihx=1;ihx<nhx;ihx++)
			{
				if((it+1+ihy+ihx)%2==1)continue;//ou
				nxtm=ihx+ihx;
				if(ihy==nhy-1||ihx==nhx-1)
				{
					data[nytm+nxtm]=0.0;
					data[nxtm+nxtm+1]=0.0;
					continue;
					
				}
				tmyre=(data[nytm-2*nhy+nxtm]-2*data[nytm+nxtm]+data[nytm+2*nhy+nxtm]);
				tmyim=(data[nytm-2*nhy+nxtm+1]-2*data[nytm+1+nxtm]+data[nytm+2*nhy+nxtm+1]);
				tmxre=(data[nytm+nxtm-2]-2*data[nytm+nxtm]+data[nytm+nxtm+2]);
				tmxim=(data[nytm+nxtm-1]-2*data[nytm+nxtm+1]+data[nytm+nxtm+3]);
				e20=(data[nytm+nxtm]*data[nytm+nxtm]+data[nytm+nxtm+1]*data[nytm+nxtm+1]);
				e21=(data[nytm+nxtm-2]*data[nytm+nxtm-2]+data[nytm+nxtm-1]*data[nytm+nxtm-1]);
				e22=(data[nytm+nxtm+2]*data[nytm+nxtm+2]+data[nytm+nxtm+3]*data[nytm+nxtm+3]);
				
				e12=(data[nytm-2*nhy+nxtm]*data[nytm-2*nhy+nxtm]+data[nytm-2*nhy+nxtm+1]*data[nytm-2*nhy+nxtm+1]);
				e11=(data[nytm+2*nhy+nxtm]*data[nytm+2*nhy+nxtm]+data[nytm+2*nhy+nxtm+1]*data[nytm+2*nhy+nxtm+1]);
				temp=(e20+(e21+e22)/2.0+(e12+e11)/2.0)/3.0;
				tmnre=temp*data[nytm+nxtm];
				tmnim=temp*data[nytm+nxtm+1];
				temp=data[nytm+nxtm]-tbhxre*tmxim-tbhyre*tmyim-tbhxim*tmxre-tbhyim*tmyre+tmnre*tnlsre-tmnim*tnlsim;
				data[nytm+nxtm+1]=data[nytm+nxtm+1]+tbhxre*tmxre+tbhyre*tmyre-tbhxim*tmxim-tbhyim*tmyim+tmnre*tnlsim+tmnim*tnlsre;
                data[nytm+nxtm]=temp;
			}
			
		}
		for(ihy=1;ihy<nhy;ihy++)
		{
			nytm=(ihy+ihy)*nhx;
			for(ihx=1;ihx<nhx;ihx++)
			{
				if((it+1+ihy+ihx)%2==0)continue;//qi
				nxtm=ihx+ihx;
				if(ihy==nhy-1||ihx==nhx-1){
					data[nytm+nxtm]=0.0;
					data[nytm+nxtm+1]=0.0;
					continue;
				}
				tmyre=(data[nytm-2*nhy+nxtm]-2*data[nytm+nxtm]+data[nytm+2*nhy+nxtm]);
				tmyim=(data[nytm-2*nhy+nxtm+1]-2*data[nytm+1+nxtm]+data[nytm+2*nhy+nxtm+1]);
				tmxre=(data[nytm+nxtm-2]-2*data[nytm+nxtm]+data[nytm+nxtm+2]);
				tmxim=(data[nytm+nxtm-1]-2*data[nytm+nxtm+1]+data[nytm+nxtm+3]);
				e20=(data[nytm+nxtm]*data[nytm+nxtm]+data[nytm+nxtm+1]*data[nytm+nxtm+1]);
				e21=(data[nytm+nxtm-2]*data[nytm+nxtm-2]+data[nytm+nxtm-1]*data[nytm+nxtm-1]);
				e22=(data[nytm+nxtm+2]*data[nytm+nxtm+2]+data[nytm+nxtm+3]*data[nytm+nxtm+3]);
				e12=(data[nytm-2*nhy+nxtm]*data[nytm-2*nhy+nxtm]+data[nytm-2*nhy+nxtm+1]*data[nytm-2*nhy+nxtm+1]);
				e11=(data[nytm+2*nhy+nxtm]*data[nytm+2*nhy+nxtm]+data[nytm+2*nhy+nxtm+1]*data[nytm+2*nhy+nxtm+1]);
				temp=(e20+(e21+e22)/2.0+(e12+e11)/2.0)/3.0;
				tmnre=temp*data[nytm+nxtm];
				tmnim=temp*data[nytm+nxtm+1];
				temp=data[nytm+nxtm]-tbhxre*tmxim-tbhyre*tmyim-tbhxim*tmxre-tbhyim*tmyre+tmnre*tnlsre-tmnim*tnlsim;
				data[nytm+nxtm+1]=data[nytm+nxtm+1]+tbhxre*tmxre+tbhyre*tmyre-tbhxim*tmxim-tbhyim*tmyim+tmnre*tnlsim+tmnim*tnlsre;
                  data[nytm+nxtm]=temp;
			}
		}
	}
	
}
void gauss(double * data, int nhx, int nhy, double hx, double hy)
{
	double x,y,tm;
	int ihy, ihx, nytm;
	for(ihy=0;ihy<nhy;ihy++)
	{
		nytm=2*ihy*nhy;
		y=hy*(ihy-nhy/2);
		for(ihx=0;ihx<nhx;ihx++)
		{
			x=hx*(ihx-nhx/2);
			tm=pow(x*x+y*y,1);
			data[nytm+2*ihx]=exp(-tm);
			data[nytm+2*ihx+1]=0.0;
		}
	}
}
bool print2D(double*data, int nn, double delta, char*filename, int pnum)
{
	//open a file for outpur the profile
	double xxx, yyy, inten;
	int jx,jy,m,k;
	ofstream outfile(filename);
	//outfile.open
	if(!outfile){
		cout<<"can not open"<<filename<<"for output!\n";
		return false;
	}
	for (jy=0;jy<=nn-1;jy++)
	{
		if(jy%pnum!=0) continue;
		m=jy*(nn+nn);
		yyy=(float)((jy-nn/2)*delta);//jy is real index
		for(jx=0;jx<=nn-1;jx++)
		{
			if(jx%pnum!=0)continue;
			k=jx+jx+m;
			xxx=(float)((jx-nn/2)*delta);//jx is real index
			inten=(float)(data[k+2]*data[k+2])+data[k+1]*data[k+1];
			//phase=(float)(atan2(data[k],data[k+1]))
			outfile<<xxx<<" "<<yyy<<" ";
			outfile<<inten<<" "<<'\n';//<<phase

		}
	}
	outfile.close();
	return true;

}
main()
{
	int nt=120;
	int nhx=60,nhy=60;
	double * data;
	data=new double[2*nhx*nhy];
	double hx=0.1,hy=0.1,tau=0.1;
	double reb=0.0, nlsre=0.0,imb=0.005,nlsim=imb*1.0e1;
	gauss(data,nhx,nhy,hx,hy);
	print2D(data,nhx,hx,"Dim2d0.dat",1);
	chafen2D(data,nt,nhx,nhy,tau,hx,hy,reb,imb,nlsre,nlsim);
	print2D(data,nhx,hx,"Dim2d4.dat",1);
	return 0;

}