170第9章 偏微分方程的差分方法含有偏导数的微分方程称为偏微分方程。

由于变量的增多和区域的复杂性，求偏微分方程的精确解一般是不可能的，经常采用数值方法求方程的近似解。

偏微分方程的数值方法种类较多，最常用的方法是差分方法。

差分方法具有格式简单，程序易于实现，计算量小等优点，特别适合于规则区域上偏微分方程的近似求解。

本章将以一些典型的偏微分方程为例，介绍差分方法的基本原理和具体实现方法。

9.1椭圆型方程边值问题的差分方法9.1.1 差分方程的建立最典型的椭圆型方程是Poisson （泊松）方程G y x y x f yux u u ∈=∂∂+∂∂-≡∆-),(),,()(2222 （9.1）G 是x ,y 平面上的有界区域，其边界Γ为分段光滑的闭曲线。

当f (x ,y )≡0时，方程（9.1）称为Laplace(拉普拉斯）方程。

椭圆型方程的定解条件主要有如下三种边界条件第一边值条件 ),(y x u α=Γ （9.2） 第二边值条件),(y x nuβ=∂∂Γ （9.3） 第三边值条件 ),()(y x ku nuγ=+∂∂Γ （9.4） 这里，n 表示Γ上单位外法向，α(x,y ),β(x,y ),γ(x,y )和k (x,y )都是已知的函数，k (x,y )≥0。

满足方程（9.1）和上述三种边值条件之一的光滑函数u (x ,y )称为椭圆型方程边值问题的解。

用差分方法求解偏微分方程，就是要求出精确解u (x ,y )在区域G 的一些离散节点（x i ，y i ）上的近似值u i ,j ≈(x i ，y i ）。

差分方法的基本思想是，对求解区域G 做网格剖分，将偏微分方程在网格节点上离散化，导出精确解在网格节点上近似值所满足的差分方程，最终通过求解差分方程，通常为一个线性方程组，得到精确解在离散节点上的近似值。

设G ={0<x <a , 0<y <b }为矩形区域，在x ,y 平面上用两组平行直线x =ih 1, i =0,1,…,N 1, h 1=a /N 1 y =jh 2, j =0,1,…,N 2, h 2=b /N 2将G 剖分为网格区域，见图9-1。

h 1，h 2分别称为x 方向和y 方向的剖分步长，网格交点（x i ，y i ）称为剖分节点（区域内节点集合记为G h ={（x i ，y i ）; （x i ，y i ）∈G }）,网格线与边界Γ的交点称为边界点，边界点集合记为Γh 。

171现在将微分方程（9.1）在每一个内节点（x i ，y i ）上进行离散。

在节点（x i ，y i ）处，方程（9.1）为h i i i i i i i i G y x y x f y x yuy x x u ∈=∂∂+∂∂-),(),,()],(),([2222 （9.5） 需进一步离散（9.5）中的二阶偏导数。

为简化记号，简记节点（x i ，y i ）=(i ,j )，节点函数值u （x i ，y i ）=u (i ,j )。

利用一元函数的Taylor 展开公式，推得二阶偏导数的差商表达式)(0)]1,(),(2)1,([1),()(0)],1(),(2),1([1),(222222212122h j i u j i u j i u h j i y u h j i u j i u j i u h j i x u +-+-++=∂∂+-+-++=∂∂代入（9.5）式中，得到方程（9.1）在节点(i ,j )处的离散形式h j i G j i h h f j i u j i u j i u h j i u j i u j i u h ∈++=-+-+--+-+-),(),(0)]1,(),(2)1,([1)],1(),(2),1([12221,2221其中),(,i i j i y x f f =。

舍去高阶小项)(02221h h +，就导出了u (i ,j )的近似值u i ,j 所满足的差分方程h j i j i j i j i j i j i j i G j i f u u u h u u u h ∈=+--+---+-+),(,]2[1]2[1,1,,1,22,1,,121 （9.6） 在节点(i ,j )处方程（9.6）逼近偏微分方程（9.1）的误差为)(2221h h O +，它关于剖分步长是二阶的。

这个误差称为差分方程逼近偏微分方程的截断误差，它的大小将影响近似解的精度。

在差分方程（9.6）中，每一个节点(i ,j )处的方程仅涉及五个节点未知量u i ,j ，u i +1,j ，u i -1,j ，u i ,j +1，u i ,j -1，因此通常称（9.6）式为五点差分格式，当h 1= h 2=h 时，它简化为172h j i j i j i j i j i j i G j i f u u u u u h∈=-+++--+-+),(,]4[1,,1,1,,1,12 差分方程（9.6）中，方程个数等于内节点总数，但未知量除内节点值u i ,j ,(i ,j )∈G h 外，还包括边界点值。

例如，点(1,j )处方程就含有边界点未知量u 0,j 。

因此，还要利用给定的边值条件补充上边界点未知量的方程。

对于第一边值条件式（9.2），可直接取u i ,j =α（x i ，y i ）， (i ,j )∈Γh （9.7） 对于第三（k =0时为第二）边值条件式（9.4）， 以左边界点(1,j )为例，见图9-2, 利用一阶差商公式)(),1(),0(),0(11h O h j u j u j n u +-=∂∂ 则得到边界点(0,j )处的差分方程j j j jj r u k h u u ,0,0,01,1,0=+- （9.8）联立差分方程（9.6）与（9.7）或（9.8）就形成了求解Poisson 方程边值问题的差分方程组，它实质上是一个关于未知量{u i ,j }的线性代数方程组，可采用第2，3章介绍的方法进行求解。

这个方程组的解就称为偏微分方程的差分近似解，简称差分解。

考虑更一般形式的二阶椭圆型方程G y x y x f Eu yu D x u C y u B y x u A x ∈=+∂∂+∂∂+∂∂∂∂+∂∂∂∂-),(),,(])()([（9.9） 其中A (x ,y )≥A m in >0, B (x ,y ) ≥B m in >0, E(x ,y ) ≥0。

引进半节点,12121h x xi i ±=±,22121h y yi i ±=±利用一阶中心差商公式，在节点（i ,j ）处可有)(2),1(),1(),()(]),1(),(),(),1([1)()],21)((),21)([(1),)((211211,211,211211h O h j i u j i u j i x u h O h j i u j i u A h j i u j i u A h h O j i x u A j i x u A h j i x u A x j i j i +--+=∂∂+----+=+-∂∂-+∂∂=∂∂∂∂-+173对yu y u B y ∂∂∂∂∂∂),(类似处理，就可推得求解方程（9.9）的差分方程 hj i j i j i j i j i j i j i j i j i j i G j i j i f u a u a u a u a u a ∈=-+++---++---+),(),,(][,,1,1,1,1,,1,1,1,1 （9.10）其中⎪⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎪⎨⎧++++=-=+=-=+=-+--+----+-+---+-+ji j i j i j i j i j i j i j i j i j i j i j i j i j i j i j i j i j i E B B h A A h a D h B h a D h B h a C h A h a C h A h a ,21,21,22,21,2121,,221,221,,221,221,,1,2121,1,1,2121,1)()()2()2()2()2( （9.11） 显然，当系数函数A (x ,y )=B (x ,y )=1, C (x ,y )=D (x ,y )=E (x ,y )=0时，椭圆型方程（9.9）就成为Poisson 方程（9.1），而差分方程（9.10）就成为差分方程（9.6）。

容易看出，差分方程（9.10）的截断误差为)(2221h h O +阶。

9.1.2 一般区域的边界条件处理前面已假设G 为矩形区域，现在考虑G 为一般区域情形，这里主要涉及边界条件的处理。

考虑Poisson 方程第一边值问题⎩⎨⎧Γ∈=∈=∆-),(),,(),(),,(y x y x u Gy x y x f u α （9.12） 其中G 可为平面上一般区域，例如为曲边区域。

仍然用两组平行直线：x =x 0+ih 1,y =y 0+jh 2,i ,j =0,±1,…,对区域G 进行矩形网格剖分，见图9-3。

如果一个内节点（i ,j ）的四个相邻节点（i +1,j ），（i -1,j ），（i ,j +1）和（i ,j -1）属于Γ⋃=G G ，则称其为正则内点，见图9-3中打“。

”号者；如果一个节点（i ,j ）属于G 且不为正则内点，则称其为非正则内点，见图9-3中打“.”号者。

记正则内点集合为hG '，非正则内点集合为h Γ'。

显然，当G 为矩形区域时，174h h h hG G Γ=Γ'=',成立。

在正则内点（i ,j ）处，完全同矩形区域情形，可建立五点差分格式h j i j i j i j i j i j i j i G j i f u u u h u u u h '∈=+--+---+-+),(,]2[1]2[1,1,,1,22,1,,121 （9.13） 在方程（9.13）中，当（i ,j ）点临近边界时，将出现非正则内点上的未知量，因此必须补充非正则内点处的方程。

若非正则内点恰好是边界点，如图9-4中 D 点，则利用边界条件可取u D =α(D)对于不是边界点的非正则内点，如图9-4中B 点，一般可采用如下两种处理方法。

a.直接转移法.取与点B 距离最近的边界点（如图9-4中E 点）上的u 的值作为u (B )的近似值u B ，即u B =u (E)=α(E)直接转移法的优点是简单易行，但精度较低，只为一阶近似。

b .线性插值法.取B 点的两个相邻点（如图9-4中边界点A 和正则内点C 作为插值节点对u (B )进行线性插值)()()()(21h O C u x x x x A u x x x x B u AC AB AC B C +--+--=则得到点B 处的方程 A B C B x x u h A h h u -=+++=δδδαδ,)(111 线性插值法精度较高，为二阶近似。