第二十章偏微分方程的数值解 自然科学与工程技术中种种运动发展过程与平衡现象各自遵守一定的规律。这些规 律的定量表述一般地呈现为关于含有未知函数及其导数的方程。我们将只含有未知多元 函数及其偏导数的方程,称之为偏微分方程。 方程中出现的未知函数偏导数的最高阶数称为偏微分方程的阶。如果方程中对于未 知函数和它的所有偏导数都是线性的,这样的方程称为线性偏微分方程,否则称它为非 线性偏微分方程 初始条件和边界条件称为定解条件,未附加定解条件的偏微分方程称为泛定方程 对于一个具体的问题,定解条件与泛定方程总是同时提出。定解条件与泛定方程作为 个整体,称为定解问题 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.-∞00.-∞0<x<+0 l(x,0)=q(x) 0<x<+0 (11) a=(x) 0<X<+00 边界条件一般也有三类,最简单的初边值问题为
-241- ⎪ ⎩ ⎪ ⎨ ⎧ = − ∞ − ∞ − ∞ < < +∞ ∂ ∂ = ∂ ∂ = 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）
T 其中n 对第二、三类边界条件则需用差商近似。下面介绍两种较简单的处理方法。 (i)在左边界(x=0)处用向前差商近似偏导数一,在右边界(x=D)处用向后差 商近似偏导数一,即 3,s-(0,0+O(h) (n,j)-l(n-1,j) +O(h) (n,j) h 即得边界条件(8)的差分近似为 h ,40y= (=0,…,m) g (i)用中心差商近似一,即 l(1,j)-l(-1,j) O(h2) 2h (=0,…,m) l(n+1,j)-l(n-1,) +O(h2) n,) 2h 则得边界条件的差分近似为 l---10y=8 2h (=0,1,…,m) 2h 这样处理边界条件,误差的阶数提高了,但式(24)中出现定解区域外的节点(-1,j)和 (n+1,j),这就需要将解拓展到定解区域外。可以通过用内节点上的a值插值求出u-1 和n1y,也可以假定热传导方程(5)在边界上也成立,将差分方程扩展到边界节点 上,由此消去u1,和n+1 223几种常用的差分格式 下面我们以热传导方程的初边值问题(7)为例给出几种常用的差分格式 (i)古典显式格式 为便于计算,令"/式(18)改写成以下形式
-245- 其中 τ T m h l n = , = 。 对第二、三类边界条件则需用差商近似。下面介绍两种较简单的处理方法。 （i）在左边界(x = 0) 处用向前差商近似偏导数 x u ∂ ∂ ，在右边界(x = l) 处用向后差 商近似偏导数 x u ∂ ∂ ，即 ( ) ( , ) ( 1, ) ( ) (1, ) (0, ) ( , ) (0, ) O h h u n j u n j x u O h h u j u j x u n j j + − − = ∂ ∂ + − = ∂ ∂ ( j = 0,1,L,m) 即得边界条件（8）的差分近似为 ⎪ ⎪ ⎩ ⎪ ⎪ ⎨ ⎧ + = − − = − − j n j j n j n j j j j j j u g h u u u g h u u 2 , 2 , 1, 1 0, 1 1, 0, λ λ ( j = 0,1,L,m) （23） （ii）用中心差商近似 x u ∂ ∂ ，即 ( ) 2 ( 1, ) ( 1, ) ( ) 2 (1, ) ( 1, ) 2 ( , ) 2 (0, ) O h h u n j u n j x u O h h u j u j x u n j j + + − − = ∂ ∂ + − − = ∂ ∂ ( j = 0,1,L,m) 则得边界条件的差分近似为 ⎪ ⎪ ⎩ ⎪ ⎪ ⎨ ⎧ + = − − = − + − − j n j j n j n j j j j j j u g h u u u g h u u 2 , 2 1, 1, 1 0, 1 1, 1, 2 2 λ λ ( j = 0,1,L,m) （24） 这样处理边界条件，误差的阶数提高了，但式（24）中出现定解区域外的节点(−1, j) 和 (n +1, j) ，这就需要将解拓展到定解区域外。可以通过用内节点上的u 值插值求出u−1, j 和un+1, j ，也可以假定热传导方程（5）在边界上也成立，将差分方程扩展到边界节点 上，由此消去u−1, j 和un+1, j 。 2.2.3 几种常用的差分格式 下面我们以热传导方程的初边值问题（7）为例给出几种常用的差分格式。 (i) 古典显式格式 为便于计算，令 2 h a r τ = ，式（18）改写成以下形式
lkA=41,+(1-2r)uk,+Plk-1 将式(18)与(21),(22)结合,我们得到求解问题(7)的一种差分格式 =Pk41+(1-2r)uk +ru (k=1,2…,n-1,j=0,1,…,m-1) (k=1,2,…,n) 10o1=g1;ln=8: (=1,2 由于第0层(=0)上节点处的a值已知(u10=k),由式(25)即可算出u在第一层 (=1)上节点处的近似值uk。重复使用式(25),可以逐层计算出各层节点的近似值。 (ⅱi)古典隐式格式 将(19)整理并与式(21),(22)联立,得差分格式如下 e lkj +r(u k-1,j+1 1)(k=1,2,…,n-1,j=0,1,…,m-1) Pu (k=0,…,n) 10=g1 g (=01,…,m) aT Bh。虽然第0层上的a值仍为已知,但不能由式(30)直接计算以上各层节 点上的值uk,,故差分格式(26)称为古典隐式格式。 (ⅲi)杜福特一弗兰克尔( Do Fort- Frankel)格式 DoFort- Frankel格式是三层显式格式,它是由式(24)与(25),(26)结合得到 的。具体形式如下 lyk1(k=1,2,…,n-1,=1,2 1) 1+2 (k=0,1,…,n) uo,=g1, un,j=2j (=0.1,…,m) 用这种格式求解时,除了第0层上的值uk由初值条件(21)得到,必须先用二层格式 求出第1层上的值k1,然后再按格式(27)逐层计算,(=2,3,…,m) 23双曲型方程的差分解法 对二阶波动方程(10) l =a 如果令v1 ,则方程(10)可化成一阶线性双曲型方程组 at- ax (28) 记v=(v1,n2),则方程组(28)可表成矩阵形式
-246- k j k j k j k j u , 1 ru 1, r u , ru 1, (1 2 ) + = + + − + − 将式（18）与（21），（22）结合，我们得到求解问题（7）的一种差分格式 ⎪ ⎩ ⎪ ⎨ ⎧ = = = = = + = + + − + − = − = − , ( 1,2, , ) ( 1,2, , ) (1 2 ) ( 1,2, , 1, 0,1, , 1) 0, 1 , 2 ,0 , 1 1, , 1, u g u g j m u k n u ru r u ru k n j m j j n j j k k k j k j k j k j L L L L ϕ （25） 由于第 0 层( j = 0) 上节点处的u 值已知( ) uk ,0 = ϕ k ，由式（25）即可算出u 在第一层 ( j = 1) 上节点处的近似值uk ,1 。重复使用式（25），可以逐层计算出各层节点的近似值。 （ii）古典隐式格式 将（19）整理并与式（21），（22）联立，得差分格式如下 ⎪ ⎩ ⎪ ⎨ ⎧ = = = = = + = + + + − + + − + = − = − , ( 0,1, , ) ( 0,1, , ) ( 2 )( 1,2, , 1, 0,1, , 1) 0, 1 , 2 ,0 , 1 , 1, 1 , 1 1, 1 u g u g j m u k n u u r u u u k n j m j j n j j k k k j k j k j k j k j L L L L ϕ （26） 其中 2 h a r τ = 。虽然第 0 层上的u 值仍为已知，但不能由式（30）直接计算以上各层节 点上的值uk , j ，故差分格式（26）称为古典隐式格式。 （iii）杜福特—弗兰克尔（DoFort—Frankel）格式 DoFort—Frankel 格式是三层显式格式，它是由式（24）与（25），（26）结合得到 的。具体形式如下： ⎪ ⎪ ⎩ ⎪ ⎪ ⎨ ⎧ = = = = = = − = − + − + + + + = + − − , ( 0,1, , ) ( 0,1, , ) ( 1,2, , 1, 1,2, , 1) 1 2 1 2 ( ) 1 2 2 0, 1 , 2 ,0 , 1 1, 1, , 1 u g u g j m u k n u k n j m r r u u r r u j j n j j k k k j k j k j k j L L L L ϕ (27) 用这种格式求解时，除了第 0 层上的值uk ,0 由初值条件（21）得到，必须先用二层格式 求出第 1 层上的值uk ,1 ，然后再按格式（27）逐层计算 ( 2,3, , ) uk , j j = L m 。 2.3 双曲型方程的差分解法 对二阶波动方程（10） 2 2 2 2 2 x u a t u ∂ ∂ = ∂ ∂ 如果令 x u v t u v ∂ ∂ = ∂ ∂ 1 = 2 , ，则方程（10）可化成一阶线性双曲型方程组 ⎪ ⎪ ⎩ ⎪ ⎪ ⎨ ⎧ ∂ ∂ = ∂ ∂ ∂ ∂ = ∂ ∂ x v t v x v a t v 2 1 1 2 2 （28） 记 T v (v ,v ) = 1 2 ，则方程组（28）可表成矩阵形式
阵A有两个不同的特征值A=±a,故存在非奇异矩阵P,使得 PAP- A 作变换w=Pv=(,2),方程组(29)可化成 (30) 方程组(30)由两个独立的一阶双曲型方程联立而成。因此下面主要讨论一阶双曲型方 程的差分解法。 考虑一阶双曲型方程的初值问题 0 t>0 (31) l(x,0)=(x) 000,则a必等于b,也就是说其具有圆柱或球体的对称关系。同时, 式中f(x,1,u,u/ax)一项为流通量(fux),而S(x,t,u,au/ax)为来源( source))项。 c(x,Lu,a/ax)为偏微分方程的对角线系数矩阵。若某一对角线元素为0,则表示该偏 微分方程为椭圆型偏微分方程,若为正值(不为0),则为拋物型偏微分方程。请注意c的 对角线元素一定不全为0。偏微分方程初始值可表示为
-247- x v A x a v t v ∂ ∂ = ∂ ∂ ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ = ∂ ∂ 1 0 0 2 （29） 矩阵 A 有两个不同的特征值λ = ±a ，故存在非奇异矩阵 P ，使得 ⎥ = Λ ⎦ ⎤ ⎢ ⎣ ⎡ − = − a a PAP 0 0 1 作变换 T w Pv (w ,w ) = = 1 2 ，方程组（29）可化成 x w t w ∂ ∂ = Λ ∂ ∂ （30） 方程组（30）由两个独立的一阶双曲型方程联立而成。因此下面主要讨论一阶双曲型方 程的差分解法。 考虑一阶双曲型方程的初值问题 ⎪ ⎩ ⎪ ⎨ ⎧ = − ∞ − ∞ 0 ，则a 必等于b ，也就是说其具有圆柱或球体的对称关系。同时， 式中 f (x,t,u,∂u ∂x) 一项为流通量(flux)，而 s(x,t,u,∂u ∂x) 为来源(source)项。 c(x,t,u,∂u ∂x) 为偏微分方程的对角线系数矩阵。若某一对角线元素为 0，则表示该偏 微分方程为椭圆型偏微分方程，若为正值(不为 0)，则为拋物型偏微分方程。请注意c 的 对角线元素一定不全为 0。偏微分方程初始值可表示为
l(x,0)=v0(x) (36) 而边界条件为 p(x,t,)+q(x,1)f(x,1,u,ou/ax)=0 其中x为两端点位置,即a或b 用以解含上述初始值及边界值条件的偏微分方程的 MATLAB命令 pdepe的用法如 sol= pdepe(m, pdepe, icfun, bcfun, xmesh, tspan, options) 其中 m为问题之对称参数 xmesh为空间变量x的网格点(mesh)位置向量,即 xmesh=[x0,x1…,xx],其中 x0=a(起点),xN=b(终点) Ispan为时间变量t的向量,即spon=[,12…,lM],其中b为起始时间,t为 终点时间。 pdlm为使用者提供的pde函数文件。其函数格式如下 Ic,f, s]=pdefun(x, t, u, dudx) 亦即,使用者仅需提供偏微分方程中的系数向量。c,∫和s均为行( column)向量,而 向量c即为矩阵c的对角线元素。 ic/im提供解l的起始值,其格式为u= icfun(x)值得注意的是u为行向量。 bc/m使用者提供的边界条件函数,格式如下 pl, ql, pr, q]=bcfun(xl,ul, xr, ur;, 1) 其中和r分别表示左边界(xl=b)和右边界(xr=a)u的近似解。输出变量中,pl 和ql分别表示左边界p和q的行向量,而p和q则为右边界P和q的行向量。 sol为解答输出。sol为多维的输出向量,sol(;:i)为u1的输出,即u1=sol(;,) 元素1(,k)=sol(,k,1)表示在t=1spom()和x= xmesh(k)时u1之答案。 options为求解器的相关解法参数。详细说明参见 odeset的使用方法。 1. MATLAB PDE求解器 pdepe的算法,主要是将原来的椭圆型和拋物线型偏微分 方程转化为一组常微分方程。此转换的过程是基于使用者所指定的mesh点,以二阶空 间离散化( spatial discretization)技术为之 Keel and berzins,1990),然后以ode5的指令 求解。采用odel5s的ode解法,主要是因为在离散化的过程中,椭圆型偏微分方程被 转化为一组代数方程,而拋物线型的偏微分方程则被转化为一组联立的微分方程。因而 原偏微分方程被离散化后,变成一组同时伴有微分方程与代数方程的微分代数方程组 故以ode5s便可顺利求解。 2.x的取点(mesh)位置对解的精确度影响很大,若 pdepe求解器给出“…has difficulty finding consistent initial considition”的讯息时,使用者可进一步将mesh点取密 点,即增加mesh点数。另外,若状态u在某些特定点上有较快速的变动时,亦需将 此处的点取密集些,以增加精确度。值得注意的是 pdepe并不会自动做 xmesh的自动取 点,使用者必须观察解的特性,自行作取点的操作。一般而言,所取的点数至少需大于 3以上。 3. tspan的选取主要是基于使用者对那些特定时间的状态有兴趣而选定。而间距
-248- ( , ) ( ) 0 0 u x t = v x （36） 而边界条件为 p(x,t,u) + q(x,t) f (x,t,u,∂u ∂x) = 0 （37） 其中 x 为两端点位置，即a 或b 用以解含上述初始值及边界值条件的偏微分方程的 MATLAB 命令 pdepe 的用法如 下： sol = pdepe(m, pdepe,icfun,bcfun, xmesh,tspan,options) 其中 m 为问题之对称参数； xmesh为空间变量 x 的网格点(mesh)位置向量，即 [ , , , ] 0 1 N xmesh = x x L x ，其中 x0 = a (起点)， xN = b (终点)。 tspan 为时间变量t 的向量，即 [ , , , ] 0 1 M tspan = t t L t ，其中 0t 为起始时间， Mt 为 终点时间。 pdefun 为使用者提供的 pde 函数文件。其函数格式如下： [ ] c, f ,s = pdefun(x,t,u, dudx) 亦即，使用者仅需提供偏微分方程中的系数向量。c ， f 和 s 均为行(column)向量，而 向量c 即为矩阵c 的对角线元素。 icfun 提供解u 的起始值，其格式为u = icfun(x) 值得注意的是u 为行向量。 bcfun 使用者提供的边界条件函数，格式如下： [ ] pl, ql, pr, ql = bcfun(xl,ul, xr,ur,t) 其中ul 和ur 分别表示左边界(xl = b)和右边界(xr = a) u 的近似解。输出变量中，pl 和 ql 分别表示左边界 p 和 q 的行向量，而 pr 和 qr 则为右边界 p 和 q 的行向量。 sol 为解答输出。 sol 为多维的输出向量， sol(:,: i) 为ui 的输出，即u sol(:,:,i) i = 。 元素u ( j, k) sol( j, k,i) i = 表示在t = tspan( j) 和 x = xmesh(k) 时ui 之答案。 options 为求解器的相关解法参数。详细说明参见 odeset 的使用方法。 注： 1. MATLAB PDE 求解器 pdepe 的算法，主要是将原来的椭圆型和拋物线型偏微分 方程转化为一组常微分方程。此转换的过程是基于使用者所指定的 mesh 点，以二阶空 间离散化(spatial discretization)技术为之(Keel and Berzins,1990)，然后以 ode15s 的指令 求解。采用 ode15s 的 ode 解法，主要是因为在离散化的过程中，椭圆型偏微分方程被 转化为一组代数方程，而拋物线型的偏微分方程则被转化为一组联立的微分方程。因而， 原偏微分方程被离散化后，变成一组同时伴有微分方程与代数方程的微分代数方程组， 故以 ode15s 便可顺利求解。 2. x 的取点(mesh)位置对解的精确度影响很大，若 pdepe 求解器给出“…has difficulty finding consistent initial considition”的讯息时，使用者可进一步将 mesh 点取密 一点，即增加 mesh 点数。另外，若状态u 在某些特定点上有较快速的变动时，亦需将 此处的点取密集些，以增加精确度。值得注意的是 pdepe 并不会自动做 xmesh 的自动取 点，使用者必须观察解的特性，自行作取点的操作。一般而言，所取的点数至少需大于 3 以上。 3. tspan 的选取主要是基于使用者对那些特定时间的状态有兴趣而选定。而间距
( step size)的控制由程序自动完成 4.若要获得特定位置及时间下的解,可配合以 deval命令。使用格式如下: Luout, duoutdx]= deval(m, xmesh, ui, xout) 其中 m代表问题的对称性。m=0表示平板;m=1表示圆柱体;m=2表示球体。其意 义同 pdepe中的自变量m。 xme$h为使用者在pepe中所指定的输出点位置向量。 xmesh=[x0,x1…,x] u即sol(;)。也就是说其为 pdepe输出中第i个输出a在各点位置 xmesh处, 时间固定为1=1spun()下的解。 xolt为所欲内插输出点位置向量。此为使用者重新指定的位置向量 为基于所指定位置xout,固定时间t下的相对应输出 duoutdx为相对应的 du dx输出值 ref. Keel, R D. and M. Berzins, "A Method for the Spatial Discritization of Parabolic Equations in One Space Variable", SIAM J. Sci. and Sat Comput., Vol. ll, pp 1-32, 1990 以下将以数个例子,详细说明 pdepe的用法 3.2求解一维偏微分方程 例2试解以下之偏微分方程式 a2u 其中0≤x≤1,且满足以下之条件限制式 (i)起始值条件 IC:(x,0)=sin(丌x) (i)边界条件 BC1:l(0,t)=0 BC2: Te- a(,) =0 注:本问题的解析解为(x,1)=e-sin(m) 解下面将叙述求解的步骤与过程。当完成以下各步骤后,可进一步将其汇总为 主程序ex201m,然后求解 步骤1将欲求解的偏微分方程改写成如式的标准式。 此即 c(x,t, u, u/ax) f(x, t, u, au/ax) s(x,t, u, Ou/ ax)=0 和m=0 步骤2编写偏微分方程的系数向量函数。 function [c, f, s]=ex20 lpdefun(x, t, u, dudx)
-249- (step size)的控制由程序自动完成。 4. 若要获得特定位置及时间下的解，可配合以 pdeval 命令。使用格式如下： [ ] uout, duoutdx = pdeval(m, xmesh,ui, xout) 其中 m 代表问题的对称性。m =0 表示平板；m =1 表示圆柱体；m =2 表示球体。其意 义同 pdepe 中的自变量m 。 xmesh为使用者在 pdepe 中所指定的输出点位置向量。 0 1 [,, , ] N xmesh x x x = L 。 ui 即 sol( j,:,i) 。也就是说其为 pdepe 输出中第 i 个输出ui 在各点位置 xmesh 处， 时间固定为t tspan( j) j = 下的解。 xout 为所欲内插输出点位置向量。此为使用者重新指定的位置向量。 uout 为基于所指定位置 xout ，固定时间 f t 下的相对应输出。 duoutdx 为相对应的 du dx 输出值。 ref. Keel,R.D. and M. Berzins,“A Method for the Spatial Discritization of Parabolic Equations in One Space Variable”，SIAM J. Sci. and Sat. Comput.,Vol.11,pp.1-32,1990. 以下将以数个例子，详细说明 pdepe 的用法。 3.2 求解一维偏微分方程 例 2 试解以下之偏微分方程式 2 2 2 x u t u ∂ ∂ = ∂ ∂ π 其中0 ≤ x ≤ 1，且满足以下之条件限制式 (i)起始值条件 IC：u(x,0) = sin(π x) (ii)边界条件 BC1：u(0,t) = 0 BC2： 0 (1, ) = ∂ ∂ + − x u t e t π 注：本问题的解析解为u(x,t) e sin( x) t π − = 解 下面将叙述求解的步骤与过程。当完成以下各步骤后，可进一步将其汇总为一 主程序 ex20_1.m，然后求解。 步骤 1 将欲求解的偏微分方程改写成如式的标准式。 0 2 0 0 ⎟ + ⎠ ⎞ ⎜ ⎝ ⎛ ∂ ∂ ∂ ∂ = ∂ ∂ x u x x x t u π 此即 ( ) 2 c x,t,u,∂u ∂x = π ( ) x u f x t u u x ∂ ∂ , , ,∂ ∂ = s( ) x,t,u,∂u ∂x = 0 和m = 0 。 步骤 2 编写偏微分方程的系数向量函数。 function [c,f,s]=ex20_1pdefun(x,t,u,dudx)