偏微分方程讲义65283.ppt

偏微分方程讲义建模 数值解和Matlab工具箱 温罗生2013 7 提纲 几个偏微分方程的典型问题建模偏微分方程的基本概念偏微分方程的建模偏微分方程的数值方法偏微分方程求解的数值方法偏微分方程的Matlab求解工具Pdetool的使用方法练习题 1 1偏微分方程的基本概念 含有未知多元函数及其偏导数的方程 称之为偏微分方程 方程中出现的未知函数偏导数的最高阶数称为偏微分方程的阶 如果方程中对于未知函数和它的所有偏导数都是线性的 这样的方程称为线性偏微分方程 否则称它为非线性偏微分方程 一些常见的偏微分方程 拉普拉斯方程 uxx uyy uzz 0适用于重力场或静电场 波动方程式 未知函数u x y z t 热传导方程式 其中k代表该材料的热导率 泊松方程 适用于所有物质或电荷的重力场或静电场 初始条件和边界条件称为定解条件 未附加定解条件的偏微分方程称为泛定方程 对于一个具体的问题 定解条件与泛定方程总是同时提出 定解条件与泛定方程作为一个整体 称为定解问题 1 2偏微分方程的建模 例子1 弦的振动方程 一根长为L的均匀细线 线密度为 x 放于x轴上拉紧后 让它离开平衡位置做微小振动 求任意时刻t弦线的位移u x t 使用元素法的思想 在弦上取 x x x 一段进行分析 因为弦的质量相对于张力非常小 因此可以看成弦仅仅受到两个张力 T x 和T x x 利用力学基本定律 可以得到水平和竖直方向两个方程 这里 a分别是两个力和水平方向的夹角 以及弦线在竖直方向的加速度 注意到弦仅仅在接近水平位置振动 所以 和 都是很小的量 于是前一个方程可以近似为 这代表张力和位置无关 直观的理解这是显然的 同时因为 和 都是很小 也有如下的近似关系 于是竖直方向的方程可以改写为 注意到弧长ds近似等于dx 因此上述方程又可以写成 当弦受竖直方向的外力F x t 作用时 竖直方向的方程可以写成 因此方程可以写成 仅仅通过上面的方程并不能完全确定弦的状态 我们必须先给出在初始时刻t 0时的状态 也就是所谓的初始时刻 即 也就是给出初始位置和初始速度 根据方程和上面的初始条件 仍然不能决定 还需要给出边界条件 更加具体的情况 一般有三类边界条件 第一类是两个端点固定 也就是 第二类边界条件描述端点受到外力作用的情况 有 第三类边界条件描述端点受到弹性支撑的情况 结合虎克定理 能得到条件 例子2 热传导方程 热传导定理的推导依据两个定理 一个是傅里叶定理和能量守恒定理 傅里叶定理说热量从高温区域向低温区域传播 热流强度满足 所谓热流强度是单位时间通过单位横截面积的热量 上面dt是时间段 k是热传导系数 是温度场梯度方向面积的投影 u x y z t 是t时刻在空间点 x y z 的温度 取点 x y z 的一个小邻域 从时间t1到t2流入该区域的热量为 因为 该区域的问题升高需要的热量满足 再根据能量守恒定理 传入的热量等于温度升高所需要的热量 于是有 后一方程进一步使用高斯公式 所以有 由t和V的任意性知道 被积函数恒为0 特别的 如果k是一个常量时 得到方程 当物体有内部热源的时候 方程为 因为 左边是升高温度需要的热量 右边第一项是传入的热量 第二项是内部热源的效果 类似与弦振动方程 要确定温度 也需要初始条件和边界条件 初始条件 第二边界条件 表示表面各点热流已知 第一边界条件 表示表面各点温度已知 第三边界条件 表示外界温度为u1 表面的热量和温度差成正比 2 1一些常见的偏微分方程 Poisson方程 带有稳定热源或内部无热源的稳定温度场的温度分布 不可压缩流体的稳定无旋流动及静电场的电势等均满足这类方程 下面的方程是Poisson方程的第一边值问题 第二第三边界条件 在研究热传导过程 气体扩散现象及电磁场的传播等随时间变化的非定常物理问题时 常常会遇到抛物型方程 下面是初边值问题 Cauchy问题 常见的一维振动与波动问题可用该方程表示 下面是初边值问题 双曲型方程 2 2偏微分方程数值解 差分方法又称为有限差分方法或网格法 是求偏微分方程定解问题的数值解中应用最广泛的方法之一 它的基本思想是 先对求解区域作网格剖分 将自变量的连续变化区域用有限离散点 网格点 集代替 将问题中出现的连续变量的函数用定义在网格点上离散变量的函数代替 通过用网格点上函数的差商代替导数 将含连续变量的偏微分方程定解问题化成只含有限个未知数的代数方程组 称为差分格式 如果差分格式有解 且当网格无限变小时其解收敛于原微分方程定解问题的解 则差分格式的解就作为原问题的近似解 数值解 因此 用差分方法求偏微分方程定解问题一般需要解决以下问题 i 选取网格 ii 对微分方程及定解条件选择差分近似 列出差分格式 iii 求解差分格式 椭圆型方程第一边值问题的差分解法 首先将区域通过一系列横纵平行线划分 考察方程在交点处的值 注意到在边界的交点值是已知的 而在内部的交点的值是要求的 为计算u在内部交点处的值 需要列出方程 使用二阶中心差商代替方程中的二阶导数 可以得到方程组 这里h和 是横纵方向上的步长 右端是已知函数f在交点kj处的值 左边分子是待求的未知量 一些内点的四个邻居都属于集合 这种点称为正则点 另一些内点的部分邻居不属于集合 这种点属于非正则点 对于每个正则点 可以得到一个方程 但是非正则点不能得到方程 这样最后得到的方程个数少于未知量个数 因此 对于非正则点要使用一些方法得到方程 常用的方法是直接转移和线性插值的方法 前者就是直接用 内临近的点的值代替 后者是用 内附近点的线性组合代替 下面以例子讲解具体的求解过程 例3 用五点菱形格式求解Laplace方程第一边值问题 首先画出网格 如右图 可以看到 内点一共4个 都是正则点 其他的是边界点 其值已知 使用五点菱形格式 得到方程组如下 代人边界条件的值 得到方程组 借助Matlab 可以轻松得到问题的解 下面是程序 clc clearf1 x 2 log 1 x f2 x log 1 x 2 1 f3 y log 1 y 2 f4 y log 4 y 2 u zeros 4 m 4 n 4 h 1 3 u 1 1 m feval f3 0 h m 1 h u n 1 m feval f4 0 h m 1 h u 1 n 1 feval f1 0 h n 1 h u 1 n m feval f2 0 h n 1 h b u 2 1 u 1 2 u 4 2 u 3 1 u 2 4 u 1 3 u 3 4 u 4 3 a 4110 1 401 10 41 011 4 x a b 当然 这里网格非常稀疏 得到结果比较粗糙 为得到更精确的值 可以将网格加密 比如两个方向各取1000个点 试自己修改上面的程序实现之 functionsol mainf1 x 2 log 1 x f2 x log 1 x 2 1 f3 y log 1 y 2 f4 y log 4 y 2 a 1 b 1 h 1 3 tol 0 00001 max1 1000 sol dirich f1 f2 f3 f4 a b h tol max1 functionU dirich f1 f2 f3 f4 a b h tol max1 Input f1 f2 f3 f4areboundaryfunctionsinputasstrings aandbrightendpointsof 0 a and 0 b hstepsize tolisthetolerance Output Usolutionmatrix Iff1 f2 f3andf4areM filefunctionscallU dirich f1 f2 f3 f4 a b h tol max1 iff1 f2 f3andf4areanonymousfunctionscallU dirich f1 f2 f3 f4 a b h tol max1 InitializeparametersandUn fix a h 1 m fix b h 1 ave a feval f1 0 feval f2 0 b feval f3 0 feval f4 0 2 a 2 b U ave ones n m BoundaryconditionsU 1 1 m feval f3 0 h m 1 h U n 1 m feval f4 0 h m 1 h U 1 n 1 feval f1 0 h n 1 h U 1 n m feval f2 0 h n 1 h 逐次超松弛迭代法 SOR 参数w 4 2 sqrt 4 cos pi n 1 cos pi m 1 2 Refineapproximationsandsweepoperatorthroughoutthegrid 366 err 1 cnt 0 while err tol 上面的程序是使用SOR方法计算线性方程组 因为偏微分方程离散化之后得到的线性方程组比较特殊 同时规模非常大 所以往往选择比较有效的线性方程组数值算法 当问题规模较小 初学时可以不考虑这些 抛物型方程的差分解法 使用差分进行近似 得到如下的方程 这里h和 是横纵方向上的步长 对初始条件及第一类边界条件 可直接得到 结合上页第一个差分格式 得到 计算的时候 由于第1层的数据已知 利用这些数据可以得到第二层的数据 再利用第二层的数据计算第三层 依次类推就能计算出全部的未知量 这个格式称为古典显格式 结合第二个差分格式 得到 这是一个三对角方程组 可以使用追赶法求解 这个格式称为古典隐格式 必须通过解下列线性方程组 例4 用古典显示格式求初边值问题 的数值解 取h 1 0 5 双曲型方程的差分解法 可以使用下面的差分进行近似 上面的格式都是显格式 当然 也可以采用不同的格式组合得到其他的格式 这里不再讲解 2 3偏微分方程的Matlab求解工具 PDEPE命令的使用方法及例子 PDEPE initial boundaryvalueproblemsforparabolic ellipticPDEsin1 D 是求解抛物线和椭圆型一维偏微分方程的命令 为求解这些类型的各种各样的问题 在求解方程的形式 初始条件以及边界条件上 PDEPE使用下面的统一结构 方程形式 函数f用来描述 通量 函数s用来描述 源 初始条件为 边界条件为 可以通过取不同的函数 对应全部的三类边界条件 Matlab的命令格式为pdepe M是参数 可以取0 1 2 pdefun是上面的函数f icfun描述初始条件 bcfun描述边界条件 xmesh和tspan是x和t方向上的划分 options是选择项 下面的两个例子是matlab帮助文件中自带的例子 这里做一个详细的讲解 例子5 求偏微分方程 其中x在 0 1 之间 且满足以下初边值条件 步骤1 将欲求解的偏微分方程改写成如下的标准形式 步骤2 编写偏微分方程的系数向量函数 function c f s ex20 1pdefun x t u dudx c pi 2 f dudx s 0 步骤3 编写起始值条件 functionu0 ex20 1ic x u0 sin pi x 步骤4 编写边界条件 在编写之前 先将边界条件改写成标准形式 找出相对应的函数pa 和pb 然后写出MATLAB的边界条件函数 例如 原边界条件可写成 因而 边界条件函数可编写成function pa qa pb qb ex20 1bc xa ua xb ub t pa ua qa 0 pb pi exp t qb 1 步骤5 取点 例如x linspace 0 1 20 x取20点t linspace 0 2 5 时间取5点输出步骤6 利用pdepe求解 m 0 依步骤1之结果sol pdepe m ex20 1pdefun ex20 1ic ex20 1bc x t 这里sol实际上是二维矩阵步骤7显示结果 surf x t sol title pde数值解 xlabel 位置 ylabel 时间 zlabel u 例6 试解以下联立的偏微分方程系统 步骤1 改写偏微分方程为标准形式 其中m 0和 初始条件 在x 0处 在x 1处 步骤2 编写偏微分方程的系数向量函数 function c f s ex20 2pdefun x t u dudx c 11 f 0 0240 170 dudx y u 1 u 2 F exp 5 73 y exp 11 47 y s FF 步骤3 编写初始条件函数functionu0 ex20 2ic x u0 10 步骤4 编写边界条件函数function pa qa pb qb ex20 2bc xa ua xb ub t pa 0ua 2 qa 10 pb ub 1 10 qb 01 步骤5 取点 m 0 x 00 0050 010 050 10 20 50 70 90 950 990 9951 t 00 0050 010 050 10 511 52 利用pdepe求解 382 sol pdepe m ex20 2pdefun ex20 2ic ex20 2bc x t u1 sol 1 第一个状态之数值解输出u2 sol 2 第二个状态之数值解输出 绘图输出 figure 1 surf x t u1 title u1之数值解 xlabel x ylabel t figure 2 surf x t u2 title u2之数值解 xlabel x ylabel t pdefun函数 function c f s ex20 2pdefun x t u dudx c 11 f 0 0240 170 dudx y u 1 u 2 F exp 5 73 y exp 11 47 y s FF 初始条件函数 functionu0 ex20 2ic x u0 10 边界条件函数 function pa qa pb qb ex20 2bc xa ua xb ub t pa 0ua 2 qa 10 pb ub 1 10 qb 01 例7 求下列偏微分方程组的数值解 将原方程改写成标准式 参数取值 m 1 圆柱 边界条件 在r 0处和r 2处 functionsol main20 6m 1 取点 r linspace 0 2 20 L linspace 0 5 40 利用pdepe求解 sol pdepe m ex20 6pdefun ex20 6ic ex20 6bc r L T sol 1 温度f sol 2 反应率 绘图输出 figure 1 surf L r T title temp xlabel L ylabel r zlabel temp 0C figure 2 surf L r f title reactionrate xlabel L ylabel r zlabel reactionrate 偏微分方程函数 function c f s ex20 6pdefun r L u DuDr T u 1 f u 2 c 11 f 0 590 94 DuDr s 9 611 0 047 T 初始条件函数 functionu0 ex20 6ic x u0 1250 边界条件函数 function pa qa pb qb ex20 6bc ra ua rb ub L pa 00 qa 11 pb 112 ub 1 273 0 qb 0 65 0 591 二维状态空间的偏微分方程的MATLAB解法 MATLAB中的偏微分方程 PDE 工具箱是用有限元法寻求典型偏微分方程的数值近似解 该工具箱求解偏微分方程具体步骤与用有限元方法求解偏微分方程的过程是一致的 包括几个步骤 即几何描述 边界条件描述 偏微分方程类型选择 有限元划分计算网格 初始化条件输入 最后给出偏微分方程的数值解 包括画图 下面我们讨论的方程是定义在平面上的有界区域 上 区域的边界记作 Matlab有下面的一些函数能够求解二维偏微分方程问题 其中最基本的是assempde命令 也能使用pdetool MATLAB工具箱可求解问题的基本类型 i 椭圆型偏微分方程 ii 抛物型偏微分方程 iii 双曲型偏微分方程 iv 特征值问题 v 非线性椭圆偏微分方程 vi 偏微分方程组 边界条件有如下三种 例8 求解泊松方程求解区域为单位圆盘 边界条件为在圆盘边界上u 0 1 问题定义g circleg 单位圆b circleb1 边界上为零条件c 1 a 0 f 1 2 产生初始的三角形网格 p e t initmesh g 3 迭代直至得到误差允许范围内的合格解error err 1 whileerr 0 01 p e t refinemesh g p e t u assempde b p e t c a f 求得数值解exact 1 p 1 2 p 2 2 4 err norm u exact inf error errorerr end 结果显示subplot 2 2 1 pdemesh p e t subplot 2 2 2 pdesurf p t u subplot 2 2 3 pdesurf p t u exact 例9考虑最小表面问题 在圆盘边界上 x u2 这是椭圆型方程 其中 g circleg b circleb2 c 1 sqrt 1 ux 2 uy 2 rtol 1e 3 p e t initmesh g p e t refinemesh g p e t u pdenonlin b p e t c 0 0 Tol rtol pdesurf p t u 例10求解正方形区域 x y 1 x y 1 上的热传导方程 这里是抛物型方程 其中d 1 f 0 a 0 c 1 1 问题定义g squareg 定义正方形区域b squareb1 边界上为零条件c 1 a 0 f 0 d 1 2 产生初始的三角形网格 p e t initmesh g 3 定义初始条件u0 zeros size p 2 1 ix find sqrt p 1 2 p 2 2 0 4 u0 ix 1 4 在时间段为0到0 1的20个点上求解nframe 20 tlist linspace 0 0 1 nframe u1 parabolic u0 tlist b p e t c a f d 5 动画图示结果forj 1 nframepdesurf p t u1 j mv j getframe endmovie mv 10 例11 求解正方形区域 x y 1 x y 1 上的波方程 这里是双曲型方程 其中d 1 f 0 a 0 c 1 1 问题定义g squareg 定义正方形区域b squareb3 定义边界c 1 a 0 f 0 d 1 2 产生初始的三角形网格 p e t initmesh g 3 定义初始条件x p 1 y p 2 u0 atan cos pi x ut0 3 sin pi x exp cos pi y 4 在时间段为0到5的31个点上求解n 31 tlist linspace 0 5 n uu hyperbolic u0 ut0 tlist b p e t c a f d 5 动画图示结果forj 1 npdesurf p t uu j mv j getframe endmovie mv 10 2 4pdetool的使用方法 图形界面解法步骤大致上为 1 定义PDE问题 包括二维空间范围 边界条件以及PDE系数等 2 产生离散化点 并将原PDE方程离散化 3 利用有限元素法 finiteelementmethod FEM 求解并显示答案 前五个按钮为PDE系统的边界范围绘制功能 由左至右的用法为 1 以对角绘制矩形或正方形 按住鼠标左键可绘制矩形 而正方形需以按住右键的方式绘制 2 从中心点至某一角边的方式绘制矩形或正方形 同样地 鼠标左键绘矩形 右键绘正方形 3 由周围界线的方式绘制椭圆或圆形区域 鼠标左键用以绘制椭圆 而右键用来绘制圆形图形 4 以中心点向外的方式绘制椭圆或圆 同样地 鼠标左键及右键 分别用以绘制椭圆及圆形的区域 5 用以绘制多边型等不规则区域 欲关闭此功能需按鼠标右键 在这些绘制按钮之后的7个按钮功能依序如下 6 用以给定边界条件 在此功能选定后 使用者可在任一图形边界上按住鼠标左键双击 然后在对话框中输入边界条件 7 用以指定PDE问题及相关参数 8 产生图形区域内离散化的网点 9 用以进一步将离散化的网点再取密一点 refinemesh 10 在指定PDE系统 边界条件及区域后 按此钮即开始解题 11 用以指定显示结果绘制方式 12 放大缩小功能 便于图形绘制及显示 图形界面解法的使用步骤要利用pdetool接口求解之前 需先定义PDE问题 其包含三大部份 1 利用绘图 draw 模式 定义欲解问题的空间范围 domain 2 利用boundary模式 指定边界条件 3 利用PDE模式 指定PDE系数 即输入c a f和d等PDE模式中的系数 在定义PDE问题之后 可依以下两个步骤求解 1 在mesh模式下 产生mesh点 以便将原问题离散化 2 在solve模式下 求解 3 最后 在Plot模式下 显示答案 以Poisson s方程式 u 10的求解为例 详细说明pdetool的用法 此问题的几何图形及相关边界条件 将于求解过程中加以说明 步骤1 在命令窗口中键入pdetool以进入GUI graphicaluserinterface 界面 选取Options中之Grid功能 以显示网格线 步骤2 利用Draw功能 画出问题之几何图形 请注意 使用者可利用内定对象 多边型 矩形 正方形 圆形 及 椭圆型 予以组合 步骤3 选取PDE功能项 以输入PDE方程的系数及类型 因问题为 u f 故此为椭圆型的问题 且其标准形式为 c u au f 比较得知 c 1 a 0和f 10 所以对话框输入的情况如图3 步骤4 选取Boundary功能 以输入边界条件 假设弧形部分边界条件为Neumann形 且为 u n 5 与标准式比较知 g 5且q 0 直线部分其边界条件则在Dirichlettype使h 0 r 0 对话框输入情况见图4 步骤5 选取Mesh功能 产生网点 使用者亦可进一步利用将网点取得密一点 refinemesh 见图5 步骤6 选取solve功能 解此PDE 见图6 注意 1 MATLAB会以图形的方式展示结果 使用者亦可点选plot下的 parameters 功能 选择适当的方式显示图形及数据 2 另外 若使用者欲将结果输出到命令窗口中 以供后续处理 可利用solve功能项下的 exportsolution 指定变量名称来完成 3 如果求抛物型或双曲型方程的数值解 还需要通过 solve 菜单下的 parameters 选项设置初值条件 4 在上面定义边界条件和初始条件时 可以使用一些内置变量 1 在边界条件输入框中 可以使用如下变量 二维坐标x和y 边界线段长度参数s 外法向矢量的分量nx和ny 如果需要边界的切线方向 可以通过tx ny和ty nx表示 解u 2 在初值条件的输入框中 也可以输入用户定义的MATLAB可接受变量 p e t x y 的函数 例13使用PDETOOL重新求例10的数值解 1 定义PDE问题 包括二维空间范围 边界条件以及PDE系数等 我们这里就省略了 2 区域剖分以后 通过 Mesh 菜单下的 ExportMesh 选项可以把p e