《最新二维波动方程第一类吸收边界条件c++实现代码.doc》由会员分享，可在线阅读，更多相关《最新二维波动方程第一类吸收边界条件c++实现代码.doc（5页珍藏版）》请在三一文库上搜索。
1、二维波动方程第一类吸收边界条件C + +实现代码#include "stdafx.h”#include <iostreani>#include <fstreani>#include <cmath>#include<string>using namespace std;const double pi=4*atan(1.0);double freq=45;double sb=7.45;double tl=2*pi/(sb*4);double source(double t)/double t2=0.0;if(t<=tl) return
2、 (sin(sb*4*t-pi/2)+1)/10;else double tep=O.O; 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 upxl,double lowx 1)intconst int Csize=300;double des=0;double stepx 1 =abs(upx 1 lowx 1 )/(Csize-l);/double te=sqrt(static
3、_cast<double>(3.0/8.0);double stept=sqrt(static_cast<double>( 1.0/2.0)*stepx 1/2.0;/ int tn=static_cast<int>(upt/stept);double r=stept/stepx 1;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;i<Csize;i+)u_curren
4、ti=new double Csize;u_oldi=new doubleCsize; u_pasti=new doubleCsize;for(i=0;i<Csize;i+)for(j=0;j<Csize;j+)u-currentij=O;u.oldiUJ=0;u.pastij=O;1double ckCsizeCsize;int flag=O;for(j=0:j<Csize;j+)for(i=0:i<Csize;i+)if(j<i)ckiUl=4; else ckij=l;1string str;cout«"n输入保存文件名:”；cin
5、87;str;ofstream fout(str); if(!fout) cout«Hn 不能打开文件H«str«endl;exit( 1);Im=0;double fO=2.O/(stept*3O);double t0=4.0/f0;while(m< 1500&&(m+)*stept+lowt)<upt) for(i=0;i<Csize;i+) for(j=0;j<Csize;j+)if(i!=0&&(i!=Csize-l)&&(j!=0&&j !=(Csize-1 )/co
6、pe with the internal points of the interesting domain/if(i=(Csize/2)&&j=(Csize/2)u_currentij=(卢 r/cki|j)*(u_oldi+lj+u_oldij+l+u_oldilj+u_oldiH) 2*u_oldij*(2*i升/cki卜 lO)-u_pastij/+stept*stept*source(m*stept+lowt);/stept*stept*exp(-fO*fO*(m*stept-tO)*(m*stept- tO);/stept*stept*source(m*stept+l
7、owt);/elseu_currenti j =(r*r/cki j)*(u_oldi+1 j+u_oldi j+1 +u_oldi- lj+u_oldijJ)-2*u_oldi|j*(2*i叫ycki|j 卜 l0) u_pastij;/iimij二卢卢 ckij*ckij*(umli+lj+um-lij+l+um-/-um-2iUJ;u_currentij=r* 卢 ckij*ckij*(u_oldi+llj+u_okiij+l+u_oldi|j - l+u_oldilj)2*(2打*卢 ckij*ckij 卜 1 尸u_oldij/-u_pastij;/u_currentij=(powe
8、r(r,2.0)/cki|j)*(u_oldi+lj+u_oldi|j+l+u_oldi- lU+u_oldiU-l)-2*u_oldiUl*(2*pow(r,2.0)/ckiU-l)-u_pasti|j;u_currentij=(r*r/ckij)*(u_oldi+l|j+u_oldi|j+l+u_oldi-lj+u_oldij-l)- 2*u_oldij*(2水卢 i7ckij 卜 l0)7_pastij;1/umO|j=umlU;/u m Csize-1 j=um Csize-2 j; /umiO=umil;/ u m i Csize-1 =u m i Csize-2; if(i=Csi
9、ze/2)u_currentCsize/2 1 =u_currentCsize/20+stepx 1 *exp(-fO*fO*(m*stept- tO)*(m*stept-tO);/stepx 1 *source(m*stept+lowt)+u m Csize/2 0 ;stepx 1 *source(m*stept+lowt)elseu_currenti0=u_oldi0+u_oldil-u_pastil+(1.0/sqrt(cki0)*r*(u_oldil- u_oldi0-u_pasti2+u_pasti 1 )；/top side absorobing boundary uhni0=u
10、mil;/u_currenti0=u_oldi0+u_oldil-u_pastil+cki0*r*(u_oldil- u_oldi0-u_pasti2+u_pasti 1 )；/top side absorobing boundary uhni0=umil;/i】_cuirentiCsizel=ii_oldiCsize1 +u_oldiCsize-2-u_pastiCsize-2-( 1.0/sqrt(ckiCsize-l )*r*(u_oldi Csize-1 - u_oldiCsize-2-u_pastiCsize-2+ u_pastiCsize-3);/ down side absorb
11、ing boundaryu_current0j=u_old0Ij+u.oldlj-u_past 1 j+( 1.0/sqrt(cki Csize-1 )*r*(u_old 1 j-u_old0 j-u_past2j+u_past 1 j); /Left side absorbing boundaryu_currentCsize-l j =u_oldCsize-lj+u_oldCsize-2j 卜 u_pastCsize2j 卜(1.0/sqrt(ckiCsizel)*r*(u_oldCsiz少 lU-u_oldCsize-2Ul-u_pastCsize-2j+u_pastCsize 3j);/
12、Right side absorbing boundaryIif(m%10=0)for(i=0;i<Csize;i+)for(j=0:j<Csize;j+)fout«H H«u_currentlji; /fout«H M«ckij； fout«endl;fout«endl;1for(i=0;i<Csize;i+)for(j=0;j<Csize:j+) u.pastij=u_oldi|j; u_oldij二 u_cuiTentij; u_currentij=O;/ m+; fout.close();for(i=0;i<Csize;i+) delete u_currenti,u_oldi,u_pasti; delete u_current,u_old,u_past;1void main ()主函数double upt=17.8*tl*10;double lowt=0.upx 1 =3.0,lowx 1 =-3.0;/ const int Csize=100;/ const int tn=100;/double upt=l;update_Vn(upt, lowt,upx 1 ,lowx 1);1