偏微分方程的数值解方法及源程序

-240-第二十章偏微分方程的数值解自然科学与工程技术中种种运动发展过程与平衡现象各自遵守一定的规律。这些规律的定量表述一般地呈现为关于含有未知函数及其导数的方程。我们将只含有未知多元函数及其偏导数的方程，称之为偏微分方程。方程中出现的未知函数偏导数的昀高阶数称为偏微分方程的阶。如果方程中对于未知函数和它的所有偏导数都是线性的，这样的方程称为线性偏微分方程，否则称它为非线性偏微分方程。初始条件和边界条件称为定解条件，未附加定解条件的偏微分方程称为泛定方程。对于一个具体的问题，定解条件与泛定方程总是同时提出。定解条件与泛定方程作为一个整体，称为定解问题。§1偏微分方程的定解问题各种物理性质的定常（即不随时间变化）过程，都可用椭圆型方程来描述。其昀典型、昀简单的形式是泊松(Poisson)方程),(2222yxfyuxuu=∂∂+∂∂=Δ（1）特别地，当0),(≡yxf时，即为拉普拉斯(Laplace)方程，又称为调和方程02222=∂∂+∂∂=Δyuxuu（2）带有稳定热源或内部无热源的稳定温度场的温度分布，不可压缩流体的稳定无旋流动及静电场的电势等均满足这类方程。Poisson方程的第一边值问题为⎪⎩⎪⎨⎧Ω∂=Γ=Ω∈=∂∂+∂∂Γ∈),(|),(),(),(),(2222yxyxuyxyxfyuxuyxϕ（3）其中Ω为以Γ为边界的有界区域，Γ为分段光滑曲线，ΓΩU称为定解区域，),(),,(yxyxfϕ分别为ΓΩ,上的已知连续函数。第二类和第三类边界条件可统一表示成),(),(yxunuyxϕα=⎟⎠⎞⎜⎝⎛+∂∂Γ∈（4）其中n为边界Γ的外法线方向。当0=α时为第二类边界条件，0≠α时为第三类边界条件。在研究热传导过程，气体扩散现象及电磁场的传播等随时间变化的非定常物理问题时，常常会遇到抛物型方程。其昀简单的形式为一维热传导方程)0(022=∂∂−∂∂axuatu（5）方程（5）可以有两种不同类型的定解问题：初值问题（也称为Cauchy问题）-241-⎪⎩⎪⎨⎧+∞∞−=+∞∞−=∂∂−∂∂xxxuxtxuatu)()0,(,0022ϕ（6）初边值问题⎪⎪⎪⎩⎪⎪⎪⎨⎧≤≤====∂∂−∂∂TttgtlutgtuxxulxTtxuatu0),(),(),(),0()()0,(0,002122ϕ（7）其中)(),(),(21tgtgxϕ为已知函数，且满足连接条件)0()(),0()0(21glg==ϕϕ问题（7）中的边界条件)(),(),(),0(21tgtlutgtu==称为第一类边界条件。第二类和第三类边界条件为TttgutxuTttgutxulxx≤≤=⎥⎦⎤⎢⎣⎡+∂∂≤≤=⎥⎦⎤⎢⎣⎡−∂∂==0),()(0),()(22101λλ（8）其中0)(,0)(21≥≥ttλλ。当0)()(21≡=ttλλ时，为第二类边界条件，否则称为第三类边界条件。双曲型方程的昀简单形式为一阶双曲型方程0=∂∂+∂∂xuatu（9）物理中常见的一维振动与波动问题可用二阶波动方程22222xuatu∂∂=∂∂（10）描述，它是双曲型方程的典型形式。方程（10）的初值问题为⎪⎪⎪⎩⎪⎪⎪⎨⎧+∞∞−=∂∂+∞∞−=+∞∞−∂∂=∂∂=xxtuxxxuxtxuatut)()()0,(,0022222φϕ（11）边界条件一般也有三类，昀简单的初边值问题为-242-⎪⎪⎪⎩⎪⎪⎪⎨⎧≤≤==≤≤=∂∂=∂∂=∂∂=Tttgtlutgtulxxtuxxulxtxuatut0)(),(),(),0(0)(),()0,(0,021022222φϕ如果偏微分方程定解问题的解存在，唯一且连续依赖于定解数据（即出现在方程和定解条件中的已知函数），则此定解问题是适定的。可以证明，上面所举各种定解问题都是适定的。§2偏微分方程的差分解法差分方法又称为有限差分方法或网格法，是求偏微分方程定解问题的数值解中应用昀广泛的方法之一。它的基本思想是：先对求解区域作网格剖分，将自变量的连续变化区域用有限离散点（网格点）集代替；将问题中出现的连续变量的函数用定义在网格点上离散变量的函数代替；通过用网格点上函数的差商代替导数，将含连续变量的偏微分方程定解问题化成只含有限个未知数的代数方程组（称为差分格式）。如果差分格式有解，且当网格无限变小时其解收敛于原微分方程定解问题的解，则差分格式的解就作为原问题的近似解（数值解）。因此，用差分方法求偏微分方程定解问题一般需要解决以下问题：（i）选取网格；（ii）对微分方程及定解条件选择差分近似，列出差分格式；（iii）求解差分格式；（iv）讨论差分格式解对于微分方程解的收敛性及误差估计。下面我们只对偏微分方程的差分解法作一简要的介绍。2.1椭圆型方程第一边值问题的差分解法以Poisson方程（1）为基本模型讨论第一边值问题的差分方法。考虑Poisson方程的第一边值问题（3）⎪⎩⎪⎨⎧Ω∂=Γ=Ω∈=∂∂+∂∂Γ∈),(|),(),(),(),(2222yxyxuyxyxfyuxuyxϕ取τ,h分别为x方向和y方向的步长，以两族平行线τjyykhxxjk====,),2,1,0,(L±±=jk将定解区域剖分成矩形网格。节点的全体记为},,,|),{(为整数jijykhxyxRjkjkτ===。定解区域内部的节点称为内点，记内点集ΩIR为τhΩ。边界Γ与网格线的交点称为边界点，边界点全体记为τhΓ。与节点),(jkyx沿x方向或y方向只差一个步长的点),(1jkyx±和),(1±jkyx称为节点),(jkyx的相邻节点。如果一个内点的四个相邻节点均属于ΓΩU，称为正则内点，正则内点的全体记为)1(Ω，至少有一个相邻节点不属于ΓΩU的内点称为非正则内点，非正则内点的全体记为)2(Ω。我们的问题是要求出问题（3）在全体内点上的数值解。为简便记，记),(),,(),(),,(),(,jkjkjkjkyxffyxujkuyxjk===。对正则内点-243-)1(),(Ω∈jk，由二阶中心差商公式)(),1(),(2),1(22),(22hOhjkujkujkuxujk+−+−+=∂∂)()1,(),(2)1,(22),(22ττOjkujkujkuyujk+−+−+=∂∂Poisson方程（1）在点),(jk处可表示为)()1,(),(2)1,(),1(),(2),1(22,22ττ++=−+−++−+−+hOfjkujkujkuhjkujkujkujk（12）在式（12）中略去)(22τ+hO，即得与方程（1）相近似的差分方程jkjkjkjkjkjkjkfuuuhuuu,21,,1,2,1,,122=+−++−−+−+τ(13)式（13）中方程的个数等于正则内点的个数，而未知数jku,则除了包含正则内点处解u的近似值，还包含一些非正则内点处u的近似值，因而方程个数少于未知数个数。在非正则内点处Poisson方程的差分近似不能按式（13）给出，需要利用边界条件得到。边界条件的处理可以有各种方案，下面介绍较简单的两种。(i)直接转移(ii)线性插值由式（13）所给出的差分格式称为五点菱形格式，实际计算时经常取τ=h，此时五点菱形格式可化为jkjkjkjkjkjkfuuuuuh,,1,1,,1,12)4(1=−+++−+−+(14)简记为jkjkfuh,,21=◊（15）其中jkjkjkjkjkjkuuuuuu,1,1,,1,1,4−+++=◊−+−+。求解差分方程组昀常用的方法是同步迭代法，同步迭代法是昀简单的迭代方式。除边界节点外，区域内节点的初始值是任意取定的。例1用五点菱形格式求解Laplace方程第一边值问题⎪⎩⎪⎨⎧Ω∂=Γ++=Ω∈=∂∂+∂∂Γ∈])1lg[(|),(),(022),(2222yxyxuyxyuxuyx其中}1,0|),{(≤≤=Ωyxyx。取31==τh。当τ=h时，利用点)1,1(),1,1(),,(+±−±jkjkjk构造的差分格式jkjkjkjkjkjkfuuuuuh,,1,11,11,11,12)4(21=−+++−−+−−+++（16）-244-称为五点矩形格式，简记为221hٱjkjkfu,,=（17）其中ٱjkjkjkjkjkjkuuuuuu,1,11,11,11,1,4−+++=−−+−−+++。2.2抛物型方程的差分解法以一维热传导方程（5）)0(022=∂∂−∂∂atuatu为基本模型讨论适用于抛物型方程定解问题的几种差分格式。首先对xt平面进行网格剖分。分别取τ,h为x方向与t方向的步长，用两族平行直线),2,1,0(L±±===kkhxxk，),2,1,0(L===jjttjτ，将xt平面剖分成矩形网格，节点为),2,1,0,,2,1,0)(,(LL=±±=jktxjk。为简便起见，记),(),(jkyxjk=，),(),(jkyxujku=，)(kkxϕϕ=，)(11jjtgg=，)(22jjtgg=，)(11jjtλλ=，)(22jjtλλ=。2.2.1微分方程的差分近似在网格内点),(jk处，对tu∂∂分别采用向前、向后及中心差商公式，对22xu∂∂采用二阶中心差商公式，一维热传导方程（5）可分别表示为)(),1(),(2),1(),()1,(22hOhjkujkujkuajkujku+=−+−+−−+ττ)(),1(),(2),1()1,(),(22hOhjkujkujkuajkujku+=−+−+−−−ττ)(),1(),(2),1(2)1,()1,(22hOhjkujkujkuajkujku+=−+−+−−−+ττ由此得到一维热传导方程的不同的差分近似022,1,,1,1,=+−−−−++huuuauujkjkjkjkjkτ（18）022,1,,11,,=+−−−−+−huuuauujkjkjkjkjkτ（19）0222,1,,11,1,=+−−−−+−+huuuauujkjkjkjkjkτ（20）2.2.2初、边值条件的处理为用差分方程求解定解问题（6），（7）等，还需对定解条件进行离散化。对初始条件及第一类边界条件，可直接得到kkkxuuϕ==)0,(0,),,1,0,1,0(nkkLL=±=或（21）jjjnijjjgtluugtuu2,,0),(),0(====)1,,1,0(−=mjL（22）-245-其中τTmhln==,。对第二、三类边界条件则需用差商近似。下面介绍两种较简单的处理方法。（i）在左边界)0(=x处用向前差商近似偏导数xu∂∂，在右边界)(lx=处用向后差商近似偏导数xu∂∂，即)(),1(),()(),0(),1(),(),0(hOhjnujnuxuhOhjujuxujnj+−−=∂∂+−=∂∂),,1,0(mjL=即得边界条件（8）的差分近似为⎪⎪⎩⎪⎪⎨⎧=+−=−−−jjnjjnjnjjjjjguhuuguhuu2,2,1,1,01,0,1λλ),,1,0(mjL=（23）（ii）用中心差商近似xu∂∂，即)(2),1(),1()(2),1(),1(2),(2),0(hOhjnujnuxuhOhjujuxujnj+−−+=∂∂+−−=∂∂),,1,0(mjL=则得边界条件的差分近似为⎪⎪⎩⎪⎪⎨⎧=+−=−−−+−jjnjjnjnjjjjjguhuuguhuu2,2,1,11,01,1,122λλ),,1,0(mjL=（24）这样处理边界条件，误差的阶数提高了，但式（24）中出现定解区域外的节点),1(j−和),1(jn+，这就需要将解拓展到定解区域外。可以通过用内节点上的u值插值求出ju,1−和jnu,1+，也可以假定热传导方程（5）在边界上也成立，将差分方程扩展到边界节点上，由此消去ju,1−和jnu,1+。2.2.3几种常用的差分格式下面我们以热传导方程的初边值问题（7）为例给出几种常用的差分格式。(i)古典显式格式为便于计算，令2harτ=，式（18）改写成以下形式-246-jkjkjkjkruurruu,1,,11,)21(−+++−+=将式（18）与（21），（22）结合，