偏微分方程—matlab

基础知识偏微分方程的定解问题各种物理性质的定常（即不随时间变化）过程，都可用椭圆型方程来描述。其最典型、最简单的形式是泊松(Poisson)方程 （1）特别地，当 f ( x, y) 0 时，即为拉普拉斯(Laplace)方程，又称为调和方程 （2）带有稳定热源或内部无热源的稳定温度场的温度分布，不可压缩流体的稳定无旋流动及静电场的电势等均满足这类方程。Poisson 方程的第一边值问题为 （3）其 中 为 以 为 边 界 的 有 界区 域 ， 为 分 段 光 滑 曲 线， U 称 为 定 解区 域 ，f (x, y),(x, y) 分别为 , 上的已知连续函数。第二类和第三类边界条件可统一表示成 （4）其中 n 为边界 的外法线方向。当 = 0 时为第二类边界条件， 0 时为第三类边界条件。在研究热传导过程，气体扩散现象及电磁场的传播等随时间变化的非定常物理问题时，常常会遇到抛物型方程。其最简单的形式为一维热传导方程 （5）方程（5）可以有两种不同类型的定解问题：初值问题（也称为 Cauchy 问题） （6）初边值问题 （7）其中 为已知函数，且满足连接条件问题（7）中的边界条件称为第一类界条件。第二类和第三类边界条件为 （8） 其中。当时，为第二类边界条件，否则称为第三类边界条件。双曲型方程的最简单形式为一阶双曲型方程 （9）物理中常见的一维振动与波动问题可用二阶波动方程 （10）描述，它是双曲型方程的典型形式。方程（10）的初值问题为 （11）边界条件一般也有三类，最简单的初边值问题为 如果偏微分方程定解问题的解存在，唯一且连续依赖于定解数据（即出现在方程和定解条件中的已知函数），则此定解问题是适定的。可以证明，上面所举各种定解问题都是适定的。2 偏微分方程的差分解法差分方法又称为有限差分方法或网格法，是求偏微分方程定解问题的数值解中应用最广泛的方法之一。它的基本思想是：先对求解区域作网格剖分，将自变量的连续变化区域用有限离散点（网格点）集代替；将问题中出现的连续变量的函数用定义在网格点上离散变量的函数代替；通过用网格点上函数的差商代替导数，将含连续变量的偏微分方程定解问题化成只含有限个未知数的代数方程组（称为差分格式）。如果差分格式有解，且当网格无限变小时其解收敛于原微分方程定解问题的解，则差分格式的解就作为原问题的近似解（数值解）。因此，用差分方法求偏微分方程定解问题一般需要解决以下问题：（i）选取网格； （ii）对微分方程及定解条件选择差分近似，列出差分格式； （iii）求解差分格式； （iv）讨论差分格式解对于微分方程解的收敛性及误差估计。 下面我们只对偏微分方程的差分解法作一简要的介绍。 2.1 椭圆型方程第一边值问题的差分解法 以 Poisson 方程（1）为基本模型讨论第一边值问题的差分方法。 考虑 Poisson 方程的第一边值问题（3） 取 h, 分别为 x 方向和 y 方向的步长，以两族平行线 将 定 解 区 域 剖 分 成 矩 形 网 格 。 节 点 的 全 体 记 为 为整数 。定解区域内部的节点称为内点，记内点集为 。边界与网格线的交点称为边界点，边界点全体记为 h。与节点 沿 x 方向或 y 方向只差一个步长的点和 称为节点 的相邻节点。如果一个内点的四个相邻节点均属于U ，称为正则内点，正则内点的全体记为(1)，至少有一个相邻节点不属于U 的内点称为非正则内点，非正则内点的全体记为(2)。我们的问题是要求出问题（3）在全体内点上的数值解。 为简便记，记。对正则内点，由二阶中心差商公式Poisson 方程（1）在点处可表示为 （12）在式（12）中略去，即得与方程（1）相近似的差分方程 （13）式（13）中方程的个数等于正则内点的个数，而未知数 , 则除了包含正则内点处解的近似值，还包含一些非正则内点处的近似值，因而方程个数少于未知数个数。在非正则内点处 Poisson 方程的差分近似不能按式（13）给出，需要利用边界条件得到。 边界条件的处理可以有各种方案，下面介绍较简单的两种。 (i) 直接转移 (ii) 线性插值 由式（13）所给出的差分格式称为五点菱形格式，实际计算时经常取h = ，此时 五点菱形格式可化为 （14）简记为 （15）其中 。求解差分方程组最常用的方法是同步迭代法，同步迭代法是最简单的迭代方式。除 边界节点外，区域内节点的初始值是任意取定的。 例 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 平面剖分成矩形网格，节点为 (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)处用向前差商近似偏导数 ，在右边界()处用向后差商近似偏导数 ，即即得边界条件（8）的差分近似为 （ii）用中心差商近似，即 则得边界条件的差分近似为 这样处理边界条件，误差的阶数提高了，但式（24）中出现定解区域外的节点(-1, j)和 (n +1, j)，这就需要将解拓展到定解区域外。可以通过用内节点上的u值插值求出 和，也可以假定热传导方程（5）在边界上也成立，将差分方程扩展到边界节点 上，由此消去和 。 2.2.3 几种常用的差分格式 下面我们以热传导方程的初边值问题（7）为例给出几种常用的差分格式。 (i) 古典显式格式 为便于计算，令，式（18）改写成以下形式 将式（18）与（21），（22）结合，我们得到求解问题（7）的一种差分格式 由于第 0 层( j = 0)上节点处的u 值已知，由式（25）即可算出u 在第一层 ( j = 1)上节点处的近似值。重复使用式（25），可以逐层计算出各层节点的近似值。 （ii）古典隐式格式 将（19）整理并与式（21），（22）联立，得差分格式如下其中。虽然第 0 层上的u 值仍为已知，但不能由式（30）直接计算以上各层节点上的值故差分格式（26）称为古典隐式格式。 （iii）杜福特弗兰克尔（DoFortFrankel）格式 DoFortFrankel 格式是三层显式格式，它是由式（24）与（25），（26）结合得到 的。具体形式如下： 用这种格式求解时，除了第 0 层上的值由初值条件（21）得到，必须先用二层格式 求出第 1 层上的值，然后再按格式（27）逐层计算 。 2.3 双曲型方程的差分解法 对二阶波动方程（10） 如果令，则方程（10）可化成一阶线性双曲型方程组 记，则方程组（28）可表成矩阵形式 矩阵 A有两个不同的特征值 = a，故存在非奇异矩阵P，使得 作变换，方程组（29）可化成 方程组（30）由两个独立的一阶双曲型方程联立而成。因此下面主要讨论一阶双曲型方 程的差分解法。 考虑一阶双曲型方程的初值问题 与抛物型方程的讨论类似，仍将 xt 平面剖分成矩形网格。取 x 方向步长为 h ，t 方向步 长为 ，网格线为。为简便起 见，记。 以不同的差商近似偏导数，可以得到方程（9）的不同的差分近似 结合离散化的初始条件，可以得到几种简单的差分格式。 3 一维状态空间的偏微分方程的 MATLAB 解法 3.1 工具箱命令介绍 MATLAB 提供了一个指令 pdepe，用以解以下的 PDE 方程式 其中时间介于之间，而位置x则介于a,b有限区域之间。m 值表示问题的对称性，其可为 0，1 或 2，分别表示平板(slab)，圆柱(cylindrical)或球体(spherical)的情形。因而，如果m 0，则a必等于b，也就是说其具有圆柱或球体的对称关系。同时，式中一项为流通量(flux)，而为来源(source)项。为偏微分方程的对角线系数矩阵。若某一对角线元素为 0，则表示该偏微分方程为椭圆型偏微分方程，若为正值(不为 0)，则为拋物型偏微分方程。请注意 c 的对角线元素一定不全为 0。偏微分方程初始值可表示为 而边界条件为 其中 x 为两端点位置，即 a 或b 用以解含上述初始值及边界值条件的偏微分方程的 MATLAB 命令 pdepe 的用法如 下： sol = pdepe(m, pdepe,icfun,bcfun, xmesh,tspan,options)其中 m 为问题之对称参数； xmesh为空间变量x的网格点(mesh)位置向量，即 ，其中 。 tspan 为时间变量t的向量，即，其中为起始时间，为 终点时间。 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 为多维的输出向量，sol(:,: i) 为的输出，即sol(:,:,i)。元素 sol( j, k,i)表示在t = tspan( j)和x = xmesh(k)时之答案。 options 为求解器的相关解法参数。详细说明参见 odeset 的使用方法。 注： 1. MATLAB PDE 求解器 pdepe 的算法，主要是将原来的椭圆型和拋物线型偏微分 方程转化为一组常微分方程。此转换的过程是基于使用者所指定的 mesh 点，以二阶空 间离散化(spatial discretization)技术为之(Keel and Berzins,1990)，然后以 ode15s 的指令 求解。采用 ode15s 的 ode 解法，主要是因为在离散化的过程中，椭圆型偏微分方程被 转化为一组代数方程，而拋物线型的偏微分方程则被转化为一组联立的微分方程。因而， 原偏微分方程被离散化后，变成一组同时伴有微分方程与代数方程的微分代数方程组， 故以 ode15s 便可顺利求解。 2. x 的取点(mesh)位置对解的精确度影响很大，若 pdepe 求解器给出“has difficulty finding consistent initial considition”的讯息时，使用者可进一步将 mesh 点取密 一点，即增加 mesh 点数。另外，若状态u 在某些特定点上有较快速的变动时，亦需将 此处的点取密集些，以增加精确度。值得注意的是 pdepe 并不会自动做 xmesh 的自动取 点，使用者必须观察解的特性，自行作取点的操作。一般而言，所取的点数至少需大于 3 以上。 3. tspan 的选取主要是基于使用者对那些特定时间的状态有兴趣而选定。而间距(step size)的控制由程序自动完成。 4. 若要获得特定位置及时间下的解，可配合以 pdeval 命令。使用格式如下： uout, duoutdx = pdeval(m, xmesh,ui, xout)其中 m 代表问题的对称性。m =0 表示平板；m =1 表示圆柱体；m =2 表示球体。其意义同 pdepe 中的自变量m 。 xmesh 为使用者在 pdepe 中所指定的输出点位置向量。 即sol( j,:,i)。也就是说其为 pdepe 输出中第 i 个输出在各点位置xmesh处，时间固定为下的解。 xout 为所欲内插输出点位置向量。此为使用者重新指定的位置向量。 uout 为基于所指定位置 xout ，固定时间下的相对应输出。 duoutdx 为相对应的du/dx 输出值。 ref. Keel,R.D. and M. Berzins,“A Method for the Spatial Discritization of Parabolic Equations in One Space Variable”，SIAM J. Sci. and Sat. Comput.,Vol.11,pp.1-32,1990. 以下将以数个例子，详细说明 pdepe 的用法。 3.2 求解一维偏微分方程 例 2 试解以下之偏微分方程式 其中0 x 1，且满足以下之条件限制式 (i)起始值条件 IC：u(x,0) = sin( x) (ii)边界条件 注：本问题的解析解为 解 下面将叙述求解的步骤与过程。当完成以下各步骤后，可进一步将其汇总为一主程序 ex20_1.m，然后求解。 步骤 1 将欲求解的偏微分方程改写成如式的标准式。此即 和m = 0。步骤 2 编写偏微分方程的系数向量函数。 function c,f,s=ex20_1pdefun(x,t,u,dudx) c=pi2; f=dudx; s=0; 步骤 3 编写起始值条件。 function u0=ex20_1ic(x) u0=sin(pi*x); 步骤 4 编写边界条件。在编写之前，先将边界条件改写成标准形式，如式（37）， 找出相对应的 p(.)和q(.)函数，然后写出 MATLAB 的边界条件函数，例如，原边界条 件可写成 即 pl = u(0,t), ql = 0,和 因而，边界条件函数可编写成 function pl,ql,pr,qr=ex20_1bc(xl,ul,xr,ur,t) pl=ul; ql=0; pr=pi*exp(-t); qr=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); 步骤 7 显示结果。 u=sol(:,:,1); surf(x,t,u) title(pde 数值解) xlabel(位置) ylabel(时间 ) zlabel(u) 若要显示特定点上的解，可进一步指定 x 或 t 的位置，以便绘图。例如，欲了解时 间为 2(终点)时，各位置下的解，可输入以下指令(利用 pdeval 指令)： figure(2); %绘成图 2 M=length(t); %取终点时间的下标 xout=linspace(0,1,100); %输出点位置 uout,dudx=pdeval(m,x,u(M,:),xout); plot(xout,uout); %绘图 title(时间为 2 时,各位置下的解) xlabel(x) ylabel(u) 综合以上各步骤，可写成一个程序求解例 2。其参考程序如下： function ex20_1 %* %求解一维热传导偏微分方程的一个综合函数程序 %* m=0; x=linspace(0,1,20); %xmesh t=linspace(0,2,20); %tspan %* %以 pde 求解 %* sol=pdepe(m,ex20_1pdefun,ex20_1ic,ex20_1bc,x,t); u=sol(:,:,1); %取出答案 %* %绘图输出 %* figure(1) surf(x,t,u) title(pde 数值解) xlabel(位置 x) ylabel(时间 t ) zlabel(数值解 u) %* %与解析解做比较 %* figure(2) surf(x,t,exp(-t)*sin(pi*x); title(解析解) xlabel(位置 x) ylabel(时间 t ) zlabel(数值解 u) %* %t=tf=2 时各位置之解 %* figure(3) M=length(t); %取终点时间的下表 xout=linspace(0,1,100); %输出点位置 uout,dudx=pdeval(m,x,u(M,:),xout); plot(xout,uout); %绘图 title(时间为 2 时,各位置下的解) xlabel(x) ylabel(u) %* %pde 函数 %* function c,f,s=ex20_1pdefun(x,t,u,dudx) c=pi2; f=dudx; s=0; %* %初始条件函数 %* function u0=ex20_1ic(x) u0=sin(pi*x); %* %边界条件函数 %* function pl,ql,pr,qr=ex20_1bc(xl,ul,xr,ur,t) pl=ul; ql=0; pr=pi*exp(-t); qr=1; 例 3 试解以下联立的偏微分方程系统 其中且0 x 1和t 0。此联立偏微分方程系统满足以下初边值条件。 (i)初值条件 (ii)边值条件 解 步骤 1：改写偏微分方程为标准式 因此 和m = 0。另外，左边界条件( x = 0处)。写成 即同理，右边界条件( x =1处)为 即步骤 2：编写偏微分方程的系数向量函数。 function c,f,s=ex20_2pdefun(x,t,u,dudx) c=1 1; f=0.024 0.170.*dudx; y=u(1)-u(2); F=exp(5.73*y)-exp(-11.47*y); s=-F F; 步骤 3：编写初始条件函数 function u0=ex20_2ic(x) u0=1 0; 步骤 4：编写边界条件函数 function pl,ql,pr,qr=ex20_2bc(xl,ul,xr,ur,t) pl=0 ul(2); ql=1 0; pr=ur(1)-1 0; qr=0 1; 步骤 5： 取点。 由于此问题的端点均受边界条件的限制，且时间 t 很小时状态的变动很大(由多次求解后的经验得知)，故在两端点处的点可稍微密集些。同时对于 t 小处亦可取密一些。例如， x=0 0.005 0.01 0.05 0.1 0.2 0.5 0.7 0.9 0.95 0.99 0.995 1; t=0 0.005 0.01 0.05 0.1 0.5 1 1.5 2; 以上几个主要步骤编写完成后，事实上就可直接完成主程序来求解。此问题的参考程序如下： function ex20_2 %*%求解一维偏微分方程组的一个综合函数程序 %* m=0; x=0 0.005 0.01 0.05 0.1 0.2 0.5 0.7 0.9 0.95 0.99 0.995 1; t=0 0.005 0.01 0.05 0.1 0.5 1 1.5 2; %* %利用 pdepe 求解 %* 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) %* %pde 函数 %* function c,f,s=ex20_2pdefun(x,t,u,dudx) c=1 1; f=0.024 0.170.*dudx; y=u(1)-u(2); F=exp(5.73*y)-exp(-11.47*y); s=-F F; %* %初始条件函数 %* function u0=ex20_2ic(x) u0=1 0; %* %边界条件函数 %* function pl,ql,pr,qr=ex20_2bc(xl,ul,xr,ur,t) pl=0 ul(2); ql=1 0; pr=ur(1)-1 0; qr=0 1; 3.3 化工应用实例 例 4 触煤反应装置内温度及转换率的分布 以外部热交换式的管形固定层触煤反应装置，进行苯加氢反应产生环己烷。此反应 系统之质量平衡及热平衡方程式如下： 其中T 为温度()， f 为反应率，L 为轴向距离，r 为径向距离。此系统的边界条件为 此外，式中之相关数据及操作条件如下： (i)反应速率式 其中 P 表示分压(atm)，而速率参数为 上式中，下标 B，H 及 C 分别代表苯，氢及环己烷。R 为理想气体常数(1.987cal/molK)。(ii)操作条件及物性数据 题意解析： 因反应速率式 A r 与分压有关，而分压又与反应率 f 有关。故需进一步将 A r 由反应 率 f 表示，方能求解偏微分方程。基于以下的反应方程 则各分压与总压之关系为 将上式，连同反应速率式，带入平衡方程式中，配合边界条件，可利用 pdepe 求解。 MATLAB 程序设计 将原方程改写成如式（35）的标准式 因此 和m = +1(圆柱)。另外，左边界条件( r = 0处)写成 即 同理右边界条件()可写成 即 根据以上的分析，可编写 MATLAB 程序求解此 PDE 问题，其参考程序如下： function ex60_3_1 %* % 触媒反应器内温度及转化率的分布 %* global Pt rw Tw G M y0 Mav rho_B Cp dHr h0 u R ke hw De %* % 给定数据 %* Pt=1.25; %总压(atm) rw=0.025; %管径(m) Tw=100+273; %壁温() G=631; %质量流率(kg/m2hr) M=30; y0=0.0323; Mav=4.47; rho_B=1200; Cp=1.74; dHr=-49250; h0=65.8; T0=125+273; Lw=1; u=8.03; R=1.987; ke=0.65; hw=112; De=0.755; %* m=1; %* % 取点 %* r=linspace(0,rw,10); L=linspace(0,Lw,10); %* % 利用 pdepe 求解 %* sol=pdepe(m,ex20_3_1pdefun,ex20_3_1ic,ex20_3_1bc,r,L); T=sol(:,:,1); %温度 f=sol(:,:,2); %反应率 %* % 绘图输出 %* figure(1) surf(L,r,T-273) title(temp) xlabel(L) ylabel(r) zlabel(temp (0C) % figure(2) surf(L,r,f) title(reaction rate) xlabel(L) ylabel(r) zlabel(reaction rate) %* % PDE 函数 %* function c1,f1,s1=ex20_3_1pdefun(r,L,u1,DuDr) global Pt rw Tw G M y0 Mav rho_B Cp dHr h0 u R ke hw De T=u1(1); f=u1(2); % k=exp(-12100/(R*T)+32.3/R); Kh=exp(15500/(R*T)-31.9/R); Kb=exp(11200/(R*T)-23.1/R); Kc=exp(8900/(R*T)-19.4/R); % a=1+M-3*f; ph=Pt*(M-3*f)/a; pb=Pt*(1-f)/a; pc=Pt*f/a; % rA=k*Kh3*Kb*ph3*pb/(1+Kh*ph+Kb*pb+Kc*pc)4; % c1=1 1; f1=ke/(G*Cp) De/u.*DuDr; %s1=ke/(G*Cp*r)*DuDr(1)-rA*rho_B*dHr/(G*Cp)-2*h0*(T-Tw)/(rw) s1=-rA*rho_B*dHr/(G*Cp);rA*rho_B*Mav/(G*y0); %* %初始条件函数 %* function u0=ex20_3_1ic(x) u0=125+273 0; %* % 边界条件%* function pl,ql,pr,qr=ex20_3_1bc(rl,ul,rr,ur,L) global Pt rw Tw G M y0 Mav rho_B Cp dHr h0 u R ke hw De pl=0 0; ql=1 1; pr=hw*(ur(1)-Tw) 0; qr=G*Cp 1; 例 5 扩散系统之浓度分布 参考如图 3 的装置。管中储放静止液体 B，高度为 L=10 ，放置于充满 A 气体的环境中。假设与 B 液体接触面之浓度为，且此浓度不随时间改变而改变，即在操作时间内( h = 10天)维持定值。气体 A 在液体 B 中之扩散系数为。试决定以下两种情况下，气体 A溶于液体 B中之流通量(flux)。 (a) A 与 B 不发生反应； (b) A 与 B 发生以下之反应 A+ BC ， 其反应速率常数。 题意解析： (a) 因气体 A 与液体 B 不发生反应，故其扩散现象的质量平衡方程如下： 依题意，其初始及边界条件为 (b) 在气体 A 与液体 B 会发生一次反应的情况下，其质量平衡方程需改写为而起始及边界条件同上。 在获得浓度分布后，即可以 Ficks law计算流通量。 MATLAB 程序设计： 此问题依旧可以利用 pdepe 迅速求解。现就各状况的处理过程简述如下。 (a)与标准式（35）比较，可得C = 1，s = 0，和m = 0。另外，经与式（37）比较后得知，左边界及右边界条件之系数分别为 左边界( z = 0 )： 右边界( z = L )： (b)与标准式（35）比较，可得m = 0，C = 1，和 。而边界条件之系数同(a)之结果。 利用以上的处理结果，可编写 MATLAB 参考程序如下： function ex20_3_2 %* % 扩散系统之浓度分布 %* clear clc global DAB k CA0 %* % 给定数据 %* CA0=0.01; L=0.1; DAB=2e-9; k=2e-7; h=10*24*3600; %* % 取点 %* t=linspace(0,h,100); z=linspace(0,L,10); %* % case (a) %* m=0; sol=pdepe(m,ex20_3_2pdefuna,ex20_3_2ic,ex20_3_2bc,z,t); CA=sol(:,:,1); for i=1:length(t) CA_i,dCAdz_i=pdeval(m,z,CA(i,:),0); NAz(i)=-dCAdz_i*DAB; end figure(1) subplot(211) surf(z,t/(24*3600),CA) title(case (a) xlabel(length (m) ylabel(time (day) zlabel(conc. (mol/m3) subplot(212) plot(t/(24*3600),NAz*24*3600) xlabel(time (day) ylabel(flux (mol/m2.day) %* % case (b) %* m=0; sol=pdepe(m,ex20_3_2pdefunb,ex20_3_2ic,ex20_3_2bc,z,t); CA=sol(:,:,1); for i=1:length(t) CA_i,dCAdz_i=pdeval(m,z,CA(i,:),0); NAz(i)=-dCAdz_i*DAB; end % figure(2) subplot(211) surf(z,t/(24*3600),CA) title(case (b) xlabel(length (m) ylabel(time (day) zlabel(conc. (mol/m3) subplot(212) plot(t/(24*3600),NAz*24*3600) xlabel(time (day) ylabel(flux (mol/m2.day) %* % PDE 函数 %* % case (a) %* function c,f,s=ex20_3_2pdefuna(z,t,CA,dCAdz) global DAB k CA0 c=1; f=DAB*dCAdz; s=0; %* % case (a) %* function c,f,s=ex20_3_2pdefunb(z,t,CA,dCAdz) global DAB k CA0 c=1; f=DAB*dCAdz; s=k*CA; %* % 初始条件函数 %* function CA_i=ex20_3_2ic(z) CA_i=0; %* % 边界条件函数 %* function pl,ql,pr,qr=ex20_3_2bc(zl,CAl,zr,CAr,t) global DAB k CA0 pl=CAl-CA0; ql=0; pr=0; qr=1/DAB; 4 二维状态空间的偏微分方程的 MATLAB 解法 MATLAB 中的偏微分方程（PDE）工具箱是用有限元法寻求典型偏微分方程的数值近似解，该工具箱求解偏微分方程具体步骤与用有限元方法求解偏微分方程的过程是一致的，包括几个步骤，即几何描述、边界条件描述、偏微分方程类型选择、有限元划分计算网格、初始化条件输入，最后给出偏微分方程的数值解（包括画图）。 下面我们讨论的方程是定义在平面上的有界区域上，区域的边界记作。 4.1 方程类型 MATLAB 工具箱可以解决下列类型的偏微分方程： （i）椭圆型偏微分方程 其中c, a, f 和未知的u可以是上的复值函数。 （ii）抛物型偏微分方程 其中c,a, f ,d 可以依赖于时间t。 （iii）双曲型偏微分方程 （iv）特征值问题 其中 是未知的特征值， d 是 上的复值函数。 （v）非线性椭圆偏微分方程 其中c, a, f 可以是u的函数。 （vi）方程组 4.2 边界条件 边界条件有如下三种： （i）Dirichlet 条件：hu = r on 。 （ii）Neumann 条件： () + qu = g on这里为区域的单位外法线，h, r, q, g是定义在上的复值函数。 对于二维方程组情形，Dirichlet 边界条件为 Neumann 边界条件为： （iii）对于偏微分方程组，混合边界条件为 这里 的计算是使得满足 Dirichlet 边界条件。 4.3 求解偏微分方程 例 6 求解泊松方程 ，求解区域为单位圆盘，边界条件为在圆盘边界上u = 0。 解 它的精确解为 下面求它的数值解，编写程序如下： %(1)问题定义 g=circleg; %单位圆 b=circleb1; %边界上为零条件 c=1;a=0;f=1; %（2）产生初始的三角形网格 p,e,t=initmesh(g); %（3）迭代直至得到误差允许范围内的合格解 error=; err=1; while err 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=error err; 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) 例7 考虑最小表面问题 在圆盘边界上u = x2 。 解 这是椭圆型方程，其中，编写程序如下：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) 例8 求解正方形区域(x, y) | x, y 1上的热传导方程初始条件为 边界条件为Dirichlet条件u = 0。 解 这里是抛物型方程，其中c = 1, a = 0, f = 0, d = 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)动画图示结果 for j=1:nframe pdesurf(p,t,u1(:,j); mv(j)=getframe; end movie(mv,10) 例9 求解正方形区域(x, y) | x, y 1上的波方程 初始条件为u(0) = arctan(cos(x)， 3sin(x) exp(cos(y )，边界条件为x =1上满足Dirichlet条件u = 0，在 y = 1上满足Neumann条件 = 0 解 这里是双曲型方程，其中c = 1, a = 0, f = 0, d = 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)动画图示结果 for j=1:n pdesurf(p,t,uu(:,j); mv(j)=getframe; end movie(mv,10) 例 10 求解泊松方程 求解区域为单位圆盘，边界条件为在圆盘边界上u = 0。 解 它的精确解为 下面求它的数值解，编写程序如下： g=circleg; b=circleb1; c=1;a=0;f=circlef; p,e,t=initmesh(g); p,e,t=refinemesh(g,p,e,t); u=assempde(b,p,e,t,c,a,f); exact=-1/(2*pi)*log(sqrt(p(1,:).2+p(2,:).2); 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) 4.4 偏微分方程的 pdetool 解法 4.4.1 图形界面解法简介 对于一般的区域，任意边界条件的偏微分方程，我们可以利用 MATLAB 中 pdetool提供的偏微分方程用户图形界面解法。 图形界面解法步骤大致上为： （1）定义 PDE 问题，包括二维空间范围，边界条件以及 PDE 系数等。 （2）产生离散化之点，并将原 PDE 方程式离散化。 （3）利用有限元素法(finite element method；FEM)求解并显示答案。 在说明此解法工具之前，先介绍此 PDE 图形界面的菜单下方的功能图标(icon)按钮。透过这些按钮，使用者可轻松地完成偏微分方程的求解。现将这些按钮的主要功能叙述如下： 前五个按钮为 PDE 系统之边界范围绘制功能，由左至右之用法为： ：以对角绘制矩形或正方形。按住鼠标左键可绘制矩形，而正方形需以按住 右键的方式绘制。 ：从中心点至某一角边的方式绘制矩形或正方形。同样地，鼠标左键绘矩形， 右键绘正方形。 ：由周围界线的方式绘制椭圆或圆形区域。鼠标左键用以绘制椭圆，而右键 用来绘制圆形图形。 ：以中心点向外的方式绘制椭圆或圆。同样地，鼠标左及右键，分别用以绘 制椭圆及圆形的区域。 ：用以绘制多边型等不规则区域，欲关闭此功能需按鼠标右键。 在这些绘制按钮之后的按钮功能依序如下： ：用以给定边界条件。在此功能选定后，使用者可在任一图形边界上按住鼠 标左键双击，然后在对话框中输入边界条件。 ：用以指定 PDE 问题及相关参数。 ：产生图形区域内离散化的网点。 ：用以进一步将离散化的网点再取密一点(refine mesh)。 ：在指定 PDE 系统，边界条件及区域后，按此钮即开始解题。 ：用以指定显示结果绘制方式。