二维波动方程第一类吸收边界条件c++实现代码.docVIP

#include stdafx.h #include iostream #include fstream #include cmath #includestring using namespace std; const double pi=4*atan(1.0); double freq=45; double sb=7.45; double t1=2*pi/(sb*4); double source(double t) { //double t2=0.0; if(t=t1) return (sin(sb*4*t-pi/2)+1)/10; else{ double tep=0.0; return tep;} //return ((1-2*pi*pi*freq*freq*t*t)*exp(-pi*pi*freq*freq*t*t)+1);//Ricker子波 } void update_Vn(double upt,double lowt,double upx1,double lowx1) { int i,j,m; const int Csize=300; double deg=0; double stepx1=abs(upx1-lowx1)/(Csize-1); //double te=sqrt(static_castdouble(3.0/8.0)); double stept=sqrt(static_castdouble(1.0/2.0))*stepx1/2.0;// int tn=static_castint(upt/stept); double r=stept/stepx1; double **u_current,**u_old,**u_past; u_current=new double *[Csize]; u_old=new double*[Csize]; u_past=new double*[Csize]; for(i=0;iCsize;i++) { u_current[i]=new double [Csize]; u_old[i]=new double[Csize]; u_past[i]=new double[Csize]; } for(i=0;iCsize;i++) for(j=0;jCsize;j++) { u_current[i][j]=0; u_old[i][j]=0; u_past[i][j]=0; } double ck[Csize][Csize]; int flag=0; for(j=0;jCsize;j++) { for(i=0;iCsize;i++) { if(ji) ck[i][j]=4; else ck[i][j]=1; } } } string str; cout\n 输入保存文件名:; cinstr; ofstream fout(str); if(!fout) { cout\n 不能打开文件strendl;exit(1); } m=0; double f0=2.0/(stept*30); double t0=4.0/f0; while(m1500((m++)*stept+lowt)upt) { for(i=0;iCsize;i++) for(j=0;jCsize;j++) { if((i!=0(i!=Csize-1))(j!=0j!=(Csize-1)))//cope with the internal points of the interesting domain { //if(i==(Csize/2)j==(Csize/2)) u_current[i][j]=(r*r/ck[i][j])*(u_old[i+1][j]+u_old[i][j+1]+u_old[i-1][j]+u_old[i][j-1])-2*u_old[i][j]*(2*r*r/ck[i][j]-1.0)-u_past[i][j] //+stept*stept*source(m*stept+lowt);//stept*stept*exp(-f0*f0*(m*stept-t0)*(m*stept-t0));//stept*stept*source(m*stept+lowt); //elseu_current[i][j]=(r*r/ck[i][j])*(u_old[i+1][