数值模拟1(FDM)课件

1、地球物理数值模拟方法地球探测与信息技术专业研究生课程一、 课程性质、目的和任务地球物理数值模拟是地球探测与信息技术专业应用地球物理方向的一门必修课。通过学习，使学生了解和掌握用有限元、边界元和有限差分法解地球物理中的边值问题。二、 教学要求、形式和方法要求学完以后，学生能够用有限元、边界元、有限差分中的一种方法编制某一实际问题的数值模拟计算程序，并撰写读书报告。在教学形式上以自学为主，课堂教学为辅，教学方法上则采用教师引导性讲学与学生结合问题的自学。 三、 教学内容第四讲：边界元结合地形影响的边界元数值模拟讲述：积分方程的导出，单元剖分，插值方法，数值积分方法，线性方程解法。第五讲：成绩评定对

2、所学知识，特别是自己所研究解决的问题作学术报告。 四、教材及参考书目1、地球物理中的有限单元法，徐世浙，科学出版社，1994；2、地球物理中的边界单元法，徐世浙，科学出版社，1995；3、电子计算机在电法勘探中的应用，罗延钟，张桂青，武汉地质学院出版社，1987；4、电法勘探数值模拟技术，周 襄，钟本 等编著，四川科学出版社，1986；5、电法勘探用边界元法，田宪谟，黄兰珍，地质出版社，1990；6、The Finite Element Method, O.C.Zienkiewicz, R.L.Taylor, Mcgraw-Hill International Editions, 1989；7

3、、阮百尧，熊彬，徐世浙，2001，三维地电断面电阻率测深有限元数值模拟，地球科学，26（1），73-77；8、阮百尧，2001，视电阻率对模型电阻率的偏导数矩阵计算方法，地质勘探，37（6），39-41；9、阮百尧，2001，三角单元剖分电导率分块连续变化点源二维电场有限元数值模拟，广西科学，8（1），1-3；10、阮百尧，徐世浙，1998，电导率分块线性二维地电断面电阻率测深有限元数值模拟，地球科学,23,303-307。第一讲：概论正演问题：已知场源和地下物性分布，计算地球物理场的空间和时间分布。对直流电阻率法而言，正演问题是已知供电电极位置、电流大小，地下介质的电阻率（或电导率）分布，求

4、地表和地下的电位。第一讲：概论正演的作用：勘探设计的依据：针对不同的地质条件、地质目标，确定用何种地球物理方法、测量的方式、测网、极距的大小，了解可能得到的异常；实测数据解释的基础：定性分析、一维、二维和三维的解释。剖面I剖面II剖面I还是剖面II？第一讲：概论正演的方法物理模拟法：水槽、土槽、导电纸、电阻网络数学模拟法：数值计算法：有限差分法、有限元法、积分方程法（边界元法）解析法：分离变量法、镜像法、类比法数学模拟法：大部分的地球物理过程都能以数学的形式加以描述。 给定模型参数，数据观测方式，建立数学物理模型，如：积分方程；微分方程、初始条件、边界条件变分方程用解析或数值模拟方法求解数学物

5、理模型，得到预测数据。数值模拟方法常用的数值模拟方法有：有限差分法（FiniteDifferenceMethodFDM）有限元法（Finite Element MethodFEM ）积分方程法 边界元法（Boundary Element MethodBEM）有限差分方法(FDM)是计算机数值模拟最早采用的方法，至今仍被广泛运用。该方法将求解域划分为差分网格，用有限个网格节点代替连续的求解域。有限差分法以Taylor级数展开等方法，把控制方程中的导数用网格节点上的函数值的差商代替进行离散，从而建立以网格节点上的值为未知数的代数方程组。该方法是一种直接将微分问题变为代数问题的近似数值解法，数学概念

6、直观，表达简单，是发展较早且比较成熟的数值方法。有限差分法的缺点是必需进行整个区域的剖分，并且要求网格比较规则，空间网格最好为直角网格。有限差分法的特点有限元法的特点有限元方法（FEM）最早应用于结构力学，其基本求解思想是借助于变分原理或加权余量法，将微分方程转变为泛函极值问题，然后把计算域划分为有限个互不重叠的单元，在每个单元内，将连续函数的泛函，离散成单元节点上函数值的泛函。根据泛函取极值的条件，得到各节点的函数值应满足的线性代数方程组。解代数方程组，得到各节点的函数值。FEM法的优点是求解规范（甚至可用软件产生FEM源码），并且可以模拟所有形体。FEM的缺点也是必需进行整个区域的剖分，与

7、FDM相比，相对比较复杂。积分方程法（边界元法）的特点边界元法（BEM）通过格林公式，将偏微分方程表示的边值问题转变成区域边界上的积分方程，然后将边界划分为有限个互不重叠的单元，进行数值积分。BEM是在有限单元法以后发展起来的一种数值计算方法，可以说是解边界积分方程的有限元法。BEM的优点是将区域问题转变为区域边界问题，使求解问题的维数降低，三维变为二维、二维变为一维。BEM的缺点是不规范，当区域内有不同数量的介质体时，需重新编制程序；另外还有格林公式中的基本解很难找，线性方程组的系数矩阵为满秩矩阵，这些都导致BEM使用不如其他两种方法广泛。数学补充：关于偏微分方程的一些基本概念几个典型方程拉

8、普拉斯算符哈密顿算符几个典型方程波动方程在一些声学、电磁学、光学和力学的波动问题中常常出现这类方程。扩散方程、热传导方程其中， 是扩散过程中某种物质的浓度，或是固体的传热过程中在x处、t时刻的温度。定解条件对于一个n阶常微分方程，常常可以将其解写成依赖于n个任意常数的通解形式。但是对于偏微分方程，情况要复杂得多，一般很难用通解的形式表示。我们都是特定条件下求方程的解，这样的条件称为定解条件。如果一个偏微分方程定解问题在某个函数集合中存在唯一的解，而且在定解条件的原始资料有微小变化时，在某种义下解也仅有微小的变化，我们说解关于定解条件是稳定的。如果一个定解问题解的存在性、唯一性和稳定性都成立，就

9、称定解问题是适定的。我们所讨论的定解问题都是适定的。第二讲：有限差分法差分一维有限差分二维有限差分时间域有限差分例题第1章 差分一阶差分xx+xx-x第2章 一维波动方程的有限差分法把上面得到的差分算子代替弹性波动方程中的微分算子，波动方程可以近似表示为：或者这个有限差分方程的截断误差是初始条件为：边界条件为：有限差分方程和离散形式的初值边界条件一起构成了波动方程的差分格式。如果已知在时间 和 ，空间 和 处的位移值，则可以根据差分波动方程中时间 和空间 处的位移值的显式表达式，通过简单数学运算实现对整个弹性波波场的时间空间逐层推进计算。这种有限差分格式称为显式格式。在解题过程中，需要利用初始

10、条件，设 时， 用中心差分可表示为：于是得 把它代入上述方程中的最后一项，并取 j0， 即对应t0 ，得出举例初位移不为零，初速为零时的例子：设初位移为：其余令 ， ，将微分方程改写成差分方程，即有其中，初始条件可以表示成此处，上标为一时为t0网格剖分有限差分方法用差分算子代替微分算子，将连续的波动方程化解为离散形式。因此，连续的波动方程求解域也要进行离散化，通常是对求解域作网格剖分。图中给出了一维波动方程的时间和空间的网格剖分。规则网格和不规则网格通常在有限差分方法中，空间步长一般取相等距离，这种剖分方法称为规则网格或者均匀网格剖分。在二维弹性波动方程有限差分离散中 ， 两个空间坐标方向的空

11、间步长，通常也取相等的距离。如果 ， 分别取不同的值，则可以称为矩形不规则网格。clearN=4010; dx=0.0024;dt=0.0005; c=dt*dt/dx/dx;x=linspace(0,1,420);u(1:420,1)=0;u(181:240,1)=0.05*sin(pi*x(181:240)*7);u(2:419,2)=u(2:419,1)+c/2*(u(3:420,1)-2*u(2:419,1)+u(1:418,1);h=plot(x,u(:,1),linewidth,3);axis(0,1,-0.05,0.05);set(h,EraseMode,xor,makersiz

12、e,18)for k=2:N set(h,Xdata,x,Ydata,u(:,2); drawnow; u(2:419,3)=2*U(2:419,2)-U(2:419,1)+C*(U(3:420,2). -2*U(2:419,2)+U(1:418,2); u(2:419,1)=u(2:419,2); u(2:419,2)=u(2:419,3);end其中，假设在t0前一时刻波场位移是 ，根据中心差分公式得到具有二阶精度 的差分格式：在上面介绍的初始条件的离散公式中，截断误差是O(t)，误差阶数低于整个差分方程的截断误差 ，误差精度降低，可以通过在时间轴上增加虚拟网格点的方法提高精度。 作业将课

13、中的例题看懂，然后将这个matlab程序改写成fortran程序。有限差分格式相容性、收敛性和稳定性第2章 一维有限差分一维泊松方程及定解问题可描述许多简单的物理现象 (例如)：弹性棒之变形张力作用下弦之变形棒之温度分布第2章 一维有限差分一维泊松方程式解的性质解 总是存在 总是较数据 “平滑”若对于所有 , , 则对于所有 , 已知 则解 为唯一有限差分离散（网格剖分）第2章 一维有限差分将区间 细分成 个相等的子区间，则其中有限差分近似（代替导数）第2章 一维有限差分3、有限差分方程（代数方程）第2章 一维有限差分则 3、有限差分方程（代数方程）第2章 一维有限差分(对称)第2章 一维有限

14、差分3、有限差分方程（代数方程） 可以证明A为对称正定矩阵，线性方程解存在且唯一。对任意因此对任意( 为对称正定)存在且唯一例题：第2章 一维有限差分其中令例题：第2章 一维有限差分有限差分收敛性？第2章 一维有限差分1. 离散解 是否仍保有连续解 之定性性质？2. 当 时解是否更精确？3.对于 是否可使任意变小？第2章 一维有限差分有限差分收敛性？离散稳定性令可证：有限差分收敛性？舍位误差第2章 一维有限差分对任意 可证明令第2章 一维有限差分有限差分收敛性？舍位误差令 为离散误差相减得且第2章 一维有限差分有限差分收敛性？舍位误差第2章 一维有限差分利用离散稳定性估计或有限差分收敛性？舍位

15、误差第2章 一维有限差分例题：第2章 一维有限差分例题：渐近于，第3章 二维有限差分二维泊松方程及定解问题在于第3章 二维有限差分有限差分离散（网格剖分）第3章 二维有限差分有限差分近似（代替导数）第3章 二维有限差分3、有限差分方程（代数方程）令 .例：第3章 二维有限差分第3章 二维有限差分例：具有带状结构，频宽：第3章 二维有限差分方程式编号为 之分量第3章 二维有限差分分块矩阵分块三对角线矩阵第3章 二维有限差分分块定义第3章 二维有限差分3、有限差分方程（代数方程） 可以证明A为对称正定矩阵，线性方程解存在且唯一。对于任何，则( 为 对称正定): 存在且唯一还可证明：第3章 二维有限差分有限差分收敛性