偏微分方程matlab

《偏微分方程matlab》由会员分享，可在线阅读，更多相关《偏微分方程matlab（33页珍藏版）》请在金锄头文库上搜索。

1、 专业资料整理 基础知识偏微分方程的定解问题各种物理性质的定常（即不随时间变化）过程，都可用椭圆型方程来描述。其最典型、最简单的形式是泊松(Poisson)方程 （1）特别地，当 f ( x, y) 0 时，即为拉普拉斯(Laplace)方程，又称为调和方程 （2）带有稳定热源或内部无热源的稳定温度场的温度分布，不可压缩流体的稳定无旋流动及静电场的电势等均满足这类方程。Poisson 方程的第一边值问题为 （3）其 中 为 以 为 边 界 的 有 界区 域 ， 为 分 段 光 滑 曲 线， U 称 为 定 解区 域 ，f (x, y),(x, y) 分别为 , 上的已知连续函数。第二类和第三类

2、边界条件可统一表示成 （4）其中 n 为边界 的外法线方向。当 = 0 时为第二类边界条件， 0 时为第三类边界条件。在研究热传导过程，气体扩散现象及电磁场的传播等随时间变化的非定常物理问题时，常常会遇到抛物型方程。其最简单的形式为一维热传导方程 （5）方程（5）可以有两种不同类型的定解问题：初值问题（也称为 Cauchy 问题） （6）初边值问题 （7）其中 为已知函数，且满足连接条件 完美WORD格式 问题（7）中的边界条件称为第一类界条件。第二类和第三类边界条件为 （8） 其中。当时，为第二类边界条件，否则称为第三类边界条件。双曲型方程的最简单形式为一阶双曲型方程 （9）物理中常见的一维

3、振动与波动问题可用二阶波动方程 （10）描述，它是双曲型方程的典型形式。方程（10）的初值问题为 （11）边界条件一般也有三类，最简单的初边值问题为 如果偏微分方程定解问题的解存在，唯一且连续依赖于定解数据（即出现在方程和定解条件中的已知函数），则此定解问题是适定的。可以证明，上面所举各种定解问题都是适定的。2 偏微分方程的差分解法差分方法又称为有限差分方法或网格法，是求偏微分方程定解问题的数值解中应用最广泛的方法之一。它的基本思想是：先对求解区域作网格剖分，将自变量的连续变化区域用有限离散点（网格点）集代替；将问题中出现的连续变量的函数用定义在网格点上离散变量的函数代替；通过用网格点上函数的

4、差商代替导数，将含连续变量的偏微分方程定解问题化成只含有限个未知数的代数方程组（称为差分格式）。如果差分格式有解，且当网格无限变小时其解收敛于原微分方程定解问题的解，则差分格式的解就作为原问题的近似解（数值解）。因此，用差分方法求偏微分方程定解问题一般需要解决以下问题：（i）选取网格； （ii）对微分方程及定解条件选择差分近似，列出差分格式； （iii）求解差分格式； （iv）讨论差分格式解对于微分方程解的收敛性及误差估计。 下面我们只对偏微分方程的差分解法作一简要的介绍。 2.1 椭圆型方程第一边值问题的差分解法 以 Poisson 方程（1）为基本模型讨论第一边值问题的差分方法。 考虑 P

5、oisson 方程的第一边值问题（3） 取 h, 分别为 x 方向和 y 方向的步长，以两族平行线 将 定 解 区 域 剖 分 成 矩 形 网 格 。 节 点 的 全 体 记 为 为整数 。定解区域内部的节点称为内点，记内点集为 。边界与网格线的交点称为边界点，边界点全体记为 h。与节点 沿 x 方向或 y 方向只差一个步长的点和 称为节点 的相邻节点。如果一个内点的四个相邻节点均属于U ，称为正则内点，正则内点的全体记为(1)，至少有一个相邻节点不属于U 的内点称为非正则内点，非正则内点的全体记为(2)。我们的问题是要求出问题（3）在全体内点上的数值解。 为简便记，记。对正则内点，由二阶中心

6、差商公式Poisson 方程（1）在点处可表示为 （12）在式（12）中略去，即得与方程（1）相近似的差分方程 （13）式（13）中方程的个数等于正则内点的个数，而未知数 , 则除了包含正则内点处解的近似值，还包含一些非正则内点处的近似值，因而方程个数少于未知数个数。在非正则内点处 Poisson 方程的差分近似不能按式（13）给出，需要利用边界条件得到。 边界条件的处理可以有各种方案，下面介绍较简单的两种。 (i) 直接转移 (ii) 线性插值 由式（13）所给出的差分格式称为五点菱形格式，实际计算时经常取h = ，此时 五点菱形格式可化为 （14）简记为 （15）其中 。求解差分方程组最常

7、用的方法是同步迭代法，同步迭代法是最简单的迭代方式。除 边界节点外，区域内节点的初始值是任意取定的。 例 1 用五点菱形格式求解 Laplace 方程第一边值问题 其中。取 。当时，利用点(k, j),(k 1, j .1),(k 1, j +1)构造的差分格式 （16）称为五点矩形格式，简记为 （17）其中。 2.2 抛物型方程的差分解法 以一维热传导方程（5） 为基本模型讨论适用于抛物型方程定解问题的几种差分格式。 首先对xt 平面进行网格剖分。分别取h, 为x方向与t方向的步长，用两族平行直 线 (k = 0,1,2,) ，k ( j = 0,1,2, )，将 xt 平面剖分成矩形网格，

8、节点为 (k = 0,1,2, , j = 0,1,2, )。为简便起见，记 。2.2.1 微分方程的差分近似 在网格内点(k, j)处，对 分别采用向前、向后及中心差商公式，对采用二 阶中心差商公式，一维热传导方程（5）可分别表示为 由此得到一维热传导方程的不同的差分近似 （18） （19） （20）2.2.2 初、边值条件的处理 为用差分方程求解定解问题（6），（7）等，还需对定解条件进行离散化。 对初始条件及第一类边界条件，可直接得到 （21） （22）其中 。对第二、三类边界条件则需用差商近似。下面介绍两种较简单的处理方法。 （i）在左边界(x = 0)处用向前差商近似偏导数 ，在右边

9、界()处用向后差商近似偏导数 ，即即得边界条件（8）的差分近似为 （ii）用中心差商近似，即 则得边界条件的差分近似为 这样处理边界条件，误差的阶数提高了，但式（24）中出现定解区域外的节点(-1, j)和 (n +1, j)，这就需要将解拓展到定解区域外。可以通过用内节点上的u值插值求出 和，也可以假定热传导方程（5）在边界上也成立，将差分方程扩展到边界节点 上，由此消去和 。 2.2.3 几种常用的差分格式 下面我们以热传导方程的初边值问题（7）为例给出几种常用的差分格式。 (i) 古典显式格式 为便于计算，令，式（18）改写成以下形式 将式（18）与（21），（22）结合，我们得到求解问

10、题（7）的一种差分格式 由于第 0 层( j = 0)上节点处的u 值已知，由式（25）即可算出u 在第一层 ( j = 1)上节点处的近似值。重复使用式（25），可以逐层计算出各层节点的近似值。 （ii）古典隐式格式 将（19）整理并与式（21），（22）联立，得差分格式如下其中。虽然第 0 层上的u 值仍为已知，但不能由式（30）直接计算以上各层节点上的值故差分格式（26）称为古典隐式格式。 （iii）杜福特弗兰克尔（DoFortFrankel）格式 DoFortFrankel 格式是三层显式格式，它是由式（24）与（25），（26）结合得到 的。具体形式如下： 用这种格式求解时，除了第

11、0 层上的值由初值条件（21）得到，必须先用二层格式 求出第 1 层上的值，然后再按格式（27）逐层计算 。 2.3 双曲型方程的差分解法 对二阶波动方程（10） 如果令，则方程（10）可化成一阶线性双曲型方程组 记，则方程组（28）可表成矩阵形式 矩阵 A有两个不同的特征值 = a，故存在非奇异矩阵P，使得 作变换，方程组（29）可化成 方程组（30）由两个独立的一阶双曲型方程联立而成。因此下面主要讨论一阶双曲型方 程的差分解法。 考虑一阶双曲型方程的初值问题 与抛物型方程的讨论类似，仍将 xt 平面剖分成矩形网格。取 x 方向步长为 h ，t 方向步 长为 ，网格线为。为简便起 见，记。

12、以不同的差商近似偏导数，可以得到方程（9）的不同的差分近似 结合离散化的初始条件，可以得到几种简单的差分格式。 3 一维状态空间的偏微分方程的 MATLAB 解法 3.1 工具箱命令介绍 MATLAB 提供了一个指令 pdepe，用以解以下的 PDE 方程式 其中时间介于之间，而位置x则介于a,b有限区域之间。m 值表示问题的对称性，其可为 0，1 或 2，分别表示平板(slab)，圆柱(cylindrical)或球体(spherical)的情形。因而，如果m 0，则a必等于b，也就是说其具有圆柱或球体的对称关系。同时，式中一项为流通量(flux)，而为来源(source)项。为偏微分方程的对

13、角线系数矩阵。若某一对角线元素为 0，则表示该偏微分方程为椭圆型偏微分方程，若为正值(不为 0)，则为拋物型偏微分方程。请注意 c 的对角线元素一定不全为 0。偏微分方程初始值可表示为 而边界条件为 其中 x 为两端点位置，即 a 或b 用以解含上述初始值及边界值条件的偏微分方程的 MATLAB 命令 pdepe 的用法如 下： sol = pdepe(m, pdepe,icfun,bcfun, xmesh,tspan,options)其中 m 为问题之对称参数； xmesh为空间变量x的网格点(mesh)位置向量，即 ，其中 。 tspan 为时间变量t的向量，即，其中为起始时间，为 终点时

14、间。 pdefun 为使用者提供的 pde 函数文件。其函数格式如下： c, f , s = pdefun(x,t,u, dudx)亦即，使用者仅需提供偏微分方程中的系数向量。 c ， f 和 s 均为行(column)向量，而向量 c 即为矩阵 c 的对角线元素。 icfun提供解u 的起始值，其格式为u = icfun(x)值得注意的是u为行向量。 bcfun 使用者提供的边界条件函数，格式如下： pl, ql, pr, ql = bcfun(xl,ul, xr,ur,t)其中ul和ur 分别表示左边界(xl = b)和右边界(xr = a) u的近似解。输出变量中，pl 和 ql 分别表示左边界 p 和 q 的行向量，而 pr 和 qr 则为右边界 p 和 q 的行向量。 sol 为解答输出。sol 为多维