第二十章偏微分方程的数值解 自然科学与工程技术中种种运动发展过程与平衡现象各自遵守一定的规律。这些规 律的定量表述一般地呈现为关于含有未知函数及其导数的方程。我们将只含有未知多元 函数及其偏导数的方程，称之 为微分万程 数称为偏微分方程的阶如果方程中对于术 导数都是线性的，这样的方程称为线性佩微分方程，否则称 和的 个整体，称为定解问题。 §1偏微分方程的定解问题 的定常（即不随时间变化）过程，都可用椭圆型方程来描述。其最典 型、最简单的形式是泊松(Poisson)方程 (1) 特别地，当f(x,y)■0时，即为拉普拉斯Laplace)方程，又称为调和方程 Au=O'u 0u 42+ (2) 带有稳定热源或内部无热源的稳定温度场的温度分布，不可压缩流体的稳定无旋流动及 静电场的电势等均满足这类方程。 Poisson方程的第一边值问题为 au+02=f0,) x2+0y2 (x,y)∈2 (3) ux,月lxr=p(x)「=2 其中2为以「为边界的有界区域，「为分段光滑曲线，2U厂称为定解区域， f(x,y),(x,y)分别为2，「上的已知连续函数。 第二类和第三类边界条件可统一表示成 (+m=o, (4) 其中n为边界「的外法线方向。当a=0时为第二类边界条件，a≠0时为第三类边界 条件。 在研究热传导过程 ,气体扩散现象及电础 题时，常常会遇 -=0a>0 (5) 方程(5)可以有两种不同类型的定解问题： 初值问题（也称为Cauchy问题） -240
-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 问题）
1>0,-0 u(0,1)=81(),0,)=g2t),0≤1≤T 其中(x,g1(),g2()为已知函数，且满足连接条件 p(0)=g,(0),pI0=g2(0) 问题(7)中的边界条件(0，)=g1(),)=g2)称为第一类边界条件。第二类和 第三类边界条件为 La-4ou =gt),0≤1≤T (8) =g,(0,0≤1≤T 其中)20,20≥0。当()=2()=0时，为第二类边界条件，否则称为第三 类边界条件 双曲型方程的最简单形式为一阶双曲型方程 (9) 物理中常见的一维振动与波动问题可用二阶波动方程 - (10) 描述，它是双曲型方程的典型形式。方程(10)的初值问题为 t>0.-0<x<+a (x,0)=p(x) -00<Y<+0 (11 =(x) -00<x<+0 边界条件一般也有三类，最简单的初边值问题为 -241
-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） 边界条件一般也有三类，最简单的初边值问题为
=a8'u [ou ? t>0,0<x<1 u(x,0)=p(x), =x)0≤x≤1 u(0,)=g1(),l,)=g20)0≤t≤T 加如果偏微分方程定解司顺的解存在，一日连续依干定解数据（即出现在方程 和定解条件中的已知函数)，则此定解问题是适定的。可以证明，上面所举各种定解问 题都是适定的。 $2偏微分方程的差分解 又称为有限 法或网格法 是微分方程定解向题的影 值解中应用 离散变 离散点的 )集 出现的 方程定解只含有限 的代数 方程组称为差分格 P) 加里差分热 解，且当网格无限变小时其解收敛于原微分方程定解问题的解，则差分格式的解就作为 原问题的近似解（数值解）。因此，用差分方法求偏微分方程定解问题一般需要解决以 下问题： ()选取网格： ()对微分方程及定解条件选择差分近似，列出差分格式： ()求解 下面我们讨论 分格 对于微分方程 的收 1 椭圆型方程第 估问题 差分解 SO方程(1)为基本模型讨论第一边值问题的差分方法。 考虑Poisson方程的第一边值问题(3) ax2+=f》Gn (x,y)lr=p(xy)「=2 取h,t分别为x方向和y方向的步长，以两族平行线x=x=k仇y=y,=江 (k,j=0,士1,2，)将定解区域剖分成矩形网格。节点的全体记为 R={(xy,)川x=仇，y,=j，为整数}。定解区域内部的节点称为内点，记内点 集R∩2为2，。边界T与网格线的交点称为边界点，边界点全体记为「：。与节点 (x,y)沿x方向或y方向只差一个步长的点(xy,)和(x,y,)称为节点 (x,y,)的相邻节点。如果一个内点的四个相邻节点均屈于2U厂，称为正则内点，正 则内点的全体记为2"，至少有一个相邻节点不属于2U厂的内点称为非正则内点， 非正则内点的全体记为2。我们的问题是要求出问题(3)在全体内点上的数值解 为简便记，记(k,)=(x,y,)4k,)=(x,y小，=fxy)。对正则内点 242
-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,)e2四，由二阶中心差商公式 8u -k+L》-2k》+k-LD+0) O'u _kj+)-2k》+k-D+0r) Poisson方程(1)在点(k,)处可表示为 k+1)-2k,)+4(k-1),ukj+)-2u(k)+(k,j-) 2 (12 =人，+0h2+x2) 在式(12)中略去Oh2+x2),即得与方程(1)相近似的差分方程 4-24+4山+4-2+w=人 13) 式(13)中方程的个数等于正则内点的个数，而未知数4，则除了包含正则内点处 解！的近似值，还包含一些非正则内点处山的近似值，因而方程个数少于未知数个数 在非正则内点处Poisson方程的差分近似不能按式(I3)给出，需要利用边界条件得到。 边界 件的处理可以有各种方案，下面介绍较简单的两种。 ①直接 回线 3)所给出的差分格式称为五点菱形格式，实际计算时经常取=t,此时 五点麦格式可化穷 u+-+m+n-4u)=f (14) 简记为 京0= (15) 其中0u,= -4u. 迭代法，同步迭代法是最简单的选代方式。除 域内节点的 意第边值问恩 (x.y)∈2 u(x.y)kver=lg[(1+x)2+y2]r=a0 其中2={(x,川0≤xy≤1}。取h=t= 当h=x时，利用点（化，），（化±1，了-1），(k±1，j+1)构造的差分格式 2+w+-wm+4-A-4)=f (16) -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）
称为五点矩形格式，简记为 元 (17) 其中i=4 挺勒型方程的差分解法一、 4 22 是-0a>0 为基本模型讨论适用于抛物型方程定解问题的几种差分格式。 首先对x1平面进行网格剖分。分别取h,x为x方向与1方向的步长，用两族平行直 线x=x=h(k=0,士12，)，1=1，=jrU=0,12,),将x1平面别分成矩形网 格，节点为(x,k=0,士1，士2，.，j=0,1,2,).为简便起见，记(k,)=(y,) (k,)=(x4,y,),p4=(x),8v=81(0,）,82,=82()，元，=元(0，)， 2.2.1微分方程的差分近似 在网格内点（化，》处，对分别采用向前、向后及中心差商公式，对 2采用二 阶中心差商公式， 一维热传导方程(5)可分别表示为 k-k_0k+-2k+k-LD.0+） k,》-k-》-ak+LD-2k刀+k-山D=0+) kj+)-kj-D-a4k+L》-2k》+k-LD=Or+的) 由此得到一维热传导方程的不同的差分近似 州-M-a"山-2w+山=0 (18) h (19) 4-a-2+-山=0 (20) h2 2.2.2初、边值条件的处理 为用差分方程求解定解问题(6)，(7)等，还需对定解条件进行离散化。 对初始条件及第一类边界条件，可直接得到 k0=xk,0)=pg(k=0,士1.或k=01,.,m） (21) 4oy=0,l,)=g U=0,1.,m-1) (22) 4n/=l,）=g2 244-
-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）
m=T 其中n= 对第二、三类边界条件则需用差商近似。下面介绍两种较简单的处理方法， ①在左边界仅=0)处用向前差商近似偏导数，在右边界（仁=D处用向后差 商近似偏导数，即 ou ul.j)-u(0.D+O(h) (U=0,1.,m) x 1 即得边界条件(8)的差分近似为 4-L-40=8 h =0,1.,m) (23) -山+y=82 h (im用中心差商近似 ,即 _LD-D+0) 2h n+-u(n-1D+O() U=0,l.,m) 2h 则得边界条件的差分近似为 4y-u-2y4w=8 uw-+,=82 j=0,l.,m) (24 2h 这样处理边界条件，误差的阶数提高了，但式(24)中出现定解区域外的节点(-1，)和 (n+1,),这就需要将解拓展到定解区域外。可以通过用内节点上的u值插值求出W 和山，也可以假定热传导方程(5)在边界上也成立，将差分方程扩展到边界节点 上，由此消去和un 223几种常用的差分格式 面我们 0古典显式格式 方程的初边值问题(7)为例给出几种常用的差分格式。 为便于计算。令r ,式(18)改写成以下形式 -245
-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）改写成以下形式
4,=ru+(1-2r)u+ug- 将式(18)与(21)，(22)结合，我们得到求解问题(7)的一种差分格式 [y4=u4y+(0-2r)4y+uk-y(k=1,2,.,n-1j=0,l.,m-) 40=pg (k=1,2,.,n) (25) (j=1,2.,m) 由于第0层(U=0)上节点处的u值已知(u0=p,),由式(25)即可算出u在第一层 U=1)上节点处的近似值·重复使用式(25)，可以逐层计算出各层节点的近似值。 (i)古典隐式格式 将(19)整理并与式(21)，(22)联立，得差分格式如下 [4m=4+r(uk1-21+u4-Hk=12,.,n-1j=01,.,m-l) .=Pu (k=01,.,n) (26) 40=81,山nJ=82) U=0,1,.,m) 其中r=。虽然第0层上的“值仍为已知，但不能由式(30)直接计算以让上各层节 点上的值，故差分格式(26)称为古典隐式格式。 ()杜特一弗卓尔(DoFort -Frankel)格式 -Frankel格式是三层显式格式，它是由式(24)与(25)，(26)结合符到 的。具体形式如下： mi+a=12-2m- 2. 44.0=pg (k=0L,.,) (27) 40y=81,4my=82 U=0,1.,m) 用这种格式求解时，除了第0层上的值山0由初值条件(21)得到，必须先用二层格式 求出第1层上的值1，然后再按格式(27)逐层计算4(=2,3，.，m)。 2.3双曲型方程的差分解法 对二阶波动方程(10) 如果 ，则方程(10)可化成一阶线性双曲型方程组 [=a (28) 记v=(,2)，则方程组(28)可表成矩阵形式 -246-
-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）可表成矩阵形式
会]度会 (29) 矩阵A有两个不同的特征值元=土，故存在非奇异矩阵P,使得 PAP-Ta0 Lo-a=A 作变换w=PY=(,w2)》',方程组(29)可化成 (30) 方程组(30)由两个独立的一阶双曲型方程联立而成。因此下面主要讨论一阶双曲型方 程的差分解法。 考虑一阶双曲型方程的初值问题 a+a=0 1>0-0<x<+00 (31) tx,0)=x) -0<X<+0 与抛物型方程的讨论类似，仍将x1平面剖分成矩形网格。取x方向步长为h,1方向步 长为t,网格线为x=x=kk=0,1,2,以，1=1，=jj=0,12,)。为简便起 见，记(k,)=(xk,y),(k,)=(x,y,),P=(x)。 以不同的差商近似偏导数，可以得到方程(9)的不同的差分近似 悲l-%L+a-=0 (32) h -a-e=0 h (33) 4-u-g+一4-=0 (34) 结合离散化的初始条件，可以得到几种简单的差分格式。 S3,】维秋态的微分方程的MAAB解法 MATLAB提供了一个指令pcpe,用以解以下的PDE方程式 0g-”u尝》+0 (35) 其中时间介于1。≤1≤t,之间，而位置x则介于[a,b]有限区域之间。m值表示问题的 对称性，其可 b) 因而，如果 0,则a必 也就是说其具有圆柱或球体 式中f(x,l,“,du/ax)一项为流通量(ux),而s(x,l,山，du/ax)为来源((source)项 c(x,l,u,ux)为偏微分方程的对角线系数矩阵。若某一对角线元素为0，则表示该偏 微分方程为椭圆型偏微分方程，若为正值（不为0），则为抛物型偏微分方程。请注意c的 对角线元素一定不全为0。偏微分方程初始值可表示为 -247
-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。偏微分方程初始值可表示为
u(x.1)=v(x) (36) 而边界条件为 x,40+g(x,0f(x,4,0u/ax)=0 (37) 下： sol=pdepe(m,pdepe,icfun,bcfun,xmesh,tspan,options) 其中 m为问题之对称参数 xmesh为空间变量x的网格点(mesh)位置向量，即xmesh=[x,x,xw],其中 x0=a(起点)，w=b(终点). span为时间变量1的向量，即tspan=o,4,.,lw],其中1。为起始时间，1w为 终点时间。 pdefun为使用者提供的pde函数文件。其函数格式如下： 即，使用者仅需提供偏微，可=pdefu(x“,dhu本 方程中的系数向量。c,∫和s均为行(column)向量，而 向量c即为矩阵c的对角线元素。 icim提供解4的起始值，其格式为u=icfx)值得注意的是u为行问量。 bcn使用者提供的边界条件函数，格式如下： [pl.ql.pr.ql]-bcfim(xl.u.xr.u.t) 其中l和w分别表 示左边界(xl=b)和右边界(xr=a)u的近似解。输出变量中，pl 和g分别表示左边界p和q的行向量，而pr和qr则为右边界p和q的行向量。 sol为解答输出。sol为多维的输出向量，sol(,:)为4，的输出，即山，=sol(::,) 元素山，(j,k)=sol(j,k,i)表示在t=tspan()j和x=xmesh(k)时，之答案。 options为求解器的相关解法参数。详细说明参见odeset的使用方法。 注 1.MATLAB PDE求解器DdDe的算法，主要是将原米的椭圆型和地物线型偏微分 方程转化为一组常微分方程。此转换的过程是基于使用者所指定的msh点，以二阶空 间离散化(spatial discretization)技术为之(Keel and Berzins,.l990),然后以odel5s的指令 求解。采用ode15s的ode解法， 主要是因为在离散化的过程中，椭圆型偏微分方程被 转化为 一组代数方程，而抛物线型的偏微 “万程 被专化为 组跃立的分万程· ,变成一组同时伴有微分方程与代数方程的微分代数方程组， 的取 点(me 位置对解的精确度影响很大，君 iti 息时使田 e求器给出 动点取 占 外，若状态“在某些特定点上有较快速的变动时，亦 此处的点取密集些，以增加精确度。值得注意的是pdepe并不会自动做xmesh的自动 点，使用者必须观察解的特性，自行作取点的操作。一般而言，所取的点数至少需大于 3以上。 3.tspan的选取主要是基于使用者对那些特定时间的状态有兴趣而选定。而间距 248-
-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.若要获得特定位置及时间下的解，可配合以pdeval命令。使用格式如下： uout,duoutdx=pdeval(m,xmesh,ui,xout) 其中 m代表问题的对称性。m=0表示平板：m=1表示圆柱体：m=2表示球体。其意 义同pdepe中的自变量m。 xmesh为使用者在pdepe中所指定的输出点位置向量。xmesh=x,x,.,xw]。 即solU,:,)。也就是说其为pdepe输出中第i个输出在各点位置xmesh处， 时间固定为1，=1pan(U)下的解。 0为所欲内插输出点位置向量。此为使用者重新指定的位置向量 1o为基于所指定位置xO,固定时间1，下的相对应输出 duoutdx为相对应的du/d输出值。 ref KeelrD.and M Berzins"A method for the Spatial Discritization of Parabolic Equations in One Space Variable",SIAMJ.Sci.and Sat.Comput.,Vol.11.pp.1-32,1990. 以下将以数个例子，详细说明pdepe的用法。 32 求 界一维偏微分方程 2试解以下之偏微分方程式 at ax? 其中0≤x≤1，且满足以下之条件限制式 ①起始值条件 IC:x,0)=sin(πx) 而边界条件 BC1:0,)=0 BC2:e+0L0=0 注：本间题的解析解为u(xt)=e-r sin) 解下面将叙述求解的步骤与过程。当完成以下各步骤后，可进一步将其汇总为 主程序cx20 状后 步骤1将欲求解的偏微分方程改写成如式的标准式。 此即 cx,l4,au/ax)=π2 f(x.t.u.ou/ox)=ou s(x.1.u,du/dx)=0 和m=0。 步骤2编写偏微分方程的系数向量函数。 function [c,f,s]-ex20 1pdefun (x,t.u,dudx) -249
-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)