第20章 偏微分方程的数值解

-240- 第二十章 偏微分方程的数值解 自然科学与 工程 路基工程安全技术交底工程项目施工成本控制工程量增项单年度零星工程技术标正投影法基本原理 技术中种种运动发展过程与平衡现象各自遵守一定的规律。这些规 律的定量表述一般地呈现为关于含有未知函数及其导数的方程。我们将只含有未知多元 函数及其偏导数的方程，称之为偏微分方程。 方程中出现的未知函数偏导数的最高阶数称为偏微分方程的阶。如果方程中对于未 知函数和它的所有偏导数都是线性的，这样的方程称为线性偏微分方程，否则称它为非 线性偏微分方程。 初始条件和边界条件称为定解条件，未附加定解条件的偏微分方程称为泛定方程。 对于一个具体的问题，定解条件与泛定方程总是同时提出。定解条件与泛定方程作为一 个整体，称为定解问题。 §1 偏微分方程的定解问题 各种物理性质的定常（即不随时间变化）过程，都可用椭圆型方程来描述。其最典 型、最简单的形式是泊松(Poisson)方程 ),(2 2 2 2 yxf y u x uu =∂ ∂+∂ ∂=Δ （1） 特别地，当 0),( ≡yxf 时，即为拉普拉斯(Laplace)方程，又称为调和方程 02 2 2 2 =∂ ∂+∂ ∂=Δ y u x uu （2） 带有稳定热源或内部无热源的稳定温度场的温度分布，不可压缩流体的稳定无旋流动及 静电场的电势等均满足这类方程。 Poisson 方程的第一边值问题为 ⎪⎩ ⎪⎨ ⎧ Ω∂=Γ= Ω∈=∂ ∂+∂ ∂ Γ∈ ),(|),( ),(),( ),( 2 2 2 2 yxyxu yxyxf y u x u yx ϕ （3） 其中Ω 为以 Γ 为边界的有界区域， Γ 为分段光滑曲线， ΓΩU 称为定解区域， ),(),,( yxyxf ϕ 分别为 ΓΩ, 上的已知连续函数。 第二类和第三类边界条件可统一表示成 ),(),( yxun u yx ϕα =⎟⎠ ⎞⎜⎝ ⎛ +∂ ∂ Γ∈ （4） 其中 n为边界Γ的外法线方向。当 0=α 时为第二类边界条件， 0≠α 时为第三类边界 条件。 在研究热传导过程，气体扩散现象及电磁场的传播等随时间变化的非定常物理问 题时，常常会遇到抛物型方程。其最简单的形式为一维热传导方程 )0(02 2 >=∂ ∂−∂ ∂ a x ua t u （5） 方程（5）可以有两种不同类型的定解问题： 初值问题（也称为 Cauchy 问题） -241- ⎪⎩ ⎪⎨ ⎧ +∞<<∞−= +∞<<∞−>=∂ ∂−∂ ∂ xxxu xt x ua t u )()0,( ,002 2 ϕ （6） 初边值问题 ⎪⎪ ⎪ ⎩ ⎪⎪ ⎪ ⎨ ⎧ ≤≤== = <<<<=∂ ∂−∂ ∂ Tttgtlutgtu xxu lxTt x ua t u 0),(),(),(),0( )()0,( 0,00 21 2 2 ϕ （7） 其中 )(),(),( 21 tgtgxϕ 为已知函数，且满足连接条件 )0()(),0()0( 21 glg == ϕϕ 问题（7）中的边界条件 )(),(),(),0( 21 tgtlutgtu == 称为第一类边界条件。第二类和 第三类边界条件为 Tttgut x u Tttgut x u lx x ≤≤=⎥⎦ ⎤⎢⎣ ⎡ +∂ ∂ ≤≤=⎥⎦ ⎤⎢⎣ ⎡ −∂ ∂ = = 0),()( 0),()( 22 1 0 1 λ λ （8） 其中 0)(,0)( 21 ≥≥ tt λλ 。当 0)()( 21 ≡= tt λλ 时，为第二类边界条件，否则称为第三 类边界条件。 双曲型方程的最简单形式为一阶双曲型方程 0=∂ ∂+∂ ∂ x ua t u （9） 物理中常见的一维振动与波动问题可用二阶波动方程 2 2 2 2 2 x ua t u ∂ ∂=∂ ∂ （10） 描述，它是双曲型方程的典型形式。方程（10）的初值问题为 ⎪⎪ ⎪ ⎩ ⎪⎪ ⎪ ⎨ ⎧ +∞<<∞−=∂ ∂ +∞<<∞−= +∞<<∞−>∂ ∂=∂ ∂ = xx t u xxxu xt x ua t u t )( )()0,( ,0 0 2 2 2 2 2 φ ϕ （11） 边界条件一般也有三类，最简单的初边值问题为 -242- ⎪⎪ ⎪ ⎩ ⎪⎪ ⎪ ⎨ ⎧ ≤≤== ≤≤=∂ ∂= <<>∂ ∂=∂ ∂ = Tttgtlutgtu lxx t uxxu lxt x ua t u t 0)(),(),(),0( 0)(),()0,( 0,0 21 0 2 2 2 2 2 φϕ 如果偏微分方程定解问题的解存在，唯一且连续依赖于定解数据（即出现在方程 和定解条件中的已知函数），则此定解问题是适定的。可以 证明 住所证明下载场所使用证明下载诊断证明下载住所证明下载爱问住所证明下载爱问 ，上面所举各种定解问 题都是适定的。 §2 偏微分方程的差分解法 差分方法又称为有限差分方法或网格法，是求偏微分方程定解问题的数值解中应用 最广泛的方法之一。它的基本思想是：先对求解区域作网格剖分，将自变量的连续变化 区域用有限离散点（网格点）集代替；将问题中出现的连续变量的函数用定义在网格点 上离散变量的函数代替；通过用网格点上函数的差商代替导数，将含连续变量的偏微分 方程定解问题化成只含有限个未知数的代数方程组（称为差分格式）。如果差分格式有 解，且当网格无限变小时其解收敛于原微分方程定解问题的解，则差分格式的解就作为 原问题的近似解（数值解）。因此，用差分方法求偏微分方程定解问题一般需要解决以 下问题： （i）选取网格； （ii）对微分方程及定解条件选择差分近似，列出差分格式； （iii）求解差分格式； （iv）讨论差分格式解对于微分方程解的收敛性及误差估计。 下面我们只对偏微分方程的差分解法作一简要的介绍。 2.1 椭圆型方程第一边值问题的差分解法 以 Poisson 方程（1）为基本模型讨论第一边值问题的差分方法。 考虑 Poisson 方程的第一边值问题（3） ⎪⎩ ⎪⎨ ⎧ Ω∂=Γ= Ω∈=∂ ∂+∂ ∂ Γ∈ ),(|),( ),(),( ),( 2 2 2 2 yxyxu yxyxf y u x u yx ϕ 取 τ,h 分别为 x 方向和 y 方向的步长，以两族平行线 τjyykhxx jk ==== , ),2,1,0,( L±±=jk 将 定 解 区 域 剖 分 成 矩 形 网 格 。 节 点 的 全 体 记 为 },,,|),{( 为整数jijykhxyxR jkjk τ=== 。定解区域内部的节点称为内点，记内点 集 ΩIR 为 τhΩ 。边界Γ与网格线的交点称为边界点，边界点全体记为 τhΓ 。与节点 ),( jk yx 沿 x 方向或 y 方向只差一个步长的点 ),( 1 jk yx ± 和 ),( 1±jk yx 称为节点 ),( jk yx 的相邻节点。如果一个内点的四个相邻节点均属于 ΓΩU ，称为正则内点，正 则内点的全体记为 )1(Ω ，至少有一个相邻节点不属于 ΓΩU 的内点称为非正则内点， 非正则内点的全体记为 )2(Ω 。我们的问题是要求出问题（3）在全体内点上的数值解。 为简便记，记 ),(),,(),(),,(),( , jkjkjkjk yxffyxujkuyxjk === 。对正则内点 -243- )1(),( Ω∈jk ，由二阶中心差商公式 )(),1(),(2),1( 22 ),( 2 2 hO h jkujkujku x u jk +−+−+=∂ ∂ )()1,(),(2)1,( 22 ),( 2 2 ττ O jkujkujku y u jk +−+−+=∂ ∂ Poisson 方程（1）在点 ),( jk 处可表示为 )( )1,(),(2)1,(),1(),(2),1( 22 , 22 τ τ ++= −+−++−+−+ hOf jkujkujku h jkujkujku jk （12） 在式（12）中略去 )( 22 τ+hO ，即得与方程（1）相近似的差分方程 jk jkjkjkjkjkjk f uuu h uuu ,2 1,,1, 2 ,1,,1 22 =+−++− −+−+ τ (13) 式（13）中方程的个数等于正则内点的个数，而未知数 jku , 则除了包含正则内点处 解u的近似值，还包含一些非正则内点处u的近似值，因而方程个数少于未知数个数。 在非正则内点处 Poisson 方程的差分近似不能按式（13）给出，需要利用边界条件得到。 边界条件的处理可以有各种 方案 气瓶 现场处置方案 .pdf气瓶 现场处置方案 .doc见习基地管理方案.doc关于群访事件的化解方案建筑工地扬尘治理专项方案下载 ，下面介绍较简单的两种。 (i) 直接转移 (ii) 线性插值 由式（13）所给出的差分格式称为五点菱形格式，实际计算时经常取 τ=h ，此时 五点菱形格式可化为 jkjkjkjkjkjk fuuuuuh ,,1,1,,1,12 )4(1 =−+++ −+−+ (14) 简记为 jkjk fuh ,,2 1 =◊ （15） 其中 jkjkjkjkjkjk uuuuuu ,1,1,,1,1, 4−+++=◊ −+−+ 。 求解差分方程组最常用的方法是同步迭代法，同步迭代法是最简单的迭代方式。除 边界节点外，区域内节点的初始值是任意取定的。 例 1 用五点菱形格式求解 Laplace 方程第一边值问题 ⎪⎩ ⎪⎨ ⎧ Ω∂=Γ++= Ω∈=∂ ∂+∂ ∂ Γ∈ ])1lg[(|),( ),( 0 22 ),( 2 2 2 2 yxyxu yx y u x u yx 其中 }1,0|),{( ≤≤=Ω yxyx 。取 3 1== τh 。 当 τ=h 时，利用点 )1,1(),1,1(),,( +±−± jkjkjk 构造的差分格式 jkjkjkjkjkjk fuuuuuh ,,1,11,11,11,12 )4( 2 1 =−+++ −−+−−+++ （16） -244- 称为五点矩形格式，简记为 22 1 h ٱ jkjk fu ,, = （17） 其中ٱ jkjkjkjkjkjk uuuuuu ,1,11,11,11,1, 4−+++= −−+−−+++ 。 2.2 抛物型方程的差分解法 以一维热传导方程（5） )0(02 2 >=∂ ∂−∂ ∂ a t ua t u 为基本模型讨论适用于抛物型方程定解问题的几种差分格式。 首先对 xt 平面进行网格剖分。分别取 τ,h 为 x方向与 t方向的步长，用两族平行直 线 ),2,1,0( L±±=== kkhxx k ， ),2,1,0( L=== jjtt j τ ，将 xt 平面剖分成矩形网 格，节点为 ),2,1,0,,2,1,0)(,( LL =±±= jktx jk 。为简便起见，记 ),(),( jk yxjk = ， ),(),( jk yxujku = ， )( kk xϕϕ = ， )(11 jj tgg = ， )(22 jj tgg = ， )(11 jj tλλ = ， )(22 jj tλλ = 。 2.2.1 微分方程的差分近似 在网格内点 ),( jk 处，对 t u ∂ ∂ 分别采用向前、向后及中心差商公式，对 2 2 x u ∂ ∂ 采用二 阶中心差商公式，一维热传导方程（5）可分别表示为 )(),1(),(2),1(),()1,( 22 hOh jkujkujkuajkujku +=−+−+−−+ ττ )(),1(),(2),1()1,(),( 22 hOh jkujkujkuajkujku +=−+−+−−− ττ )(),1(),(2),1( 2 )1,()1,( 2 2 hOh jkujkujkuajkujku +=−+−+−−−+ ττ 由此得到一维热传导方程的不同的差分近似 0 2 2 ,1,,1,1, =+−−− −++ h uuu a uu jkjkjkjkjk τ （18） 0 2 2 ,1,,11,, =+−−− −+− h uuu a uu jkjkjkjkjk τ （19） 0 2 2 2 ,1,,11,1, =+−−− −+−+ h uuu a uu jkjkjkjkjk τ （20） 2.2.2 初、边值条件的处理 为用差分方程求解定解问题（6），（7）等，还需对定解条件进行离散化。 对初始条件及第一类边界条件，可直接得到 kkk xuu ϕ== )0,(0, ),,1,0,1,0( nkk LL =±= 或 （21） jjjn ijjj gtluu gtuu 2, ,0 ),( ),0( == == )1,,1,0( −= mj L （22） -245- 其中 τ Tm h ln == , 。 对第二、三类边界条件则需用差商近似。下面介绍两种较简单的处理方法。 （i）在左边界 )0( =x 处用向前差商近似偏导数 x u ∂ ∂ ，在右边界 )( lx = 处用向后差 商近似偏导数 x u ∂ ∂ ，即 )(),1(),( )(),0(),1( ),( ),0( hO h jnujnu x u hO h juju x u jn j +−−=∂ ∂ +−=∂ ∂ ),,1,0( mj L= 即得边界条件（8）的差分近似为 ⎪⎪⎩ ⎪⎪⎨ ⎧ =+− =−− − jjnj jnjn jjj jj gu h uu gu h uu 2,2 ,1, 1,01 ,0,1 λ λ ),,1,0( mj L= （23） （ii）用中心差商近似 x u ∂ ∂ ，即 )( 2 ),1(),1( )( 2 ),1(),1( 2 ),( 2 ),0( hO h jnujnu x u hO h juju x u jn j +−−+=∂ ∂ +−−=∂ ∂ ),,1,0( mj L= 则得边界条件的差分近似为 ⎪⎪⎩ ⎪⎪⎨ ⎧ =+− =−− −+ − jjnj jnjn jjj jj gu h uu gu h uu 2,2 ,1,1 1,01 ,1,1 2 2 λ λ ),,1,0( mj L= （24） 这样处理边界条件，误差的阶数提高了，但式（24）中出现定解区域外的节点 ),1( j− 和 ),1( jn + ，这就需要将解拓展到定解区域外。可以通过用内节点上的u值插值求出 ju ,1− 和 jnu ,1+ ，也可以假定热传导方程（5）在边界上也成立，将差分方程扩展到边界节点 上，由此消去 ju ,1− 和 jnu ,1+ 。 2.2.3 几种常用的差分格式 下面我们以热传导方程的初边值问题（7）为例给出几种常用的差分格式。 (i) 古典显式格式 为便于计算，令 2h ar τ= ，式（18）改写成以下形式 -246- jkjkjkjk ruurruu ,1,,11, )21( −++ +−+= 将式（18）与（21），（22）结合，我们得到求解问题（7）的一种差分格式 ⎪⎩ ⎪⎨ ⎧ === == −=−=+−+= −++ ),,2,1( , ),,2,1( )1,,1,0,1,,2,1( )21( 2,1,0 0, ,1,,11, mjgugu nku mjnkruurruu jjnjj kk jkjkjkjk L L LL ϕ （25） 由于第 0 层 )0( =j 上节点处的u值已知 )( 0, kku ϕ= ，由式（25）即可算出u在第一层 )1( =j 上节点处的近似值 1,ku 。重复使用式（25），可以逐层计算出各层节点的近似值。 （ii）古典隐式格式 将（19）整理并与式（21），（22）联立，得差分格式如下 ⎪⎩ ⎪⎨ ⎧ === == −=−=+−+= +−++++ ),,1,0( , ),,1,0( )1,,1,0,1,,2,1)(2( 2,1,0 0, 1,11,1,1,1, mjgugu nku mjnkuuuruu jjnjj kk jkjkjkjkjk L L LL ϕ （26） 其中 2h ar τ= 。虽然第 0 层上的u值仍为已知，但不能由式（30）直接计算以上各层节 点上的值 jku , ，故差分格式（26）称为古典隐式格式。 （iii）杜福特—弗兰克尔（DoFort—Frankel）格式 DoFort—Frankel 格式是三层显式格式，它是由式（24）与（25），（26）结合得到 的。具体形式如下： ⎪⎪⎩ ⎪⎪⎨ ⎧ === == −=−=+ −+++= −−++ ),,1,0( , ),,1,0( )1,,2,1,1,,2,1( 21 21)( 21 2 2,1,0 0, 1,,1,11, mjgugu nku mjnku r ruu r ru jjnjj kk jkjkjkjk L L LL ϕ (27) 用这种格式求解时，除了第 0 层上的值 0,ku 由初值条件（21）得到，必须先用二层格式 求出第 1 层上的值 1,ku ，然后再按格式（27）逐层计算 ),,3,2(, mju jk L= 。 2.3 双曲型方程的差分解法 对二阶波动方程（10） 2 2 2 2 2 x ua t u ∂ ∂=∂ ∂ 如果令 x uv t uv ∂ ∂=∂ ∂= 21 , ，则方程（10）可化成一阶线性双曲型方程组 ⎪⎪⎩ ⎪⎪⎨ ⎧ ∂ ∂=∂ ∂ ∂ ∂=∂ ∂ x v t v x va t v 12 221 （28） 记 Tvvv ),( 21= ，则方程组（28）可表成矩阵形式 -247- x vA x va t v ∂ ∂=∂ ∂⎥⎦ ⎤⎢⎣ ⎡=∂ ∂ 01 0 2 （29） 矩阵 A有两个不同的特征值 a±=λ ，故存在非奇异矩阵P，使得 Λ=⎥⎦ ⎤⎢⎣ ⎡ −= − a a PAP 0 01 作变换 TwwPvw ),( 21== ，方程组（29）可化成 x w t w ∂ ∂Λ=∂ ∂ （30） 方程组（30）由两个独立的一阶双曲型方程联立而成。因此下面主要讨论一阶双曲型方 程的差分解法。 考虑一阶双曲型方程的初值问题 ⎪⎩ ⎪⎨ ⎧ +∞<<∞−= +∞<<∞−>=∂ ∂+∂ ∂ xxxu xt x ua t u )()0,( 0 0 ϕ (31) 与抛物型方程的讨论类似，仍将 xt 平面剖分成矩形网格。取 x方向步长为 h，t方向步 长为τ ，网格线为 ),,2,1,0( L±±=== kkhxx k ),2,1,0( L=== jjtt j τ 。为简便起 见，记 ),(),( jk yxjk = ， ),(),( jk yxujku = ， )( kk xϕϕ = 。 以不同的差商近似偏导数，可以得到方程（9）的不同的差分近似 0,,1,1, =−+− ++ h uu a uu jkjkjkjk τ （32） 0,1,,1, =−−− −+ h uu a uu jkjkjkjk τ （33） 0 2 ,1,1,1, =−−− −++ h uu a uu jkjkjkjk τ （34） 结合离散化的初始条件，可以得到几种简单的差分格式。 §3 一维状态空间的偏微分方程的 MATLAB 解法 3.1 工具箱命令介绍 MATLAB 提供了一个指令 pdepe，用以解以下的 PDE 方程式 ),,,()),,,((),,,( x uutxs x uutxfx x x t u x uutxc mm ∂ ∂+∂ ∂ ∂ ∂=∂ ∂ ∂ ∂ − （35） 其中时间介于 fttt ≤≤0 之间，而位置 x则介于 ],[ ba 有限区域之间。m值表示问题的 对称性，其可为 0，1 或 2，分别表示平板(slab)，圆柱(cylindrical)或球体(spherical)的情 形。因而，如果 0>m ，则a必等于b，也就是说其具有圆柱或球体的对称关系。同时， 式中 ),,,( xuutxf ∂∂ 一项为流通量 (flux)，而 ),,,( xuutxs ∂∂ 为来源 (source)项。 ),,,( xuutxc ∂∂ 为偏微分方程的对角线系数矩阵。若某一对角线元素为 0，则表示该偏 微分方程为椭圆型偏微分方程，若为正值(不为 0)，则为拋物型偏微分方程。请注意 c的 对角线元素一定不全为 0。偏微分方程初始值可表示为 -248- )(),( 00 xvtxu = （36） 而边界条件为 0),,,(),(),,( =∂∂+ xuutxftxqutxp （37） 其中 x为两端点位置，即a或b 用以解含上述初始值及边界值条件的偏微分方程的 MATLAB 命令 pdepe 的用法如 下： ),,,,,,( optionstspanxmeshbcfunicfunpdepempdepesol = 其中 m为问题之对称参数； xmesh为空间变量 x的网格点(mesh)位置向量，即 ],,,[ 10 Nxxxxmesh L= ，其中 ax =0 (起点)， bxN = (终点)。 tspan 为时间变量 t的向量，即 ],,,[ 10 Mttttspan L= ，其中 0t 为起始时间， Mt 为 终点时间。 pdefun为使用者提供的 pde 函数文件。其函数格式如下： [ ] ),,,(,, dudxutxpdefunsfc = 亦即，使用者仅需提供偏微分方程中的系数向量。 c， f 和 s均为行(column)向量，而 向量 c即为矩阵 c的对角线元素。 icfun提供解u的起始值，其格式为 )(xicfunu = 值得注意的是u为行向量。 bcfun使用者提供的边界条件函数，格式如下： [ ] ( )turxrulxlbcfunqlprqlpl ,,,,,,, = 其中ul和ur 分别表示左边界 )( bxl = 和右边界 )( axr = u的近似解。输出变量中，pl 和ql分别表示左边界 p和q的行向量，而 pr 和qr则为右边界 p和q的行向量。 sol 为解答输出。 sol 为多维的输出向量， ):(:, isol 为 iu 的输出，即 ):,(:, isolui = 。 元素 ),,(),( ikjsolkjui = 表示在 )( jtspant = 和 )(kxmeshx = 时 iu 之 答案 八年级地理上册填图题岩土工程勘察试题省略号的作用及举例应急救援安全知识车间5s试题及答案 。 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 的选取主要是基于使用者对那些特定时间的状态有兴趣而选定。而间距 -249- (step size)的控制由程序自动完成。 4. 若要获得特定位置及时间下的解，可配合以 pdeval 命令。使用格式如下： [ ] ),,,(, xoutuixmeshmpdevalduoutdxuout = 其中 m代表问题的对称性。m =0 表示平板；m =1 表示圆柱体；m =2 表示球体。其意 义同 pdepe 中的自变量m。 xmesh为使用者在 pdepe 中所指定的输出点位置向量。 0 1[ , , , ]Nxmesh x x x= L 。 ui即 ):,,( ijsol 。也就是说其为 pdepe 输出中第 i个输出ui在各点位置 xmesh处， 时间固定为 )( jtspant j = 下的解。 xout 为所欲内插输出点位置向量。此为使用者重新指定的位置向量。 uout为基于所指定位置 xout ，固定时间 ft 下的相对应输出。 duoutdx 为相对应的 dxdu 输出值。 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 ∂ ∂=∂ ∂π 其中 10 ≤≤ x ，且满足以下之条件限制式 (i)起始值条件 IC： )sin()0,( xxu π= (ii)边界条件 BC1： 0),0( =tu BC2： 0),1( =∂ ∂+− x tue tπ 注：本问题的解析解为 )sin(),( xetxu t π−= 解 下面将叙述求解的步骤与过程。当完成以下各步骤后，可进一步将其汇总为一 主程序 ex20_1.m，然后求解。 步骤 1 将欲求解的偏微分方程改写成如式的 标准 excel标准偏差excel标准偏差函数exl标准差函数国标检验抽样标准表免费下载红头文件格式标准下载 式。 0002 +⎟⎠ ⎞⎜⎝ ⎛ ∂ ∂ ∂ ∂=∂ ∂ x ux x x t uπ 此即 ( ) 2,,, π=∂∂ xuutxc ( ) x uxuutxf ∂ ∂=∂∂,,, ( ) 0,,, =∂∂ xuutxs 和 0=m 。 步骤 2 编写偏微分方程的系数向量函数。 function [c,f,s]=ex20_1pdefun(x,t,u,dudx) -250- c=pi^2; f=dudx; s=0; 步骤 3 编写起始值条件。 function u0=ex20_1ic(x) u0=sin(pi*x); 步骤 4 编写边界条件。在编写之前，先将边界条件改写成标准形式，如式（37）， 找出相对应的 )(⋅p 和 )(⋅q 函数，然后写出 MATLAB 的边界条件函数，例如，原边界条 件可写成 BC1： 0,0),0(0),0( ==∂ ∂⋅+ xt x tu BC2： 1,0),1(1 ==∂ ∂⋅+− xt x ue tπ 即 (0, ), 0,pl u t ql= = 和 , 1tpr e qrπ −= = 因而，边界条件函数可编写成 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') -251- 综合以上各步骤，可写成一个程序求解例 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=pi^2; f=dudx; s=0; %****************** -252- %初始条件函数 %****************** 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 试解以下联立的偏微分方程系统 )(024.0 212 1 2 1 uuF x u t u −−∂ ∂=∂ ∂ )(170.0 212 2 2 2 uuF x u t u −+∂ ∂=∂ ∂ 其中 ))(46.11exp())(73.5exp()( 212121 uuuuuuF −−−−=− ，且 10 ≤≤ x 和 0≥t 。 此联立偏微分方程系统满足以下初边值条件。 (i)初值条件 1)0,(1 =xu 0)0,(2 =xu (ii)边值条件 0),0(1 =∂ ∂ t x u 0),0(2 =tu 1),1(1 =tu 0),1(2 =∂ ∂ t x u 解 步骤 1：改写偏微分方程为标准式 ⎥⎦ ⎤⎢⎣ ⎡ − −−+ ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢ ⎣ ⎡ ∂ ∂∂ ∂ ∂ ∂=⎥⎦ ⎤⎢⎣ ⎡ ∂ ∂⎥⎦ ⎤⎢⎣ ⎡ • )( )( 170.0 024.0 * 1 1 21 21 2 1 2 1 uuF uuF x u x u xu u t 因此 ⎥⎦ ⎤⎢⎣ ⎡= 1 1 c ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢ ⎣ ⎡ ∂ ∂∂ ∂ = x u x u f 2 1 170.0 024.0 -253- ⎥⎦ ⎤⎢⎣ ⎡ − −−= )( )( 21 21 uuF uuF s 和 0=m 。另外，左边界条件( 0x = 处)。写成 ⎥⎦ ⎤⎢⎣ ⎡= ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢ ⎣ ⎡ ∂ ∂∂ ∂ ∗⎥⎦ ⎤⎢⎣ ⎡+⎥⎦ ⎤⎢⎣ ⎡ • 0 0 170.0 024.0 0 10 2 1 2 x u x u u 即 ⎥⎦ ⎤⎢⎣ ⎡= 2 0 u pl ， ⎥⎦ ⎤⎢⎣ ⎡= 0 1 ql 同理，右边界条件( 1x = 处)为 ⎥⎦ ⎤⎢⎣ ⎡= ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢ ⎣ ⎡ ∂ ∂∂ ∂ ∗⎥⎦ ⎤⎢⎣ ⎡+⎥⎦ ⎤⎢⎣ ⎡ − • 0 0 170.0 024.0 1 0 0 1 2 1 1 x u x u u 即 ⎥⎦ ⎤⎢⎣ ⎡ −= 0 11upr ， ⎥⎦ ⎤⎢⎣ ⎡= 1 0 qr 步骤 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 %*************************************** -254- %求解一维偏微分方程组的一个综合函数程序 %*************************************** 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 触煤反应装置内温度及转换率的分布 -255- 以外部热交换式的管形固定层触煤反应装置，进行苯加氢反应产生环己烷。此反应 系统之质量平衡及热平衡方程式如下： 0)(12 2 =Δ−+⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ ∂ ∂+∂ ∂+∂ ∂− p rBA p e GC Hr r T rr T GC k L T ρ 01 0 2 2 =+⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ ∂ ∂+∂ ∂+∂ ∂− Gy Mr r f rr f u D L f arBAe ρ 其中T 为温度(℃)， f 为反应率，L为轴向距离，r为径向距离。此系统的边界条件为 0=L , )(0 rTT = , )(0 rff = 0=r , 0=∂ ∂=∂ ∂ r f r T wrr = , )( wwe TThr Tk −=∂ ∂− wrr = , 0=∂ ∂ r f 此外，式中之相关数据及操作条件如下： (i)反应速率式 23H 苯 氢 环己烷 图 1 反应示意图 4 33 )1( CCBBHH BHBH A PKPKPK PPKkKr +++= 其中P表示分压(atm)，而速率参数为 RRTK 3.3212100ln +−= RRTKH 9.3115500ln −= RRTKB 1.2311200ln −= RRTKC 4.198900ln −= 上式中，下标 B，H 及 C 分别代表苯，氢及环己烷。R 为理想气体常数(1.987cal/mol·K)。 (ii)操作条件及物性数据 总压 atmPt 25.1= 反应管管径 cmrw 5.2= 壁温 100=wT ℃ 质量速度 hrmkgG ⋅= 2631 苯对氢之莫耳流量比 30=m 反应管入口的苯之莫耳分率 0323.00 =y 反应气体之平均分子量 47.4=avM -256- 触煤层密度 31200 mkgPB = 流体平均比热 molkgkcalCP −= 74.1 反应热 molkgkcalH r −−=Δ 49250 整体传热系数 Chrmkcalh 020 8.65 ⋅⋅= 进料温度 125)0( =T ℃ 反应管管长 cmL 45= 流速 hrmu 03.8= 有效热传导系数 ChrmkcalKe 065.0 ⋅⋅= 壁境膜传热系