第9章 偏微分方程的差分方法
含有偏导数的微分方程称为偏微分方程。由于变量的增多和区域的复杂性,求偏微分方程的精确解一般是不可能的,经常采用数值方法求方程的近似解。偏微分方程的数值方法种类较多,最常用的方法是差分方法。差分方法具有格式简单,程序易于实现,计算量小等优点,特别适合于规则区域上偏微分方程的近似求解。本章将以一些典型的偏微分方程为例,介绍差分方法的基本原理和具体实现方法。
9.1椭圆型方程边值问题的差分方法
9.1.1 差分方程的建立
最典型的椭圆型方程是Poisson(泊松)方程
?2u?2u ??u??(2?2)?f(x,y),(x,y)?G (9.1)
?x?yG是x,y平面上的有界区域,其边界Γ为分段光滑的闭曲线。当f(x,y)≡0时,方程
(9.1)称为Laplace(拉普拉斯)方程。椭圆型方程的定解条件主要有如下三种边界条件
第一边值条件 u???(x,y) (9.2) 第二边值条件
?u???(x,y) (9.3) ?n?u?ku)???(x,y) (9.4) ?n 第三边值条件 (这里,n表示Γ上单位外法向,α(x,y),β(x,y),γ(x,y)和k(x,y)都是已知的函数,k(x,y)≥0。满足方程(9.1)和上述三种边值条件之一的光滑函数u(x,y)称为椭圆型方程边值问题的解。
用差分方法求解偏微分方程,就是要求出精确解u(x,y)在区域G的一些离散节点(xi,yi)上的近似值ui,j≈(xi,yi)。差分方法的基本思想是,对求解区域G做网格剖分,将偏微分方程在网格节点上离散化,导出精确解在网格节点上近似值所满足的差分方程,最终通过求解差分方程,通常为一个线性方程组,得到精确解在离散节点上的近似值。
设G={0 x=ih1, i=0,1,…,N1, h1=a/N1 y=jh2, j=0,1,…,N2, h2=b/N2 将G剖分为网格区域,见图9-1。h1,h2分别称为x方向和y方向的剖分步长,网格交点(xi,yi)称为剖分节点(区域内节点集合记为Gh={(xi,yi); (xi,yi)∈G}),网格线与边界Γ的交点称为边界点,边界点集合记为Γh。 . 现在将微分方程(9.1)在每一个内节点(xi,yi)上进行离散。在节点(xi,yi)处,方程(9.1)为 ?2u?2u?[2(xi,yi)?2(xi,yi)]?f(xi,yi),(xi,yi)?Gh (9.5) ?x?y需进一步离散(9.5)中的二阶偏导数。为简化记号,简记节点(xi,yi)=(i,j),节 点函数值u(xi,yi)=u(i,j)。利用一元函数的Taylor展开公式,推得二阶偏导数的差商表达式 ?2u1(i,j)??[u(i?1,j)?2u(i,j)?u(i?1,j)]?0(h12)22?xh1?u12(i,j)??[u(i,j?1)?2u(i,j)?u(i,j?1)]?0(h2)22?yh2代入(9.5)式中,得到方程(9.1)在节点(i,j)处的离散形式 2 11[u(i?1,j)?2u(i,j)?u(i?1,j)]?[u(i,j?1)?2u(i,j)?u(i,j?1)]22h2 h1 2?fi,j?0(h12?h2),(i,j)?Gh?2其中fi,j?f(xi,yi)。舍去高阶小项0(h12?h2),就导出了u(i,j)的近似值ui,j所满 足的差分方程 ?11 [u?2u?u]?[ui,j?1?2ui,j?ui,j?1]?fi,j,(i,j)?Gh (9.6)i?1,ji,ji?1,j22h1h222在节点(i,j)处方程(9.6)逼近偏微分方程(9.1)的误差为O(h1?h2),它关于剖分步长是二阶的。这个误差称为差分方程逼近偏微分方程的截断误差,它的大小将影响近似解的精度。 在差分方程(9.6)中,每一个节点(i,j)处的方程仅涉及五个节点未知量ui,j,ui+1,j,ui-1,j,ui,j+1,ui,j-1,因此通常称(9.6)式为五点差分格式,当h1= h2=h时,它简化为 . ?1[ui?1,j?ui?1,j?ui,j?1?ui,j?1?4ui,j]?fi,j,(i,j)?Gh 2h差分方程(9.6)中,方程个数等于内节点总数,但未知量除内节点值ui,j ,(i,j)∈Gh外,还包括边界点值。例如,点(1,j)处方程就含有边界点未知量u0,j。因此,还要利用给定的边值条件补充上边界点未知量的方程。 对于第一边值条件式(9.2),可直接取ui,j=α(xi,yi), (i,j)∈Γh (9.7) 对于第三(k=0时为第二)边值条件式(9.4), 以左边界点(1,j)为例,见图9-2, 利用一阶差商公式 ?uu(0,j)?u(1,j)(0,j)??O(h1) ?nh1则得到边界点(0,j)处的差分方程 u0,j?u1,jh1?k0,ju0,j?r0,j (9.8) 联立差分方程(9.6)与(9.7)或(9.8)就形成了求解Poisson方程边值问题的差分方程组,它实质上是一个关于未知量{ui,j}的线性代数方程组,可采用第2,3章介绍的方法进行求解。这个方程组的解就称为偏微分方程的差分近似解,简称差分解。 考虑更一般形式的二阶椭圆型方程 ?[??u??u?u?u(A)?(B)?C?D?Eu]?f(x,y),(x,y)?G (9.9) ?x?x?y?y?x?y其中A(x,y)≥Amin>0, B(x,y) ≥Bmin >0, E(x,y) ≥0。引进半节点xi?121?xi?h1, 2yi?121?yi?h2,利用一阶中心差商公式,在节点(i,j)处可有 2??u1?u1?u1(A)(i,j)?[(A)(i?,j)?(A)(i?,j)]?O(h12)?x?xh1?x2?x2?1u(i?1,j)?u(i,j)u(i,j)?u(i?1,j)[A1?A1]?O(h12)i?,jh1i?2,jh1h12?uu(i?1,j)?u(i?1,j)(i,j)??O(h12)?x2h1. 对 ??u?u类似处理,就可推得求解方程(9.9)的差分方程 (B),?y?y?y?[ai?1,jui?1,j?ai?1,jui?1,j?ai,j?1ui,j?1?ai,j?1ui,j?1?ai,jui,j]?f(i,j),(i,j)?Gh其中 (9.10) h1??2a?h(A?Ci,j)11?i?1,ji?,j22?h1?2?a?h(A?Ci,j)11?i?1,ji?,j22??h?2 (9.11) (B1?2Di,j)?ai,j?1?h2i,j?2?2?h2?2a?h(B?Di,j)?i,j?121i,j?22??2?2?ai,j?h1(A1?A1)?h2(B1?B1)?Ei,ji?,ji?,ji,j?i,j??2222?显然,当系数函数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)。容易看 2出,差分方程(9.10)的截断误差为O(h12?h2)阶。 9.1.2 一般区域的边界条件处理 前面已假设G为矩形区域,现在考虑G为一般区域情形,这里主要涉及边界条件的处理。 考虑Poisson方程第一边值问题 ???u?f(x,y),(x,y)?G (9.12) ?u??(x,y),(x,y)???其中G可为平面上一般区域,例如为曲边区域。仍然用两组平行直线: x=x0+ih1,y=y0+jh2,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中打“.”号者。记正则 ?,非正则内点集合为?h?。显然,当G为矩形区域时,内点集合为Gh. ??Gh,?h???h成立。 Gh 在正则内点(i,j)处,完全同矩形区域情形,可建立五点差分格式 ?11? (9.13)[u?2u?u]?[ui,j?1?2ui,j?ui,j?1]?fi,j,(i,j)?Gh i?1,ji,ji?1,j22h1h2在方程(9.13)中,当(i,j)点临近边界时,将出现非正则内点上的未知量,因此必须补充非正则内点处的方程。 若非正则内点恰好是边界点,如图9-4中 D点,则利用边界条件可取uD=α(D) 对于不是边界点的非正则内点,如图9-4中B点,一般可采用如下两种处理方法。 a.直接转移法.取与点B距离最近的边界点(如图9-4中E点)上的u的值作为u(B)的近似值uB,即uB=u(E)=α(E) 直接转移法的优点是简单易行,但精度较低,只为一阶近似。 b.线性插值法.取B点的两个相邻点(如图9-4中边界点A和正则内点C作为插值节点对u(B)进行线性插值 u(B)?xC?xBx?xAu(A)?Bu(C)?O(h12) xC?xAxC?xA则得到点B处的方程 uB?h1??(A)?uC,??xB?xA h1??h1??线性插值法精度较高,为二阶近似。 对每一个非正则内点进行上述处理,将所得到的方程与(9.13)式联立,就组成了方程个数与未知量个数相一致的线性代数方程组。求解此方程组就可得到一般 .