第十二讲__4.1傅里叶伪谱法(12).pptx

单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,2013/11/4,#,1.,确认修课名单；,2.,作业在下周四前交齐，可以重新提交；,3.,考试时间：,11,月,28,日下午,3,：,20-5,：,20,；,4.,答疑地点：科研楼东侧,403,室,5.,考前布置一个大作业，不用提交。,傅里叶伪谱法,有限差分法：直接从微分方程出发，用差商代替微商，,从而将微分方程直接转化为代数方程。,有限元法：首先将微分方程转化为等效积分形式，,应用变分原理或者加权余量法，,从而将定解问题转化为代数方程。,伪谱法,是一种高效、高精度计算地震波动场的数值方法，它利用,快速傅里叶变换对波动方程进行空间求导,。对介质参数平滑变化的非均匀介质，最小波长所需要的节点数比有限差分法少，因此和有限差分法相比，它的计算精度高、计算速度快，很适合于大规模模型的计算。,最初的伪谱法大多基于三角函数，随后出现了基于傅立叶函数、切比雪夫,函数的伪谱法，发展出更有效的自由边界和吸收边界的处理方式。,伪谱法是全局方法,只能用在空间域,。如果所求解的是含时微分方程，一般用简单的有限差分方法来处理时间域的计算。,傅里叶伪谱法将波场函数表示为傅里叶级数的展开形式，,空间域的求导转化为波数域的乘积计算，,时间域的求导用二阶差分格式给出。,（,Kosloff et al.,1984,）,在空间域，利用离散傅里叶变换求出位移对空间的导数。,考虑偏微分方程和边界条件,:,谱方法,其中,L,和,B,分别为线性微分算子和边界微分算子。,问题：求以上偏微分方程的数值解。,答案：数值解是满足边界条件，同时使以下残差较小的函数,如何度量残差的大小？,=,加权余量法,(Method of Weighted Residuals:MWR):,设,PN,是,Hilbert,空间的一个有限维子空间，在该空间内寻找,取该子空间的一组基函数作为试探函数,(trial functions),：,Hilbert,空间,：完备的内积空间，,希尔伯特空间为基于任意正交系上的多项式表示的傅立叶级数和,傅立叶变换提供了一种有效的表述方式，而这也是泛函分析的核心概念之一。,一般把偏微分方程形式称为强形式,(strong form),，由于以上微分方程在,U,内每一点都成立，因此总可以找到权函数,w,，使得下式成立：,上式称为偏微分方程的等效积分形式。,对于数值解，由于存在残差：,所以以上等效积分形式不能严格满足。由以上等效积分形式可以得,到以残差表示的近似的等效积分形式：,加权余量法是通过选择待定系数，强迫残差,R,在某种平均意义上等于零，残差的加权积分为零可得到一组方程，用来求解待定系数，进而得到原问题的近似解答。,谱方法,谱方法,有限差分方法：试探函数为重叠的局部,(local),低阶多项式。,有限单元方法：试探函数为局部,(local),光滑函数,(,固定阶数的多项式，只在,U,的子域内,非零,),。,谱方法：试探函数为,全局的,(global),光滑函数（如：,Fourier,级数）。,用配点法求解偏微分方程组,Fourier,伪谱方法,连续傅里叶变换,离散傅里叶变换,快速,Fourier,变换,(Fast Fourier transform,，,FFT),：,例：乘积运算次数,网格节点数 全矩阵,向量乘积,FFT,比值,1D(nx=512)2.6x10,5,9.2x10,3,28.4,1D(nx=2096)4.4x10,6,4.6x10,4,94.98,1D(nx=8384)7.0 x10,7,2.2x10,5,312.6,全矩阵,-,向量相乘的运算次数与,FFT,运算次数的比值可以看作使用,FFT,时的,加速比,。,李本文,李本文,function df=sderld(f,dx),%SDERLD(f,dx)spectral derivative of vector,nx=max(size(f);,%initialize k,kmax=pi/dx;,dk=kmax/(nx/2);,for i=1:nx/2;,k(i)=i*dk;,k(nx/2+i)=-kmax+i*dk;,end,k=sqrt(-1)*k;,%FFT and IFFT,ff=fft(f);,ff=k.*ff;,df=real(ifft(ff),;,clc;,clear;,figure(color,w);,for i=1:100;,t(i)=(i-1)*pi/100;,f(i)=sin(2*t(i);,end,nx=max(size(f);,%initialize k,dx=pi/100;,kmax=pi/dx;,dk=kmax/(nx/2);,for i=1:nx/2;,k(i)=i*dk;,k(nx/2+i)=-kmax+i*dk;,end,k=sqrt(-1)*k;,%FFT and IFFT,df=sderld(f,dx);,plot(t,f,t,df);,xlabel(time);,ylabel(Signal);,legend(f,df);,吴宝年等，,2012,吴宝年等，,2012,吴宝年等，,2012,一阶导数发生数值振荡，二阶导数没有发生数值振荡,吴宝年等，,2012,吴宝年等，,2012,一阶导数没有发生数值振荡，二阶导数发生数值振荡,采用交错网格计算一阶导数，采用中心网格计算二阶导数，该方法可以有效消除由于数值振荡产生的干扰伪波。,吴宝年等，,2012,吴宝年等，,2012,垂直分量,吴宝年等，,2012,水平分量,龙桂华等，,2009,龙桂华等，,2009,时间的离散,在数值计算过程中，介质沿水平和纵向分割成,255255,个网格，相邻格点之间的间距为,10m,震源位于地表下,190m,处，且位于模型水平方向的正中心，为主频,25,的,Ricker,子波，其中子波的时间采样间隔为,1ms.,龙桂华等，,2009,龙桂华等，,2009,龙桂华等，,2009,利用伪谱法模拟横向非均匀全球模型中的,SH,波场,中国科学,:,地球科学，,2012,年第,42,卷第,1,期,:140 148,王彦宾,Hiroshi TAKENAKA,基于傅里叶伪谱法,(Fourier pseudospectral method,PSM),，,发展了模拟具有任意横向非均匀结构的全地球模型中,SH,波传播的算法。所采用的模型为通过地球大圆选取的一个二维全地球剖面，,所求解的弹性波动方程定义在二维柱坐标系下。波动方程中空间微分的数值求解通过波数域的乘积运算实现，,空间域与波数域之间的变换通过快速傅里叶变换,(fast Fourier transformation,，,FFT),进行,.,由于,PSM,具有精度较高、计算速度较快的特点，所以对于非均匀结构可以近似为在垂直剖面方向上对称的模型，本算法能够较精确地计算较高频率的全球,SH,波场的理论地震图,.,王彦宾,Hiroshi TAKENAKA,2012,王彦宾,Hiroshi TAKENAKA,2012,王彦宾,Hiroshi TAKENAKA,2012,王彦宾,Hiroshi TAKENAKA,2012,周期性边界条件,（,Periodic Boundary Conditions,PBC,）是边界条件的一种，反映的是如何利用边界条件替代所选部分（系统）受到周边（环境）的影响。可以看作是如果去掉周边环境，保持该系统不变应该附加的条件，也可以看作是由部分的性质来推广表达全局的性质。,