m=T 其中n= 对第二、三类边界条件则需用差商近似。下面介绍两种较简单的处理方法， ①在左边界仅=0)处用向前差商近似偏导数，在右边界（仁=D处用向后差 商近似偏导数，即 ou ul.j)-u(0.D+O(h) (U=0,1.,m) x 1 即得边界条件(8)的差分近似为 4-L-40=8 h =0,1.,m) (23) -山+y=82 h (im用中心差商近似 ,即 _LD-D+0) 2h n+-u(n-1D+O() U=0,l.,m) 2h 则得边界条件的差分近似为 4y-u-2y4w=8 uw-+,=82 j=0,l.,m) (24 2h 这样处理边界条件，误差的阶数提高了，但式(24)中出现定解区域外的节点(-1，)和 (n+1,),这就需要将解拓展到定解区域外。可以通过用内节点上的u值插值求出W 和山，也可以假定热传导方程(5)在边界上也成立，将差分方程扩展到边界节点 上，由此消去和un 223几种常用的差分格式 面我们 0古典显式格式 方程的初边值问题(7)为例给出几种常用的差分格式。 为便于计算。令r ,式(18)改写成以下形式 -245
-245- 其中 τ T m h l n = , = 。 对第二、三类边界条件则需用差商近似。下面介绍两种较简单的处理方法。 （i）在左边界(x = 0) 处用向前差商近似偏导数 x u ∂ ∂ ，在右边界(x = l) 处用向后差 商近似偏导数 x u ∂ ∂ ，即 ( ) ( , ) ( 1, ) ( ) (1, ) (0, ) ( , ) (0, ) O h h u n j u n j x u O h h u j u j x u n j j + − − = ∂ ∂ + − = ∂ ∂ ( j = 0,1,L,m) 即得边界条件（8）的差分近似为 ⎪ ⎪ ⎩ ⎪ ⎪ ⎨ ⎧ + = − − = − − j n j j n j n j j j j j j u g h u u u g h u u 2 , 2 , 1, 1 0, 1 1, 0, λ λ ( j = 0,1,L,m) （23） （ii）用中心差商近似 x u ∂ ∂ ，即 ( ) 2 ( 1, ) ( 1, ) ( ) 2 (1, ) ( 1, ) 2 ( , ) 2 (0, ) O h h u n j u n j x u O h h u j u j x u n j j + + − − = ∂ ∂ + − − = ∂ ∂ ( j = 0,1,L,m) 则得边界条件的差分近似为 ⎪ ⎪ ⎩ ⎪ ⎪ ⎨ ⎧ + = − − = − + − − j n j j n j n j j j j j j u g h u u u g h u u 2 , 2 1, 1, 1 0, 1 1, 1, 2 2 λ λ ( j = 0,1,L,m) （24） 这样处理边界条件，误差的阶数提高了，但式（24）中出现定解区域外的节点(−1, j) 和 (n +1, j) ，这就需要将解拓展到定解区域外。可以通过用内节点上的u 值插值求出u−1, j 和un+1, j ，也可以假定热传导方程（5）在边界上也成立，将差分方程扩展到边界节点 上，由此消去u−1, j 和un+1, j 。 2.2.3 几种常用的差分格式 下面我们以热传导方程的初边值问题（7）为例给出几种常用的差分格式。 (i) 古典显式格式 为便于计算，令 2 h a r τ = ，式（18）改写成以下形式
4,=ru+(1-2r)u+ug- 将式(18)与(21)，(22)结合，我们得到求解问题(7)的一种差分格式 [y4=u4y+(0-2r)4y+uk-y(k=1,2,.,n-1j=0,l.,m-) 40=pg (k=1,2,.,n) (25) (j=1,2.,m) 由于第0层(U=0)上节点处的u值已知(u0=p,),由式(25)即可算出u在第一层 U=1)上节点处的近似值·重复使用式(25)，可以逐层计算出各层节点的近似值。 (i)古典隐式格式 将(19)整理并与式(21)，(22)联立，得差分格式如下 [4m=4+r(uk1-21+u4-Hk=12,.,n-1j=01,.,m-l) .=Pu (k=01,.,n) (26) 40=81,山nJ=82) U=0,1,.,m) 其中r=。虽然第0层上的“值仍为已知，但不能由式(30)直接计算以让上各层节 点上的值，故差分格式(26)称为古典隐式格式。 ()杜特一弗卓尔(DoFort -Frankel)格式 -Frankel格式是三层显式格式，它是由式(24)与(25)，(26)结合符到 的。具体形式如下： mi+a=12-2m- 2. 44.0=pg (k=0L,.,) (27) 40y=81,4my=82 U=0,1.,m) 用这种格式求解时，除了第0层上的值山0由初值条件(21)得到，必须先用二层格式 求出第1层上的值1，然后再按格式(27)逐层计算4(=2,3，.，m)。 2.3双曲型方程的差分解法 对二阶波动方程(10) 如果 ，则方程(10)可化成一阶线性双曲型方程组 [=a (28) 记v=(,2)，则方程组(28)可表成矩阵形式 -246-
-246- k j k j k j k j u , 1 ru 1, r u , ru 1, (1 2 ) + = + + − + − 将式（18）与（21），（22）结合，我们得到求解问题（7）的一种差分格式 ⎪ ⎩ ⎪ ⎨ ⎧ = = = = = + = + + − + − = − = − , ( 1,2, , ) ( 1,2, , ) (1 2 ) ( 1,2, , 1, 0,1, , 1) 0, 1 , 2 ,0 , 1 1, , 1, u g u g j m u k n u ru r u ru k n j m j j n j j k k k j k j k j k j L L L L ϕ （25） 由于第 0 层( j = 0) 上节点处的u 值已知( ) uk ,0 = ϕ k ，由式（25）即可算出u 在第一层 ( j = 1) 上节点处的近似值uk ,1 。重复使用式（25），可以逐层计算出各层节点的近似值。 （ii）古典隐式格式 将（19）整理并与式（21），（22）联立，得差分格式如下 ⎪ ⎩ ⎪ ⎨ ⎧ = = = = = + = + + + − + + − + = − = − , ( 0,1, , ) ( 0,1, , ) ( 2 )( 1,2, , 1, 0,1, , 1) 0, 1 , 2 ,0 , 1 , 1, 1 , 1 1, 1 u g u g j m u k n u u r u u u k n j m j j n j j k k k j k j k j k j k j L L L L ϕ （26） 其中 2 h a r τ = 。虽然第 0 层上的u 值仍为已知，但不能由式（30）直接计算以上各层节 点上的值uk , j ，故差分格式（26）称为古典隐式格式。 （iii）杜福特—弗兰克尔（DoFort—Frankel）格式 DoFort—Frankel 格式是三层显式格式，它是由式（24）与（25），（26）结合得到 的。具体形式如下： ⎪ ⎪ ⎩ ⎪ ⎪ ⎨ ⎧ = = = = = = − = − + − + + + + = + − − , ( 0,1, , ) ( 0,1, , ) ( 1,2, , 1, 1,2, , 1) 1 2 1 2 ( ) 1 2 2 0, 1 , 2 ,0 , 1 1, 1, , 1 u g u g j m u k n u k n j m r r u u r r u j j n j j k k k j k j k j k j L L L L ϕ (27) 用这种格式求解时，除了第 0 层上的值uk ,0 由初值条件（21）得到，必须先用二层格式 求出第 1 层上的值uk ,1 ，然后再按格式（27）逐层计算 ( 2,3, , ) uk , j j = L m 。 2.3 双曲型方程的差分解法 对二阶波动方程（10） 2 2 2 2 2 x u a t u ∂ ∂ = ∂ ∂ 如果令 x u v t u v ∂ ∂ = ∂ ∂ 1 = 2 , ，则方程（10）可化成一阶线性双曲型方程组 ⎪ ⎪ ⎩ ⎪ ⎪ ⎨ ⎧ ∂ ∂ = ∂ ∂ ∂ ∂ = ∂ ∂ x v t v x v a t v 2 1 1 2 2 （28） 记 T v (v ,v ) = 1 2 ，则方程组（28）可表成矩阵形式
会]度会 (29) 矩阵A有两个不同的特征值元=土，故存在非奇异矩阵P,使得 PAP-Ta0 Lo-a=A 作变换w=PY=(,w2)》',方程组(29)可化成 (30) 方程组(30)由两个独立的一阶双曲型方程联立而成。因此下面主要讨论一阶双曲型方 程的差分解法。 考虑一阶双曲型方程的初值问题 a+a=0 1>0-0<x<+00 (31) tx,0)=x) -0<X<+0 与抛物型方程的讨论类似，仍将x1平面剖分成矩形网格。取x方向步长为h,1方向步 长为t,网格线为x=x=kk=0,1,2,以，1=1，=jj=0,12,)。为简便起 见，记(k,)=(xk,y),(k,)=(x,y,),P=(x)。 以不同的差商近似偏导数，可以得到方程(9)的不同的差分近似 悲l-%L+a-=0 (32) h -a-e=0 h (33) 4-u-g+一4-=0 (34) 结合离散化的初始条件，可以得到几种简单的差分格式。 S3,】维秋态的微分方程的MAAB解法 MATLAB提供了一个指令pcpe,用以解以下的PDE方程式 0g-”u尝》+0 (35) 其中时间介于1。≤1≤t,之间，而位置x则介于[a,b]有限区域之间。m值表示问题的 对称性，其可 b) 因而，如果 0,则a必 也就是说其具有圆柱或球体 式中f(x,l,“,du/ax)一项为流通量(ux),而s(x,l,山，du/ax)为来源((source)项 c(x,l,u,ux)为偏微分方程的对角线系数矩阵。若某一对角线元素为0，则表示该偏 微分方程为椭圆型偏微分方程，若为正值（不为0），则为抛物型偏微分方程。请注意c的 对角线元素一定不全为0。偏微分方程初始值可表示为 -247
-247- x v A x a v t v ∂ ∂ = ∂ ∂ ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ = ∂ ∂ 1 0 0 2 （29） 矩阵 A 有两个不同的特征值λ = ±a ，故存在非奇异矩阵 P ，使得 ⎥ = Λ ⎦ ⎤ ⎢ ⎣ ⎡ − = − a a PAP 0 0 1 作变换 T w Pv (w ,w ) = = 1 2 ，方程组（29）可化成 x w t w ∂ ∂ = Λ ∂ ∂ （30） 方程组（30）由两个独立的一阶双曲型方程联立而成。因此下面主要讨论一阶双曲型方 程的差分解法。 考虑一阶双曲型方程的初值问题 ⎪ ⎩ ⎪ ⎨ ⎧ = − ∞ < < +∞ = > − ∞ < < +∞ ∂ ∂ + ∂ ∂ u x x x t x x u a t u ( ,0) ( ) 0 0 ϕ (31) 与抛物型方程的讨论类似，仍将 xt 平面剖分成矩形网格。取 x 方向步长为 h ，t 方向步 长为τ ，网格线为 x = x = kh(k = 0,±1,±2,L), k t = t = j ( j = 0,1,2,L) j τ 。为简便起 见，记( , ) ( , ) k j k j = x y ， ( , ) ( , ) k j u k j = u x y ， ( ) k k ϕ = ϕ x 。 以不同的差商近似偏导数，可以得到方程（9）的不同的差分近似 0 , 1 , 1, , = − + + − + h u u a uk j uk j k j k j τ （32） 0 , 1 , , 1, = − − + − − h u u a uk j uk j k j k j τ （33） 0 2 , 1 , 1, 1, = − − + − + − h u u a uk j uk j k j k j τ （34） 结合离散化的初始条件，可以得到几种简单的差分格式。 §3 一维状态空间的偏微分方程的 MATLAB 解法 3.1 工具箱命令介绍 MATLAB 提供了一个指令 pdepe，用以解以下的 PDE 方程式 ( , , , ) ( ( , , , )) ( , , , ) x u s x t u x u x f x t u x x t u x u c x t u m m ∂ ∂ + ∂ ∂ ∂ ∂ = ∂ ∂ ∂ ∂ − （35） 其中时间介于 f t ≤ t ≤ t 0 之间，而位置 x 则介于[a,b]有限区域之间。m 值表示问题的 对称性，其可为 0，1 或 2，分别表示平板(slab)，圆柱(cylindrical)或球体(spherical)的情 形。因而，如果 m > 0 ，则a 必等于b ，也就是说其具有圆柱或球体的对称关系。同时， 式中 f (x,t,u,∂u ∂x) 一项为流通量(flux)，而 s(x,t,u,∂u ∂x) 为来源(source)项。 c(x,t,u,∂u ∂x) 为偏微分方程的对角线系数矩阵。若某一对角线元素为 0，则表示该偏 微分方程为椭圆型偏微分方程，若为正值(不为 0)，则为拋物型偏微分方程。请注意c 的 对角线元素一定不全为 0。偏微分方程初始值可表示为
u(x.1)=v(x) (36) 而边界条件为 x,40+g(x,0f(x,4,0u/ax)=0 (37) 下： sol=pdepe(m,pdepe,icfun,bcfun,xmesh,tspan,options) 其中 m为问题之对称参数 xmesh为空间变量x的网格点(mesh)位置向量，即xmesh=[x,x,xw],其中 x0=a(起点)，w=b(终点). span为时间变量1的向量，即tspan=o,4,.,lw],其中1。为起始时间，1w为 终点时间。 pdefun为使用者提供的pde函数文件。其函数格式如下： 即，使用者仅需提供偏微，可=pdefu(x“,dhu本 方程中的系数向量。c,∫和s均为行(column)向量，而 向量c即为矩阵c的对角线元素。 icim提供解4的起始值，其格式为u=icfx)值得注意的是u为行问量。 bcn使用者提供的边界条件函数，格式如下： [pl.ql.pr.ql]-bcfim(xl.u.xr.u.t) 其中l和w分别表 示左边界(xl=b)和右边界(xr=a)u的近似解。输出变量中，pl 和g分别表示左边界p和q的行向量，而pr和qr则为右边界p和q的行向量。 sol为解答输出。sol为多维的输出向量，sol(,:)为4，的输出，即山，=sol(::,) 元素山，(j,k)=sol(j,k,i)表示在t=tspan()j和x=xmesh(k)时，之答案。 options为求解器的相关解法参数。详细说明参见odeset的使用方法。 注 1.MATLAB PDE求解器DdDe的算法，主要是将原米的椭圆型和地物线型偏微分 方程转化为一组常微分方程。此转换的过程是基于使用者所指定的msh点，以二阶空 间离散化(spatial discretization)技术为之(Keel and Berzins,.l990),然后以odel5s的指令 求解。采用ode15s的ode解法， 主要是因为在离散化的过程中，椭圆型偏微分方程被 转化为 一组代数方程，而抛物线型的偏微 “万程 被专化为 组跃立的分万程· ,变成一组同时伴有微分方程与代数方程的微分代数方程组， 的取 点(me 位置对解的精确度影响很大，君 iti 息时使田 e求器给出 动点取 占 外，若状态“在某些特定点上有较快速的变动时，亦 此处的点取密集些，以增加精确度。值得注意的是pdepe并不会自动做xmesh的自动 点，使用者必须观察解的特性，自行作取点的操作。一般而言，所取的点数至少需大于 3以上。 3.tspan的选取主要是基于使用者对那些特定时间的状态有兴趣而选定。而间距 248-
-248- ( , ) ( ) 0 0 u x t = v x （36） 而边界条件为 p(x,t,u) + q(x,t) f (x,t,u,∂u ∂x) = 0 （37） 其中 x 为两端点位置，即a 或b 用以解含上述初始值及边界值条件的偏微分方程的 MATLAB 命令 pdepe 的用法如 下： sol = pdepe(m, pdepe,icfun,bcfun, xmesh,tspan,options) 其中 m 为问题之对称参数； xmesh为空间变量 x 的网格点(mesh)位置向量，即 [ , , , ] 0 1 N xmesh = x x L x ，其中 x0 = a (起点)， xN = b (终点)。 tspan 为时间变量t 的向量，即 [ , , , ] 0 1 M tspan = t t L t ，其中 0t 为起始时间， Mt 为 终点时间。 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) 为ui 的输出，即u sol(:,:,i) i = 。 元素u ( j, k) sol( j, k,i) i = 表示在t = tspan( j) 和 x = xmesh(k) 时ui 之答案。 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中所指定的输出点位置向量。xmesh=x,x,.,xw]。 即solU,:,)。也就是说其为pdepe输出中第i个输出在各点位置xmesh处， 时间固定为1，=1pan(U)下的解。 0为所欲内插输出点位置向量。此为使用者重新指定的位置向量 1o为基于所指定位置xO,固定时间1，下的相对应输出 duoutdx为相对应的du/d输出值。 ref KeelrD.and M Berzins"A method for the Spatial Discritization of Parabolic Equations in One Space Variable",SIAMJ.Sci.and Sat.Comput.,Vol.11.pp.1-32,1990. 以下将以数个例子，详细说明pdepe的用法。 32 求 界一维偏微分方程 2试解以下之偏微分方程式 at ax? 其中0≤x≤1，且满足以下之条件限制式 ①起始值条件 IC:x,0)=sin(πx) 而边界条件 BC1:0,)=0 BC2:e+0L0=0 注：本间题的解析解为u(xt)=e-r sin) 解下面将叙述求解的步骤与过程。当完成以下各步骤后，可进一步将其汇总为 主程序cx20 状后 步骤1将欲求解的偏微分方程改写成如式的标准式。 此即 cx,l4,au/ax)=π2 f(x.t.u.ou/ox)=ou s(x.1.u,du/dx)=0 和m=0。 步骤2编写偏微分方程的系数向量函数。 function [c,f,s]-ex20 1pdefun (x,t.u,dudx) -249
-249- (step size)的控制由程序自动完成。 4. 若要获得特定位置及时间下的解，可配合以 pdeval 命令。使用格式如下： [ ] uout, duoutdx = pdeval(m, xmesh,ui, xout) 其中 m 代表问题的对称性。m =0 表示平板；m =1 表示圆柱体；m =2 表示球体。其意 义同 pdepe 中的自变量m 。 xmesh为使用者在 pdepe 中所指定的输出点位置向量。 0 1 [, , ] N xmesh x x x = L 。 ui 即 sol( j,:,i) 。也就是说其为 pdepe 输出中第 i 个输出ui 在各点位置 xmesh 处， 时间固定为t tspan( j) j = 下的解。 xout 为所欲内插输出点位置向量。此为使用者重新指定的位置向量。 uout 为基于所指定位置 xout ，固定时间 f t 下的相对应输出。 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 试解以下之偏微分方程式 2 2 2 x u t u ∂ ∂ = ∂ ∂ π 其中0 ≤ x ≤ 1，且满足以下之条件限制式 (i)起始值条件 IC：u(x,0) = sin(π x) (ii)边界条件 BC1：u(0,t) = 0 BC2： 0 (1, ) = ∂ ∂ + − x u t e t π 注：本问题的解析解为u(x,t) e sin( x) t π − = 解 下面将叙述求解的步骤与过程。当完成以下各步骤后，可进一步将其汇总为一 主程序 ex20_1.m，然后求解。 步骤 1 将欲求解的偏微分方程改写成如式的标准式。 0 2 0 0 ⎟ + ⎠ ⎞ ⎜ ⎝ ⎛ ∂ ∂ ∂ ∂ = ∂ ∂ x u x x x t u π 此即 ( ) 2 c x,t,u,∂u ∂x = π ( ) x u f x t u u x ∂ ∂ , , ,∂ ∂ = s( ) x,t,u,∂u ∂x = 0 和m = 0 。 步骤 2 编写偏微分方程的系数向量函数。 function [c,f,s]=ex20_1pdefun(x,t,u,dudx)