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