《偏微分方程的数值解方法及源程序资料》由会员分享，可在线阅读，更多相关《偏微分方程的数值解方法及源程序资料（30页珍藏版）》请在金锄头文库上搜索。
1、 -240- 第二十章第二十章 偏微分方程的数值解偏微分方程的数值解 自然科学与工程技术中种种运动发展过程与平衡现象各自遵守一定的规律。 这些规 律的定量表述一般地呈现为关于含有未知函数及其导数的方程。 我们将只含有未知多元 函数及其偏导数的方程，称之为偏微分方程。 方程中出现的未知函数偏导数的最高阶数称为偏微分方程的阶。 如果方程中对于未 知函数和它的所有偏导数都是线性的， 这样的方程称为线性偏微分方程， 否则称它为非 线性偏微分方程。 初始条件和边界条件称为定解条件，未附加定解条件的偏微分方程称为泛定方程。 对于一个具体的问题， 定解条件与泛定方程总是同时提出。 定解条件与泛定方程作为一
2、个整体，称为定解问题。 1 偏微分方程的定解问题 各种物理性质的定常（即不随时间变化）过程，都可用椭圆型方程来描述。其最典 型、最简单的形式是泊松(Poisson)方程 ),( 2 2 2 2 yxf y u x u u= + = （1） 特别地，当0),(yxf时，即为拉普拉斯(Laplace)方程，又称为调和方程 0 2 2 2 2 = + = y u x u u （2） 带有稳定热源或内部无热源的稳定温度场的温度分布， 不可压缩流体的稳定无旋流动及 静电场的电势等均满足这类方程。 Poisson 方程的第一边值问题为 = = + ),(| ),( ),(),( ),( 2 2 2 2 y
3、xyxu yxyxf y u x u yx （3） 其中为以为边界的有界区域，为分段光滑曲线，U称为定解区域， ),(),(yxyxf分别为,上的已知连续函数。 第二类和第三类边界条件可统一表示成 ),( ),( yxu n u yx = + （4） 其中n为边界的外法线方向。当0=时为第二类边界条件，0时为第三类边界 条件。 在研究热传导过程，气体扩散现象及电磁场的传播等随时间变化的非定常物理问 题时，常常会遇到抛物型方程。其最简单的形式为一维热传导方程 )0(0 2 2 = a x u a t u （5） 方程（5）可以有两种不同类型的定解问题： 初值问题（也称为 Cauchy 问题） -
4、241- += a t u a 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=，)( 11jj tgg=，)( 22jj tgg=，)( 11jj t=， )( 22jj t=。 2.2.1 微分方程的差
5、分近似 在网格内点),(jk处，对 t u 分别采用向前、向后及中心差商公式，对 2 2 x u 采用二 阶中心差商公式，一维热传导方程（5）可分别表示为 )( ), 1(),(2), 1(),() 1,( 2 2 hO h jkujkujku a jkujku += + + )( ), 1(),(2), 1() 1,(),( 2 2 hO h jkujkujku a jkujku += + )( ), 1(),(2), 1( 2 ) 1,() 1,( 2 2 hO h jkujkujku a jkujku += + + 由此得到一维热传导方程的不同的差分近似 0 2 2 , 1, 1,1,
6、= + + 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(nkkLL=或 （21） jjjn ijjj gtluu gtuu 2, , 0 ),( ), 0( = = ) 1, 1
7、 , 0(=mjL （22） -245- 其中 T m h l n=,。 对第二、三类边界条件则需用差商近似。下面介绍两种较简单的处理方法。 （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(mjL= 即得边界条件（8）的差分近似为 =+ = jjnj jnjn jjj jj gu h uu gu h uu 2,2 , 1, 1, 01 , 0, 1
8、 ), 1 , 0(mjL= （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(mjL= 则得边界条件的差分近似为 =+ = + jjnj jnjn jjj jj gu h uu gu h uu 2,2 , 1, 1 1, 01 , 1, 1 2 2 ), 1 , 0(mjL= （24） 这样处理边界条件，误差的阶数提高了，但式（24）中出现定解区域外的节点), 1(j和 ), 1(jn +， 这
9、就需要将解拓展到定解区域外。 可以通过用内节点上的u值插值求出 j u , 1 和 jn u , 1+ ，也可以假定热传导方程（5）在边界上也成立，将差分方程扩展到边界节点 上，由此消去 j u , 1 和 jn u , 1+ 。 2.2.3 几种常用的差分格式 下面我们以热传导方程的初边值问题（7）为例给出几种常用的差分格式。 (i) 古典显式格式 为便于计算，令 2 h a r =，式（18）改写成以下形式 -246- jkjkjkjk ruurruu , 1, 11, )21 ( + += 将式（18）与（21） ， （22）结合，我们得到求解问题（7）的一种差分格式 = = =+= +
10、 ), 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,kk u=，由式（25）即可算出u在第一层 ) 1( =j上节点处的近似值 1 ,k u。 重复使用式 （25） ， 可以逐层计算出各层节点的近似值。 （ii）古典隐式格式 将（19）整理并与式（21） ， （22）联立，得差分格式如下 = = =+= + ), 1 , 0( , ), 1 , 0
11、( ) 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） 其中 2 h a r =。虽然第 0 层上的u值仍为已知，但不能由式（30）直接计算以上各层节 点上的值 jk u , ，故差分格式（26）称为古典隐式格式。 （iii）杜福特弗兰克尔（DoFortFrankel）格式 DoFortFrankel 格式是三层显式格式，它是由式（24）与（25） ， （26）结合得到 的。具体形式如下： = = = + + + = + ), 1 , 0(
12、 , ), 1 , 0( ) 1, 2 , 1, 1, 2 , 1( 21 21 )( 21 2 2,1, 0 0 , 1, 1, 11, mjgugu nku mjnku r r uu r r u jjnjj kk jkjkjkjk L L LL (27) 用这种格式求解时，除了第 0 层上的值 0 ,k u由初值条件（21）得到，必须先用二层格式 求出第 1 层上的值 1 ,k u，然后再按格式（27）逐层计算), 3 , 2( , mju jk L=。 2.3 双曲型方程的差分解法 对二阶波动方程（10） 2 2 2 2 2 x u a t u = 如果令 x u v t u v = =
13、 21 ,，则方程（10）可化成一阶线性双曲型方程组 = = x v t v x v a t v 12 221 （28） 记 T vvv),( 21 =，则方程组（28）可表成矩阵形式 -247- x v A x va t v = = 01 0 2 （29） 矩阵A有两个不同的特征值a=，故存在非奇异矩阵P，使得 = = a a PAP 0 0 1 作变换 T wwPvw),( 21 =，方程组（29）可化成 x w t w = （30） 方程组（30）由两个独立的一阶双曲型方程联立而成。因此下面主要讨论一阶双曲型方 程的差分解法。 考虑一阶双曲型方程的初值问题 +m，则a必等于b，也就是说其
14、具有圆柱或球体的对称关系。同时， 式中),(xuutxf一项为流通量(flux)，而),(xuutxs为来源(source)项。 ),(xuutxc为偏微分方程的对角线系数矩阵。若某一对角线元素为 0，则表示该偏 微分方程为椭圆型偏微分方程， 若为正值(不为 0)， 则为拋物型偏微分方程。 请注意c的 对角线元素一定不全为 0。偏微分方程初始值可表示为 -248- )(),( 00 xvtxu= （36） 而边界条件为 0),(),(),(=+xuutxftxqutxp （37） 其中x为两端点位置，即a或b 用以解含上述初始值及边界值条件的偏微分方程的 MATLAB 命令 pdepe 的用法如 下： ),(optionstspanxmeshbcfunicfunpdep
