第二十章偏微分方程的数值解 自然科学与工程技术中种种运动发展过程与平衡现象各自遵守一定的规律。这些规 律的定量表述一般地呈现为关于含有未知函数及其导数的方程。我们将只含有未知多元 函数及其偏导数的方程,称之为偏微分方程。 方程中出现的未知函数偏导数的最高阶数称为偏微分方程的阶。如果方程中对于未 知函数和它的所有偏导数都是线性的,这样的方程称为线性偏微分方程,否则称它为非 线性偏微分方程 初始条件和边界条件称为定解条件,未附加定解条件的偏微分方程称为泛定方程 对于一个具体的问题,定解条件与泛定方程总是同时提出。定解条件与泛定方程作为 个整体,称为定解问题 S1偏微分方程的定解问题 各种物理性质的定常(即不随时间变化)过程,都可用椭圆型方程来描述。其最典 型、最简单的形式是泊松( Poisson)方程 △t= =f(, y) (1) 特别地,当∫(x,y)≡0时,即为拉普拉斯( Laplace)方程,又称为调和方程 △= =0 带有稳定热源或内部无热源的稳定温度场的温度分布,不可压缩流体的稳定无旋流动及 静电场的电势等均满足这类方程。 Poisson方程的第一边值问题为 2l+=f(x,y)(xy)∈9 u(x, y)Ix,per=p(x, y) I=aQ2 其中Ω为以为边界的有界区域,T为分段光滑曲线,ΩUr称为定解区域, f(x,y),(x,y)分别为Ω,r上的已知连续函数。 第二类和第三类边界条件可统一表示成 +asyer=p(x,y) (4) 其中n为边界r的外法线方向。当a=0时为第二类边界条件,a≠0时为第三类边界 条件 在研究热传导过程,气体扩散现象及电磁场的传播等随时间变化的非定常物理问 题时,常常会遇到抛物型方程。其最简单的形式为一维热传导方程 Ou a-u 方程(5)可以有两种不同类型的定解问题 初值问题(也称为 Cauchy问题)
-240- 第二十章 偏微分方程的数值解 自然科学与工程技术中种种运动发展过程与平衡现象各自遵守一定的规律。这些规 律的定量表述一般地呈现为关于含有未知函数及其导数的方程。我们将只含有未知多元 函数及其偏导数的方程，称之为偏微分方程。 方程中出现的未知函数偏导数的最高阶数称为偏微分方程的阶。如果方程中对于未 知函数和它的所有偏导数都是线性的，这样的方程称为线性偏微分方程，否则称它为非 线性偏微分方程。 初始条件和边界条件称为定解条件，未附加定解条件的偏微分方程称为泛定方程。 对于一个具体的问题，定解条件与泛定方程总是同时提出。定解条件与泛定方程作为一 个整体，称为定解问题。 §1 偏微分方程的定解问题 各种物理性质的定常（即不随时间变化）过程，都可用椭圆型方程来描述。其最典 型、最简单的形式是泊松(Poisson)方程 ( , ) 2 2 2 2 f x y y u x u u = ∂ ∂ + ∂ ∂ Δ = （1） 特别地，当 f (x, y) ≡ 0 时，即为拉普拉斯(Laplace)方程，又称为调和方程 0 2 2 2 2 = ∂ ∂ + ∂ ∂ Δ = y u x u u （2） 带有稳定热源或内部无热源的稳定温度场的温度分布，不可压缩流体的稳定无旋流动及 静电场的电势等均满足这类方程。 Poisson 方程的第一边值问题为 ⎪ ⎩ ⎪ ⎨ ⎧ = Γ = ∂Ω = ∈Ω ∂ ∂ + ∂ ∂ ∈Γ ( , ) | ( , ) ( , ) ( , ) ( , ) 2 2 2 2 u x y x y f x y x y y u x u x y ϕ （3） 其中 Ω 为以 Γ 为边界的有界区域， Γ 为分段光滑曲线， Ω U Γ 称为定解区域， f (x, y),ϕ(x, y) 分别为Ω,Γ 上的已知连续函数。 第二类和第三类边界条件可统一表示成 ( , ) ( , ) u x y n u α ⎟ x y = ϕ ⎠ ⎞ ⎜ ⎝ ⎛ + ∂ ∂ ∈Γ （4） 其中 n 为边界Γ 的外法线方向。当α = 0 时为第二类边界条件，α ≠ 0时为第三类边界 条件。 在研究热传导过程，气体扩散现象及电磁场的传播等随时间变化的非定常物理问 题时，常常会遇到抛物型方程。其最简单的形式为一维热传导方程 0 ( 0) 2 2 = > ∂ ∂ − ∂ ∂ a x u a t u （5） 方程（5）可以有两种不同类型的定解问题： 初值问题（也称为 Cauchy 问题）
Ou a-u 0t>0.-∞0<x<+ (6) x02=o(x) 0<X<+ 初边值问题 Ou a2u =00<t<7,0<x<l l(x,0)=(x) (7) 0)=g(),u(,1)=g:(),0≤t≤7 其中(x),g1(1),g2()为已知函数,且满足连接条件 q(0)=g1(0),g()=g2(0) 问题(7)中的边界条件u(0,D)=g1(m),u(l,1)=g2(1)称为第一类边界条件。第二类和 第三类边界条件为 A1()=g;(0),01≤7 (8) +2(O)=g2(),0≤t≤T 其中A1(m)≥0,A2(1)≥0。当A1(1)=2(1)≡0时,为第二类边界条件,否则称为第三 类边界条件。 双曲型方程的最简单形式为一阶双曲型方程 +a (9) at ax 物理中常见的一维振动与波动问题可用二阶波动方程 a-u l (10) t 描述,它是双曲型方程的典型形式。方程(10)的初值问题为 a2u a2u t>0.-∞0<x<+0 l(x,0)=q(x) 0<x<+0 (11) a=(x) 0<X<+00 边界条件一般也有三类,最简单的初边值问题为
-241- ⎪ ⎩ ⎪ ⎨ ⎧ = − ∞ < < +∞ = > − ∞ < < +∞ ∂ ∂ − ∂ ∂ u x x x t x x u a t u ( ,0) ( ) 0 0, 2 2 ϕ （6） 初边值问题 ⎪ ⎪ ⎪ ⎩ ⎪ ⎪ ⎪ ⎨ ⎧ = = ≤ ≤ = = < < < < ∂ ∂ − ∂ ∂ u t g t u l t g t t T u x x t T x l x u a t u (0, ) ( ), ( , ) ( ), 0 ( ,0) ( ) 0 0 , 0 1 2 2 2 ϕ （7） 其中 ( ), ( ), ( ) 1 2 ϕ x g t g t 为已知函数，且满足连接条件 (0) (0), ( ) (0) 1 2 ϕ = g ϕ l = g 问题（7）中的边界条件 (0, ) ( ), ( , ) ( ) 1 2 u t = g t u l t = g t 称为第一类边界条件。第二类和 第三类边界条件为 t u g t t T x u t u g t t T x u x l x = ≤ ≤ ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ + ∂ ∂ = ≤ ≤ ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ − ∂ ∂ = = ( ) ( ), 0 ( ) ( ), 0 2 2 1 0 1 λ λ （8） 其中 ( ) 0, ( ) 0 λ1 t ≥ λ2 t ≥ 。当 ( ) ( ) 0 λ1 t = λ2 t ≡ 时，为第二类边界条件，否则称为第三 类边界条件。 双曲型方程的最简单形式为一阶双曲型方程 = 0 ∂ ∂ + ∂ ∂ x u a t u （9） 物理中常见的一维振动与波动问题可用二阶波动方程 2 2 2 2 2 x u a t u ∂ ∂ = ∂ ∂ （10） 描述，它是双曲型方程的典型形式。方程（10）的初值问题为 ⎪ ⎪ ⎪ ⎩ ⎪ ⎪ ⎪ ⎨ ⎧ = − ∞ < < +∞ ∂ ∂ = − ∞ < < +∞ > − ∞ < < +∞ ∂ ∂ = ∂ ∂ = x x t u u x x x t x x u a t u t ( ) ( ,0) ( ) 0, 0 2 2 2 2 2 φ ϕ （11） 边界条件一般也有三类，最简单的初边值问题为
auea a-u t>0.0<x<l (x,0)=(x) =x)0x≤1 l(0,1)=g1(1),t(l,t)=g2(1)0≤t≤T 如果偏微分方程定解问题的解存在,唯一且连续依赖于定解数据(即出现在方程 和定解条件中的已知函数),则此定解问题是适定的。可以证明,上面所举各种定解问 题都是适定的。 §2偏微分方程的差分解法 差分方法又称为有限差分方法或网格法,是求偏微分方程定解问题的数值解中应用 最广泛的方法之一。它的基本思想是:先对求解区域作网格剖分,将自变量的连续变化 区域用有限离散点(网格点)集代替;将问题中出现的连续变量的函数用定义在网格点 上离散变量的函数代替;通过用网格点上函数的差商代替导数,将含连续变量的偏微分 方程定解问题化成只含有限个未知数的代数方程组(称为差分格式)。如果差分格式有 解,且当网格无限变小时其解收敛于原微分方程定解问题的解,则差分格式的解就作为 原问题的近似解(数值解)。因此,用差分方法求偏微分方程定解问题一般需要解决以 下问题 (i)选取网格 (i)对微分方程及定解条件选择差分近似,列出差分格式; (ⅲi)求解差分格式 (ⅳ)讨论差分格式解对于微分方程解的收敛性及误差估计。 下面我们只对偏微分方程的差分解法作一简要的介绍 2.1椭圆型方程第一边值问题的差分解法 以 Poisson方程(1)为基本模型讨论第一边值问题的差分方法。 考虑 Poisson方程的第一边值问题(3) ax+a2=/(x)(x,y)∈ u(x, y)laper=p(x, y) I=aQ2 取hz分别为x方向和y方向的步长,以两族平行线x=xk=Mh,y=y,=jr (k,j=0,±1,±2,…)将定解区域剖分成矩形网格。节点的全体记为 R={(x,y)x=M,y,=Jx,,为整数}。定解区域内部的节点称为内点,记内点 集R∩g为n2。边界r与网格线的交点称为边界点,边界点全体记为F2。与节点 (xk,y)沿x方向或y方向只差一个步长的点(x,y)和(x,ym)称为节点 (xk,y)的相邻节点。如果一个内点的四个相邻节点均属于ΩUr,称为正则内点,正 则内点的全体记为Ω),至少有一个相邻节点不属于ΩUT的内点称为非正则内点, 非正则内点的全体记为92)。我们的问题是要求出问题(3)在全体内点上的数值解。 为简便记,记(k,)=(xk,y),(k,D=(xk,y),f=f(xk,y)。对正则内点
-242- ⎪ ⎪ ⎪ ⎩ ⎪ ⎪ ⎪ ⎨ ⎧ = = ≤ ≤ = ≤ ≤ ∂ ∂ = > < < ∂ ∂ = ∂ ∂ = u t g t u l t g t t T x x l t u u x x t x l x u a t u t (0, ) ( ), ( , ) ( ) 0 ( ,0) ( ), ( ) 0 0, 0 1 2 0 2 2 2 2 2 ϕ φ 如果偏微分方程定解问题的解存在，唯一且连续依赖于定解数据（即出现在方程 和定解条件中的已知函数），则此定解问题是适定的。可以证明，上面所举各种定解问 题都是适定的。 §2 偏微分方程的差分解法 差分方法又称为有限差分方法或网格法，是求偏微分方程定解问题的数值解中应用 最广泛的方法之一。它的基本思想是：先对求解区域作网格剖分，将自变量的连续变化 区域用有限离散点（网格点）集代替；将问题中出现的连续变量的函数用定义在网格点 上离散变量的函数代替；通过用网格点上函数的差商代替导数，将含连续变量的偏微分 方程定解问题化成只含有限个未知数的代数方程组（称为差分格式）。如果差分格式有 解，且当网格无限变小时其解收敛于原微分方程定解问题的解，则差分格式的解就作为 原问题的近似解（数值解）。因此，用差分方法求偏微分方程定解问题一般需要解决以 下问题： （i）选取网格； （ii）对微分方程及定解条件选择差分近似，列出差分格式； （iii）求解差分格式； （iv）讨论差分格式解对于微分方程解的收敛性及误差估计。 下面我们只对偏微分方程的差分解法作一简要的介绍。 2.1 椭圆型方程第一边值问题的差分解法 以 Poisson 方程（1）为基本模型讨论第一边值问题的差分方法。 考虑 Poisson 方程的第一边值问题（3） ⎪ ⎩ ⎪ ⎨ ⎧ = Γ = ∂Ω = ∈Ω ∂ ∂ + ∂ ∂ ∈Γ ( , ) | ( , ) ( , ) ( , ) ( , ) 2 2 2 2 u x y x y f x y x y y u x u x y ϕ 取 h,τ 分别为 x 方向和 y 方向的步长，以两族平行线 x x kh y y jτ = k = , = j = (k, j = 0,±1,±2,L) 将定解区域剖分成矩形网格。节点的全体记为 R {(x , y ) | x kh, y j ,i, j为整数} k j k j = = = τ 。定解区域内部的节点称为内点，记内点 集 R I Ω 为 Ωhτ 。边界Γ 与网格线的交点称为边界点，边界点全体记为Γhτ 。与节点 ( , ) k j x y 沿 x 方向或 y 方向只差一个步长的点 ( , ) k 1 j x y ± 和 ( , ) k j±1 x y 称为节点 ( , ) k j x y 的相邻节点。如果一个内点的四个相邻节点均属于Ω U Γ ，称为正则内点，正 则内点的全体记为 (1) Ω ，至少有一个相邻节点不属于Ω U Γ 的内点称为非正则内点， 非正则内点的全体记为 (2) Ω 。我们的问题是要求出问题（3）在全体内点上的数值解。 为简便记，记( , ) ( , ), ( , ) ( , ), ( , ) k j k j k , j k j k j = x y u k j = u x y f = f x y 。对正则内点
(k,j)∈Ω",由二阶中心差商公式 a(k+1D=2ak,D)+a(k-1.D+O() a2L (k,j+1)-2(k,j)+l(k,j-1) +O(r) Poisson方程(1)在点(k,处可表示为 l(k+1,j)-2u(k,)+u(k-1,),u(k,j+1)-2(k,j+l(k,j-1) h (12) 在式(12)中略去O(h2+x2),即得与方程(1)相近似的差分方程 lk 2uri +u +1-2lk,+l k,-=jk, (13) 式(13)中方程的个数等于正则内点的个数,而未知数k则除了包含正则内点处 解l的近似值,还包含一些非正则内点处u的近似值,因而方程个数少于未知数个数 在非正则内点处 Poisson方程的差分近似不能按式(13)给出,需要利用边界条件得到。 边界条件的处理可以有各种方案,下面介绍较简单的两种。 (i)直接转移 (i)线性插值 由式(13)所给出的差分格式称为五点菱形格式,实际计算时经常取h=r,此时 五点菱形格式可化为 )=f (14) 简记为 其中Ql,=4k 求解差分方程组最常用的方法是同步迭代法,同步迭代法是最简单的迭代方式。除 边界节点外,区域内节点的初始值是任意取定的。 例1用五点菱形格式求解 Laplace方程第一边值问题 0 (x,y)∈9 u(x,y)ls,er=1g[(+x)2+y2]r=02 其中Ω={(x,y)0≤x,y≤1}。取h=r 当h=r时,利用点(k,)(k±1,j-1),(k±1,j+1)构造的差分格式 2h2 (lk+1+1+k+y-1+lk-1+1+l-1y-1-4uk,)=fk -243
-243- (1) (k, j) ∈Ω ，由二阶中心差商公式 ( ) ( 1, ) 2 ( , ) ( 1, ) 2 2 ( , ) 2 2 O h h u k j u k j u k j x u k j + + − + − = ∂ ∂ ( ) ( , 1) 2 ( , ) ( , 1) 2 2 ( , ) 2 2 τ τ O u k j u k j u k j y u k j + + − + − = ∂ ∂ Poisson 方程（1）在点(k, j) 处可表示为 ( ) ( 1, ) 2 ( , ) ( 1, ) ( , 1) 2 ( , ) ( , 1) 2 2 , 2 2 τ τ = + + + − + − + + − + − f O h u k j u k j u k j h u k j u k j u k j k j （12） 在式（12）中略去 ( ) 2 2 O h +τ ，即得与方程（1）相近似的差分方程 k j k j k j k j k j k j k j f u u u h u u u 2 , , 1 , , 1 2 1, 2 , 1, 2 = − + + + − + − + − τ (13) 式（13）中方程的个数等于正则内点的个数，而未知数uk , j 则除了包含正则内点处 解u 的近似值，还包含一些非正则内点处u 的近似值，因而方程个数少于未知数个数。 在非正则内点处 Poisson 方程的差分近似不能按式（13）给出，需要利用边界条件得到。 边界条件的处理可以有各种方案，下面介绍较简单的两种。 (i) 直接转移 (ii) 线性插值 由式（13）所给出的差分格式称为五点菱形格式，实际计算时经常取 h = τ ，此时 五点菱形格式可化为 k j k j k j k j k j k j u u u u u f h2 1, 1, , 1 , 1 , , ( 4 ) 1 + + − + + + − − = (14) 简记为 k j k j u f h2 , , 1 ◊ = （15） 其中 uk , j = uk 1, j + uk 1, j + uk , j 1 + uk , j 1 − 4uk , j ◊ + − + − 。 求解差分方程组最常用的方法是同步迭代法，同步迭代法是最简单的迭代方式。除 边界节点外，区域内节点的初始值是任意取定的。 例 1 用五点菱形格式求解 Laplace 方程第一边值问题 ⎪ ⎩ ⎪ ⎨ ⎧ = + + Γ = ∂Ω = ∈Ω ∂ ∂ + ∂ ∂ ∈Γ ( , ) | lg[(1 ) ] 0 ( , ) 2 2 ( , ) 2 2 2 2 u x y x y x y y u x u x y 其中Ω = {(x, y) | 0 ≤ x, y ≤ 1}。取 3 1 h = τ = 。 当h = τ 时，利用点(k, j),(k ±1, j −1),(k ±1, j +1) 构造的差分格式 k j k j k j k j k j k j u u u u u f h2 1, 1 1, 1 1, 1 1, 1 , , ( 4 ) 2 1 + + + + − + − + + − − − = （16）
称为五点矩形格式,简记为 uk=fr (17) 2h 其中4=4k41+1+41-1+u4-1+l4-1=1-44, 22抛物型方程的差分解法 以一维热传导方程(5) l 为基本模型讨论适用于抛物型方程定解问题的几种差分格式。 首先对xt平面进行网格剖分。分别取h,r为x方向与t方向的步长,用两族平行直 x=x=Mh(k=0,±1+2,…),t=t=jr(=0,1,2,…),将x平面剖分成矩形网 格,节点为(x21)k=0+1+2,…,j=01,2,…)。为简便起见,记(k,)=(xk,y), (k,j)=(xk,y),0k=(x),g1=81(1),g2=g2(t,),41=41() 221微分方程的差分近似 在网格内点(k,)处,对分别采用向前、向后及中心差商公式,对2采用二 阶中心差商公式,一维热传导方程 可分别表示为 k,/+)-(k,D-aa(k+1.D)-2ak,)+(k-1.D=0(x+h hr n(k,)-a(k,/-1)-(k+1)-2a(k,)+a(k-1)=0(r+h) h a(k,)+1)-a,-1)-a2k+1.)-2a(k,D)+a(k-D=0 由此得到一维热传导方程的不同的差分近似 (18) 气。2-2,+u=0 (19) h (20) h 222初、边值条件的处理 为用差分方程求解定解问题(6),(7)等,还需对定解条件进行离散化。 对初始条件及第一类边界条件,可直接得到 l(x,0) (k=0,±1…或k=0,1,…,n) uo=u(0,t)=8, l(1,)=82
-244- 称为五点矩形格式，简记为 2 2 1 h ٱ k j k j u f , = , （17） 其中ٱuk , j = uk+1, j+1 + uk+1, j−1 + uk−1, j+1 + uk−1, j−1 − 4uk , j 。 2.2 抛物型方程的差分解法 以一维热传导方程（5） 0 ( 0) 2 2 = > ∂ ∂ − ∂ ∂ a t u a t u 为基本模型讨论适用于抛物型方程定解问题的几种差分格式。 首先对 xt 平面进行网格剖分。分别取 h,τ 为 x 方向与t 方向的步长，用两族平行直 线 x = x = kh(k = 0,±1,±2,L) k ，t = t = j ( j = 0,1,2,L) j τ ，将 xt 平面剖分成矩形网 格，节点为(x ,t )(k = 0,±1,±2,L, j = 0,1,2,L) k j 。为简便起见，记( , ) ( , ) k j k j = x y ， ( , ) ( , ) k j u k j = u x y ， ( ) k k ϕ = ϕ x ， ( ) 1 j 1 j g = g t ， ( ) 2 j 2 j g = g t ， ( ) 1 j 1 j λ = λ t ， ( ) 2 j 2 j λ = λ t 。 2.2.1 微分方程的差分近似 在网格内点(k, j) 处，对 t u ∂ ∂ 分别采用向前、向后及中心差商公式，对 2 2 x u ∂ ∂ 采用二 阶中心差商公式，一维热传导方程（5）可分别表示为 ( ) ( , 1) ( , ) ( 1, ) 2 ( , ) ( 1, ) 2 2 O h h u k j u k j u k j a u k j u k j = + + − + − − + − τ τ ( ) ( , ) ( , 1) ( 1, ) 2 ( , ) ( 1, ) 2 2 O h h u k j u k j u k j a u k j u k j = + + − + − − − − τ τ ( ) ( 1, ) 2 ( , ) ( 1, ) 2 ( , 1) ( , 1) 2 2 O h h u k j u k j u k j a u k j u k j = + + − + − − + − − τ τ 由此得到一维热传导方程的不同的差分近似 0 2 2 , 1 , 1, , 1, = − + − + − + − h u u u a uk j uk j k j k j k j τ （18） 0 2 2 , , 1 1, , 1, = − + − − − + − h u u u a uk j uk j k j k j k j τ （19） 0 2 2 2 , 1 , 1 1, , 1, = − + − + − − + − h u u u a uk j uk j k j k j k j τ （20） 2.2.2 初、边值条件的处理 为用差分方程求解定解问题（6），（7）等，还需对定解条件进行离散化。 对初始条件及第一类边界条件，可直接得到 k k k u ,0 = u(x ,0) = ϕ (k = 0,±1,L或k = 0,1,L,n) （21） n j j j j j ij u u l t g u u t g , 2 0, ( , ) (0, ) = = = = ( j = 0,1,L,m −1) （22）