T 其中n 对第二、三类边界条件则需用差商近似。下面介绍两种较简单的处理方法。 (i)在左边界(x=0)处用向前差商近似偏导数一,在右边界(x=D)处用向后差 商近似偏导数一,即 3,s-(0,0+O(h) (n,j)-l(n-1,j) +O(h) (n,j) h 即得边界条件(8)的差分近似为 h ,40y= (=0,…,m) g (i)用中心差商近似一,即 l(1,j)-l(-1,j) O(h2) 2h (=0,…,m) l(n+1,j)-l(n-1,) +O(h2) n,) 2h 则得边界条件的差分近似为 l---10y=8 2h (=0,1,…,m) 2h 这样处理边界条件,误差的阶数提高了,但式(24)中出现定解区域外的节点(-1,j)和 (n+1,j),这就需要将解拓展到定解区域外。可以通过用内节点上的a值插值求出u-1 和n1y,也可以假定热传导方程(5)在边界上也成立,将差分方程扩展到边界节点 上,由此消去u1,和n+1 223几种常用的差分格式 下面我们以热传导方程的初边值问题(7)为例给出几种常用的差分格式 (i)古典显式格式 为便于计算,令"/式(18)改写成以下形式
-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）改写成以下形式
lkA=41,+(1-2r)uk,+Plk-1 将式(18)与(21),(22)结合,我们得到求解问题(7)的一种差分格式 =Pk41+(1-2r)uk +ru (k=1,2…,n-1,j=0,1,…,m-1) (k=1,2,…,n) 10o1=g1;ln=8: (=1,2 由于第0层(=0)上节点处的a值已知(u10=k),由式(25)即可算出u在第一层 (=1)上节点处的近似值uk。重复使用式(25),可以逐层计算出各层节点的近似值。 (ⅱi)古典隐式格式 将(19)整理并与式(21),(22)联立,得差分格式如下 e lkj +r(u k-1,j+1 1)(k=1,2,…,n-1,j=0,1,…,m-1) Pu (k=0,…,n) 10=g1 g (=01,…,m) aT Bh。虽然第0层上的a值仍为已知,但不能由式(30)直接计算以上各层节 点上的值uk,,故差分格式(26)称为古典隐式格式。 (ⅲi)杜福特一弗兰克尔( Do Fort- Frankel)格式 DoFort- Frankel格式是三层显式格式,它是由式(24)与(25),(26)结合得到 的。具体形式如下 lyk1(k=1,2,…,n-1,=1,2 1) 1+2 (k=0,1,…,n) uo,=g1, un,j=2j (=0.1,…,m) 用这种格式求解时,除了第0层上的值uk由初值条件(21)得到,必须先用二层格式 求出第1层上的值k1,然后再按格式(27)逐层计算,(=2,3,…,m) 23双曲型方程的差分解法 对二阶波动方程(10) l =a 如果令v1 ,则方程(10)可化成一阶线性双曲型方程组 at- ax (28) 记v=(v1,n2),则方程组(28)可表成矩阵形式
-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）可表成矩阵形式
阵A有两个不同的特征值A=±a,故存在非奇异矩阵P,使得 PAP- A 作变换w=Pv=(,2),方程组(29)可化成 (30) 方程组(30)由两个独立的一阶双曲型方程联立而成。因此下面主要讨论一阶双曲型方 程的差分解法。 考虑一阶双曲型方程的初值问题 0 t>0 (31) l(x,0)=(x) 00<X<+0 与抛物型方程的讨论类似,仍将x平面剖分成矩形网格。取x方向步长为h,t方向步 长为r,网格线为x=xk=Mh(k=0,+1+2…),t=1,=r(=0,2…)。为简便起 见,记(k,=(xk,y),以(k,)=(xk,y),Q=(x) 以不同的差商近似偏导数,可以得到方程(9)的不同的差分近似 l4 llk + a 0 h 结合离散化的初始条件,可以得到几种简单的差分格式 §3一维状态空间的偏微分方程的 MATLAB解法 3.1工具箱命令介绍 MATLAB提供了一个指令 pdepe,用以解以下的PDE方程式 c(x,t 其中时间介于≤t≤t之间,而位置x则介于[a,b]有限区域之间。m值表示问题的 对称性,其可为0,1或2,分别表示平板(slab),圆柱( cylindrical减或球体( spherical的情 形。因而,如果m>0,则a必等于b,也就是说其具有圆柱或球体的对称关系。同时, 式中f(x,1,u,u/ax)一项为流通量(fux),而S(x,t,u,au/ax)为来源( source))项。 c(x,Lu,a/ax)为偏微分方程的对角线系数矩阵。若某一对角线元素为0,则表示该偏 微分方程为椭圆型偏微分方程,若为正值(不为0),则为拋物型偏微分方程。请注意c的 对角线元素一定不全为0。偏微分方程初始值可表示为
-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。偏微分方程初始值可表示为
l(x,0)=v0(x) (36) 而边界条件为 p(x,t,)+q(x,1)f(x,1,u,ou/ax)=0 其中x为两端点位置,即a或b 用以解含上述初始值及边界值条件的偏微分方程的 MATLAB命令 pdepe的用法如 sol= pdepe(m, pdepe, icfun, bcfun, xmesh, tspan, options) 其中 m为问题之对称参数 xmesh为空间变量x的网格点(mesh)位置向量,即 xmesh=[x0,x1…,xx],其中 x0=a(起点),xN=b(终点) Ispan为时间变量t的向量,即spon=[,12…,lM],其中b为起始时间,t为 终点时间。 pdlm为使用者提供的pde函数文件。其函数格式如下 Ic,f, s]=pdefun(x, t, u, dudx) 亦即,使用者仅需提供偏微分方程中的系数向量。c,∫和s均为行( column)向量,而 向量c即为矩阵c的对角线元素。 ic/im提供解l的起始值,其格式为u= icfun(x)值得注意的是u为行向量。 bc/m使用者提供的边界条件函数,格式如下 pl, ql, pr, q]=bcfun(xl,ul, xr, ur;, 1) 其中和r分别表示左边界(xl=b)和右边界(xr=a)u的近似解。输出变量中,pl 和ql分别表示左边界p和q的行向量,而p和q则为右边界P和q的行向量。 sol为解答输出。sol为多维的输出向量,sol(;:i)为u1的输出,即u1=sol(;,) 元素1(,k)=sol(,k,1)表示在t=1spom()和x= xmesh(k)时u1之答案。 options为求解器的相关解法参数。详细说明参见 odeset的使用方法。 1. MATLAB PDE求解器 pdepe的算法,主要是将原来的椭圆型和拋物线型偏微分 方程转化为一组常微分方程。此转换的过程是基于使用者所指定的mesh点,以二阶空 间离散化( spatial discretization)技术为之 Keel and berzins,1990),然后以ode5的指令 求解。采用odel5s的ode解法,主要是因为在离散化的过程中,椭圆型偏微分方程被 转化为一组代数方程,而拋物线型的偏微分方程则被转化为一组联立的微分方程。因而 原偏微分方程被离散化后,变成一组同时伴有微分方程与代数方程的微分代数方程组 故以ode5s便可顺利求解。 2.x的取点(mesh)位置对解的精确度影响很大,若 pdepe求解器给出“…has difficulty finding consistent initial considition”的讯息时,使用者可进一步将mesh点取密 点,即增加mesh点数。另外,若状态u在某些特定点上有较快速的变动时,亦需将 此处的点取密集些,以增加精确度。值得注意的是 pdepe并不会自动做 xmesh的自动取 点,使用者必须观察解的特性,自行作取点的操作。一般而言,所取的点数至少需大于 3以上。 3. tspan的选取主要是基于使用者对那些特定时间的状态有兴趣而选定。而间距
-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.若要获得特定位置及时间下的解,可配合以 deval命令。使用格式如下: Luout, duoutdx]= deval(m, xmesh, ui, xout) 其中 m代表问题的对称性。m=0表示平板;m=1表示圆柱体;m=2表示球体。其意 义同 pdepe中的自变量m。 xme$h为使用者在pepe中所指定的输出点位置向量。 xmesh=[x0,x1…,x] u即sol(;)。也就是说其为 pdepe输出中第i个输出a在各点位置 xmesh处, 时间固定为1=1spun()下的解。 xolt为所欲内插输出点位置向量。此为使用者重新指定的位置向量 为基于所指定位置xout,固定时间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. ll, pp 1-32, 1990 以下将以数个例子,详细说明 pdepe的用法 3.2求解一维偏微分方程 例2试解以下之偏微分方程式 a2u 其中0≤x≤1,且满足以下之条件限制式 (i)起始值条件 IC:(x,0)=sin(丌x) (i)边界条件 BC1:l(0,t)=0 BC2: Te- a(,) =0 注:本问题的解析解为(x,1)=e-sin(m) 解下面将叙述求解的步骤与过程。当完成以下各步骤后,可进一步将其汇总为 主程序ex201m,然后求解 步骤1将欲求解的偏微分方程改写成如式的标准式。 此即 c(x,t, u, u/ax) f(x, t, u, au/ax) s(x,t, u, Ou/ ax)=0 和m=0 步骤2编写偏微分方程的系数向量函数。 function [c, f, s]=ex20 lpdefun(x, t, u, dudx)
-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)