资源描述：
《伪谱法求解波动方程算法的VC_实现_杨三女》由会员上传分享，免费在线阅读，更多相关内容在行业资料-天天文库。
1、第26卷第2期华东地质学院学报Vol126No122003年6月JOURNALOFEASTCHINAGEOLOGICALINSTITUTEJun.2003伪谱法求解波动方程算法的VC++实现杨三女陈永清陈建国(中国地质大学资源学院,武汉430074)摘要:伪谱法利用傅立叶变换和有限差分求解波动方程,是三维波场正演模拟的重要方法之一。用VC++语言实现这一系统必须解决二阶空间偏导数、二阶时间偏导数的离散化问题。我们利用实信号的傅立叶变换和反傅立叶变换的基本性质即共轭性解决了空间偏导数的离散化问题,对于时间偏导数,我们则采用
2、二阶中心差分的方法解决。最后模拟结果证明是有效的。关键词:伪谱法;傅立叶变换;有限差分中图分类号:P628+.3文献标识码:A文章编号:1000-2251(2003)02-183-04222波动方程的解法有解析法和数值法两种。一25u5u5u22C2+2+2$t+s$t(2)般说来,波动方程的解很难用解析的形式表示出5x5y5z来,于是人们便寻求用数值方法来求解,随着计算(2)式中,对于时间导数,我们采用二阶中心差机的广泛应用,这种方法的实现已成为可能。波动分法求解;对于空间导数采用快速傅立叶变换方法方程的数值解法有多种
3、,如有限差分法、有限元法、求解。克希霍夫积分法和伪谱法。有限差分法的算法简1.2.2傅立叶变换法求空间导数的原理单快速,但难以克服频散效应(马永旺等,2001),而设函数g(x1,x2,,,,xm)的傅立叶正变换为要解决频散问题,需加密数值计算的网格,但这势G(kx,kx,,,,kx)则g(x1,x2,,,,xm)=12m必导致计算量增加,效率下降。有限元法模拟是目+]+]+]i(kx+kx+,kx)前最精确的一种正演模拟方法,但它对计算机内存,G(kx,kx,,kx)exi1x22xmmQ-]Q-]Q-]12m要求高,
4、计算量大,让人难以承受。因此,选择一种dkxdkx,dkxm(3)(刘洋等,2000)。§式两边对x112既能精确计算,又有较高计算效率的方法就显得非取n1阶导数,对xm取nm阶导数,,,,对x2取n2常必要。伪谱法正好符合这种要求,伪谱法是一种阶导数,则有:有效的数值模拟方法,权衡精度和效率,有其它2(n,n,,n)种方法不可替代的优点。dg12m(x1,x2,,,xm)nnn=dx11dx22,dxmm1伪谱法的基本原理+]+]+]nnn,(-ikx)1(-ikx)2,(-ikx)m*Q-]Q-]Q-]12m1.1三
5、维声波波动方程i(kx+kx+,+kx)2222G(kxkx,,kx)ex11x22xmmdkxdkx,dkx15u5u5u5u12m12m22=2+2+2+s(1)C5t5x5y5z(4)¥式中,t为声波传播时间;c=c(x,y,z)为若令x=(x1,x2,,,xm),K=(kx,kx,12声波传播速度;u=u(x,yx,z,t)为波场压力或位,kx),则用傅立叶变换求解空间导数的过程为m移;S=S(x,y,z,t)为震源项(庄东海等,1998;FFTnnnResheFM,1998)。g(x)G(K)v(-ikx)1(
6、-ikx)2(-ikx)m12mFFT-1(n,n,,,n)1.2波动方程伪谱法的数值解法dg12m(x)nnn。1.2.1声波传播伪谱法数值模拟的递推方程dx11dx21,dxmmu(t+$t)=2u(t)-u(t-$t)+1.2.3空间导数的计算根据式§、¨中提到的原理,我们知道:收稿日期:2002-10-08作者简介:杨三女(1979)),男,在读研究生,现主要从事地球探测与信息技术研究。184华东地质学院学报2003年2的二阶空间偏导数ui2Pixx(i,j,k)如下^uxx(i,j,k)=^u(i,j,k)(5
7、)Dxg_Uxx(,)//实现uxx(i,j,k)的函数2i2Pj{^uyy(i,j,k)=^u(i,j,k)(6)Dy,,//缩略2i2Pk^uzz(i,j,k)=^u(i,j,k)(7)T1=Real(fft3D(u(i,j,k)));//将u(i,j,k)经过Dz2三维傅立叶变换后得到的实部赋给T15u(i,j,k)其中^uxx(i,j,k)表示2的离散值,T2=Img(fft3D(u(i,j,k)));//将u(i,j,k)经过5x52u(i,j,k)三维傅立叶变换后得到的实部赋给T2^uyy(i,j,k)表示2
8、的离散值,^uzz(i,j,k)5yif(iFM/2)//M为u(i,j,k)在x方向上的抽25u(i,j,k)样点数表示2的离散值,^u(i,j,k)为波场分量5z{的离散值,i为虚部单位,Dx为模型中x方向上的长P2221=(2*P)*(i)T1/((M-1)*$x);//度,Dy为模型中y方向上的长度,Dz为模
