http://www.ma-xy.com
目录
第一章 偏微分方程 1
1.1 偏微分方程建模 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.1.1 热传导方程 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.1.2 弦振动方程 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.1.3 膜振动方程 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.2 偏微分方程基本理论 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.2.1
三大基本方程
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.2.2 偏微分基本定义 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.2.3 解的存在唯一性 * . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.2.4 偏微分方程组 * . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.2.5 偏微分方程外问题 * . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.2.6 偏微分方程稳定性 * . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.3 偏微分方程数值方法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.3.1 有限差分法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.3.2 变分原理 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
1.3.3 Rizt 变分原理 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
1.3.4 Galerkin 虚功原理 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
1.3.5 Ritz-Galerkin 方法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
1.3.6 有限元方法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
1.3.7 谱分析方法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
1.4 符号和函数空间注记 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
1.4.1 符号注记 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
1.4.2 函数空间注记 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
1.5 MATLAB 解偏微分方程 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
1.5.1 有限元解法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
1.5.2 MATLAB 应用实例 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
1
http://www.ma-xy.com
第一章 偏微分方程
1.1 偏微分方程建模
1.1.1 热传导方程
我们仅研究均匀介质的固体热传导问题，如果是关于液体和气体的热传导，则需要重新推导
热传导方程的形式。下面，我们考虑这样一个问题：对一物体 (固) 进行加热，问不同时间 t 时，
物体各部分的温度是多少？
可以想象，在物体内的一点产生初始热量 (最初物体温度均匀)，热量向各个方向 (x, y, z) 扩
散应该是一致的。不妨将问题简化到一维情况，比如仅考虑 x 方向上的。考虑一个极细杆 (一维)
的热传导过程，如图 (1.1) 所示
图 1.1: 一维热传导示意图
设杆的横截面积为 A，杆的长度为 L，假设同介质且密度均匀，密度为 ρ。我们用 u(x, t) 表
示 t 时刻 x 点处的温度，用 e(x, t) 表示热能密度，即单位体积上的热能，则 [x, x + ∆x] 内具有
的热能为 e(x, t)A∆x。
考虑 [x, x + ∆x] 段内热能的瞬时变化率，由热量守恒，我们知道热能的瞬时变化是瞬时从
x 面内进入的能量减去瞬时从 x + ∆x 面内出去的热量，再加上瞬时内部 ([x, x + ∆x] 段自身) 产
生的热量，有
∂
∂t
e(x, t)A∆x ≈ ϕ(x, t)A − ϕ(x + ∆x, t)A + Q(x, t)A∆x (1.1)
其中：ϕ(x, t) 为热通量，表示单位时间 (瞬时) 通过 x 处的热量；Q(x, t) 为热源，表示单位时间
物体自身所产生的热量。
将上式 (1.1) 两边同除 ∆x，并令 ∆x → 0，有
∂e
∂t
= lim
∆x→0
ϕ(x, t) −ϕ(x + ∆x, t)
∆x
+ Q(x, t )
=
∂ϕ
∂x
+ Q
1
http://www.ma-xy.com
1.1 偏微分方程建模 第一章 偏微分方程
下面，我们来讨论 e(x, t) 和 ϕ(x, t)。e(x, t) 为热能密度，其和温度 u(x, t) 的关系为
e(x, t) = c(x)ρ(x)u(x, t)
其中：c(x) 为物体比热容，ρ(x) 为物体密度。
ϕ(x, t) 表示单位时间 (瞬时) 进入的热量，它和温度 u(x, t) 有关，两侧温差越大，热能流动
速度越快，热能由高温物体流向低温物体。由傅里叶热传导定律，我们有
ϕ(x, t) = −k
0
∂u
∂x
其中：k
0
为导热系数，与物体的材料有关；
∂u
∂x
可以表示单位温差，与 ϕ 成比例。
整理方程，得到
cρ
∂u
∂t
=
∂
∂x
(k
0
∂u
∂x
) + Q
= k
0
∂
2
u
∂x
2
+ Q (1.2)
我们将 Q 视为恒定热源，c, ρ, k
0
皆与物体有关，是物体的属性，在给定具体的物体后，不妨将
其视为常量 (均匀分布)。上述方程 (1.2) 描述了热能是如何展开 (扩散)，是一个通用的扩散过程，
如化学污染物，大气污染物等等的扩散，皆可描述为上述扩散过程 (1.2)，因此，方程 (1.2) 也被
成为扩散方程。
由于各个方向 (x, y, z) 的相似性，我们很容易将上述方程 (1.2) 扩展到三维空间
cρ
∂u
∂t
= k
0
∂
2
u
∂x
2
+
∂
2
u
∂y
2
+
∂
2
u
∂z
2
+ Q (1.3)
不妨记
△u = ∇
2
u =
∂
2
u
∂x
2
+
∂
2
u
∂y
2
+
∂
2
u
∂z
2
称 ∇
2
u 为 u 的双哈密顿算子形式，△u 为 u 的拉普拉斯算子形式。于是上式 (1.3) 写为算子形
式有
cρ
∂u
∂t
= k
0
∇
2
u + Q
我们定义热传导方程 (1.3) 的稳态：当温度 u(x, t) 不再随时间而改变，即
∂u
∂t
= 0
此时的状态为稳态。于是，我们有
∇
2
u = −
Q
k
0
(1.4)
称上述方程 (1.4) 为泊松方程。∇
2
u = 0 为拉普拉斯方程 (拉普拉斯算子 ∇
2
的方程)。
http://www.ma-xy.com 2 http://www.ma-xy.com
http://www.ma-xy.com
第一章 偏微分方程 1.1 偏微分方程建模
上面的热传导方程 (1.3) 比较适用于物体形状是长方体的情况，当物体是圆柱或者球体时，
可以进行相应的坐标变换，例如：对 x, y, z 进行柱面坐标变换，令
x = r cos θ
y = r sin θ
z = z
有
∇
2
u =
1
r
∂
∂r
r
∂u
∂r
+
1
r
2
∂
2
u
∂θ
2
+
∂
2
u
∂z
2
对 x, y, z 进行球面坐标变换，令
x = ρ sin ϕ cos θ
y = ρ sin ϕ sin θ
z = ρ cos ϕ
有
∇
2
u =
1
ρ
2
∂
∂ρ
ρ
2
∂u
∂ρ
+
1
ρ
2
sin ϕ
∂
∂ϕ
sin ϕ
∂u
∂ϕ
+
1
ρ
2
sin
2
ϕ
∂
2
u
∂θ
2
1.1.2 弦振动方程
考虑柔软均匀的细弦，将弦两端固定，当受到与平衡位置垂直的外作用力时，弦开始振动。
假设运动发生在同一平面 (仅上下振动)，求不同时间 t 弦上各点的坐标，也就是求不同时间 t 的
弦曲线。
¬建系
以弦左端点处为坐标原点 O 点，以弦平衡 (静止) 位置为 x 轴，垂直向上为 u 轴，建立坐标
系 xOu，如图 (1.2) 所示
图 1.2: 弦振动示意图
定量
设弦长为 l，外力为 f(x)，重力为 G，求 t 时刻各点位置，即 u(x, t)，(可以尝试忽略 G)。
®微元法进行受力分析
在微元段 [x, x + ∆x] 上进行受力分析，微元段长度为 ∆S ≈ ∆x(当 ∆x → 0 时)，弦受重力、
外力和张力的作用，其受力分析如图 (1.3) 所示
http://www.ma-xy.com 3 http://www.ma-xy.com
http://www.ma-xy.com
1.1 偏微分方程建模 第一章 偏微分方程
图 1.3: 弦受力分析示意图
由于弦很柔软，所以张力 T 是沿切线方向的，并且由胡克定律可知，T 与时间 t 无关，否则
要写成 T (x, t)。由短时间内受力平衡，我们有
T (x + ∆x
0
) cos α
2
= T (x) cos α
1
T (x + ∆x
0
) sin α
2
+ f
0
∆x = T (x) sin α
1
+ G
0
+ ρ∆x
∂
2
u
∂t
2
其中：ρ 为单位弦长的密度；
∂
2
u
∂t
2
表示 u 方向上的加速度，f
0
(x, t) 为力 F 单位长度的力。
在 ∆x → 0 时，有 G
0
→ 0，并且有 cos α
1
≈ cos α
2
≈ 1，以及
sin α
1
≈ tan α
1
=
∂u(x, t)
∂x
sin α
2
≈ tan α
2
=
∂u(x + ∆x, t)
∂x
我们可以假定 T (x) ≈ T (∆x + x) = T
0
，即张力 T 为常量，将其带入到 (1.5) 方程当中，有
ρ∆x
∂
2
u
∂t
2
≈ T
0
∂u(x + ∆x, t)
∂x
−
∂u(x, t)
∂x
+ f
0
∆x (1.5)
将上式两边同除 ∆x，并且令 ∆x → 0，有
ρ
∂
2
u
∂t
2
= T
0
∂
2
u
∂x
2
+ f
0
(1.6)
由弦均匀，我们知道密度 ρ 为常数，记
T
0
ρ
= a
2
,
f
0
ρ
= f(x, t)
带入式 (1.6)，整理后有
ρ
∂
2
u
∂t
2
= a
2
∂
2
u
∂x
2
+ f (x, t) (1.7)
上式 (1.7) 即为弦振动的一般方程。
1.1.3 膜振动方程
我们应该在二维情况下考虑薄膜的振动，三维情况下考虑“体”的振动。这里，我们仅对振
动的一维情况 (1.6) 做一个简单的推广：在二维情况下，有薄膜振动方程
ρ
∂
2
u
∂t
2
= T
0
∂
2
u
∂x
2
+
∂
2
u
∂y
2
+ f (x, y, t) (1.8)
http://www.ma-xy.com 4 http://www.ma-xy.com
http://www.ma-xy.com
第一章 偏微分方程 1.2 偏微分方程基本理论
在三维情况下 (声波、电磁波等在空间中的传播)，有如下振动方程
∂
2
u
∂t
2
= a
2
∂
2
u
∂x
2
+
∂
2
u
∂y
2
+
∂
2
u
∂z
2
+ f (x, y, z, t) (1.9)
1.2 偏微分方程基本理论
1.2.1 三大基本方程
在介绍经典偏微分的三大方程之间，我们先来介绍一下哈密顿算子和拉普拉斯算子：定义哈
密顿算子 ∇ 为
∇ = i
∂
∂x
+ j
∂
∂y
+ k
∂
∂z
= grid
其中：i, j, k 是 x, y, z 方向上的单位向量，哈密顿算子是具有方向的算子。
定义拉普拉斯算子 △ 为
△ =
∂
∂x
2
+
∂
∂y
2
+
∂
∂z
2
拉普拉斯算子是不具有方向的标量算子，并且有
△ = ∇∇ = ∇
2
下面，来介绍经典偏微分的三大基本方程：波动方程 (双曲型方程)、扩散方程 (抛物型方程)
和稳态方程 (椭圆型方程)。
波动方程 前面引入的弦振动和膜振动等方程即为波动方程，波动方程的特点是方程中含有未知
量 u(待求函数) 关于时间 t 的二阶偏导数，一般形式为
∂
2
u
∂t
2
= a
2
∂
2
u
∂x
2
+
∂
2
u
∂y
2
+
∂
2
u
∂z
2
+ f (x, y, z, t)
此方程用于描述弦、膜、波、传输线上电流电压等物理量的变化规律。可以看出，由于 u 是时变
的，所以，求解得到的 u 是一个随时间 t 变化而变化的三维立体 (曲面)u(x, y, z, t)，即是一个动
图。
扩散方程 前面引入的热传导方程即为扩散方程，扩散方程的特点是方程中含有 u 对时间 t 的
一阶偏导，其一般形式为
∂u
∂t
= a
2
∂
2
u
∂x
2
+
∂
2
u
∂y
2
+
∂
2
u
∂z
2
+ f (x, y, z, t)
此方程用于描述热量扩散、液体渗透、气体和半导体中杂质的扩散过程。
http://www.ma-xy.com 5 http://www.ma-xy.com
http://www.ma-xy.com
1.2 偏微分方程基本理论 第一章 偏微分方程
稳态方程 未知函数 u 不随时间 t 变化而变化，即对 ∀n = 1, 2, . . . ,，有
∂u
∂t
n
= 0
稳态方程的典型代表是泊松方程和拉普拉斯方程
0 =
∂
2
u
∂x
2
+
∂
2
u
∂y
2
+
∂
2
u
∂z
2
+ f (x, y, z)
0 =
∂
2
u
∂x
2
+
∂
2
u
∂y
2
+
∂
2
u
∂z
2
可以看出，这里的 u 不是时变的，与时间 t 无关，所以解为一固定的三维曲面。
1.2.2 偏微分基本定义
1. 偏微分方程的一般形式为
F
x,
∂u
∂x
1
, . . . ,
∂u
∂x
n
, . . . ,
∂
m
u
∂x
m
1
1
∂x
m
2
2
, . . . , ∂x
m
n
n
= 0 (1.10)
其中：x = (x
1
, x
2
, . . . , x
n
) 为自变量，u 为未知函数，F 为已知的函数形式，
m
i
= m。
2. 偏微分方程的阶数。方程中最高阶偏导数的阶即为方程的阶，对上述方程 (1.10)，其阶数
为 m。
3. 显隐性。偏微分方程的显隐性由方程的主函数 F 决定，如果 F 为显函数，则偏微分方程
也为显式方程。
4. 方程的线性与非线性：
a
∂
2
u
∂x
2
+ b
∂u
∂x
+ cu = f(x, t)
为线性偏微分方程
a
sin
∂
2
u
∂x
2
+
b
cos
∂u
∂x
+
cu
= 0
为非线性偏微分方程。
5. 常系数方程和变系数方程。如果方程的系数是常量，则称为常系数微分方程
a
∂
2
u
∂x
2
+ b
∂u
∂x
+ cu = f(x, t)
为常系数方程
a(x, y, t)
∂
2
u
∂x
2
+ b(x, t )
∂u
∂x
+ c(t)u = f(x, t)
为变系数方程。
6. 齐次与非齐次方程。
a
∂
2
u
∂x
2
+ b
∂u
∂x
+ cu = f(x, t)
为非齐次方程
a
∂
2
u
∂x
2
+ b
∂u
∂x
+ cu = 0
为齐次方程。
http://www.ma-xy.com 6 http://www.ma-xy.com
http://www.ma-xy.com
第一章 偏微分方程 1.2 偏微分方程基本理论
7. 边值条件
初值条件 像常微分方程有初值条件和边值条件那样，偏微分方程也有相应的条件。比如，我们
可以要求偏微分方程的初始状态为某一具体的函数 (曲面) 形式。不失为一般性，考虑二维平面
xOy，在区域 Ω ∈ R
2
上求一个函数 u(x, y, t)，这个函数 u 是时变的，要求这个函数在初始时刻
t
0
的状态为 φ(x, y)，即
u(x, y, t)
t=t
0
= u(x, y, t
0
) = φ(x, y)
上面介绍的是偏微分方程的初值条件，可以想象，求一个时变的曲面 u，要求其在初始时刻
t
0
的曲面形状就为 φ(x, y)。当然，这个条件在微分方程当中是最简单的，实际问题会有许许多
多的条件 (约束)。下面，我们再介绍两类常见的边值条件：第一边值条件 (Dirichlet 条件)、第二
边值条件 (Neumann 条件) 以及它们的混合。
第一边值条件 不失为一般性，考虑二维平面 xOy，在区域 Ω ∈ R
2
的边界 ∂Ω 上，有一个连续
函数 g (x, y)。我们求一个函数 u(x, y, t)，要求 u(x, y, t) 不仅满足某一偏微分方程，而且在边界
∂Ω 上等于 g(x, y)。可以想象，在三维空间中一直有一个曲面 u(x, y, t)，每个时刻 t 这个曲面是
不同的，但要求每个时刻 t 它在边界 ∂Ω 上等于 g(x, y)，即
u(x, y)
∂Ω
= g(x, y)
第二边值条件 在区域 Ω ∈ R
2
的边界 ∂Ω 上给定一个连续函数，求 u(x, y, t)，要求 u(x, y, t) 不
仅在 Ω 内满足偏微分方程，而且在 ∂Ω 上的任意一个沿 ∂Ω 的外法线方向 n 的方向导数
∂u
∂n
存
在，并且等于 g 在该点的函数值，即
∂u
∂n
∂Ω
= g
混合边值条件 我们将前面的第一边值条件和第二边值条件混合，可以得到混合边值问题
∂u
∂n
+ bu
∂Ω
= g
限于篇幅，关于偏微分方程解的存在唯一性、偏微分方程组、偏微分方程外问题和稳定性这
里不做介绍，有兴趣的可以参考相关的偏微分书籍。
若 PDE 在附加条件 ∂Ω 和求解域 Ω 的一定要求下，它的解在已知度量的某函数类中存在并
且唯一，并且关于附加条件是稳定的，那么就称该定解问题在相应的函数类中为适定的。在偏微
分方程的解法上，我们主要介绍 PDE 的数值解法，而其古典解析方法：傅里叶变换法、拉普拉
斯变换法、分离变量法和格林函数法等在一般理论书籍上均由介绍。后面，我们将主要讨论数值
方法：有限差分法、变分法、有限元方法 (基于变分法) 和谱方法 (基于变分法)。
http://www.ma-xy.com 7 http://www.ma-xy.com
http://www.ma-xy.com
1.3 偏微分方程数值方法 第一章 偏微分方程
1.2.3 解的存在唯一性 *
1.2.4 偏微分方程组 *
1.2.5 偏微分方程外问题 *
1.2.6 偏微分方程稳定性 *
1.3 偏微分方程数值方法
1.3.1 有限差分法
二维泊松方程第一边值条件
由于泊松方程的解是静态的 (与时间 t 无关)，因此，我们用泊松方程来介绍差分法的思想。
差分法的思想与常微分方程中的欧拉法是相似的，考虑如下二维情况下的含第一边值条件的泊松
方程
0 =
∂
2
u
∂x
2
+
∂
2
u
∂y
2
+ f (x, y), (x, y) ∈ Ω
u(x, y)
∂Ω
= g(x, y), (x, y) ∈ ∂Ω
将其写为一般形式
− △u = −
∂
2
u
∂x
2
+
∂
2
u
∂y
2
= f(x, y)
u(x, y)
∂Ω
= g(x, y), (x, y) ∈ ∂Ω
其中：Ω 为自变量 x, y 的定义域，f (x, y) 和 g(x, y) 为已知的函数形式，u(x, y) 为待求函数，边
界 ∂Ω 为分段光滑的曲线。注：我们有必要讨论函数 f, g, u 的类型 (函数类)。
我们记 Γ = ∂Ω，记 Lu = −△u。这样，我们就可以避免写 − 符号了。有限差分法的基本思
想是通过对求解区域网格化离散化，然后建立差分方程求解函数 u 的近似值。我们在区域 Ω 上
取两组平行于 x, y 轴的直线，不妨设为等距
x : x
i
= ih i = 0, ±1, ±2, . . .
y : y
j
= jk j = 0, ±1, ±2, . . .
其中：h, k 为步长，有 x
i+1
− x
i
= h 以及 y
j+1
− y
j
= k。
图
1.4:
差分示意图
http://www.ma-xy.com 8 http://www.ma-xy.com
http://www.ma-xy.com
第一章 偏微分方程 1.3 偏微分方程数值方法
如图 (1.4)，我们设四点所构成的区域为 ∆Ω，Ω = {(x, y)|f }。用 Ω
h
表示所有属于 Ω 网格
内的网格点/节点，用 ∂Ω
h
表示所有位于 ∂Ω 上的网格点或者节点。
为了方便，不妨考虑 Ω = (x, y)|a < x < b, c < y < d 的矩形区域，如图 (1.5) 所示
图 1.5: 矩形差分示意图
将 [a, b] 进行 N 等分，将 [c, d] 进行 M 等分，有 h =
b−a
N
, k =
d−c
M
，并且，
x : x
i
= a + ih i = 0, 1, . . . , N
y : y
j
= c + jk j = 0, 1, . . . , M
两簇指点的交点 (x
i
, y
j
) 就是网格节点，记为
Ω
h
= {(x
i
, y
j
)|0 < i < N, 0 < j < M}
∂
Ω
h
=
{
(
x
i
, y
j
)
|
i
= 0
, N, j
= 0
,
1
, . . . , M
;
j
= 0
, M, i
= 0
,
1
, . . . , N
}
边值条件为：在 (x, y) ∈ ∂Ω 处有 u = g(x, y)。像一般的数值解法那样，我们自然希望得到待
求函数 u 在各个网格点 (x
i
, y
j
) 处的函数值 u(x
i
, y
j
)，这样就可以用它来估计待求函数 u。现在
的问题是，如果我们知道了 u 在 (x
1
, y
1
) 处的值 u
1
，如何求 u
2
？一般的，如果我们知道了 (x
i
, y
j
)
处的函数值 u
i,j
，如何求 u
i+1,j
处的函数值？
像常微分方程那样，我们用差分方程组来进行求解。我们仅就 x 方向来看，已知 u 在 (x
i
, y
j
)
的函数值 u
i
满足偏微分方程
−△u(x
i
, y
j
) = −
∂
2
u
∂x
2
(x
i
, y
j
) +
∂
2
u
∂y
2
(x
i
, y
j
)
= f(x
i
, y
j
) (1.11)
类似
f
′
(x
i
) =
f(x
i+1
) − f (x
i
)
h
−
h
2
y
′′
(ξ
i
)
我们将 u 在 (x
i
± h, y
j
) 处进行泰勒展开，有
u(x
i
± h, y
j
) = u(x
i
h, y
j
) ± h
∂u
∂x
(x
i
, y
j
) +
h
2
2
∂
2
u
∂x
2
(x
i
, y
j
)
±
h
3
6
∂
3
u
∂x
3
(x
i
, y
j
) +
h
4
24
∂
4
u
∂x
4
(ξ
±
, y
j
) (1.12)
其中：x
i
⩽ ξ
+
⩽ x
i+1
，x
i−1
⩽ ξ
−
⩽ x
i
。
http://www.ma-xy.com 9 http://www.ma-xy.com
http://www.ma-xy.com
1.3 偏微分方程数值方法 第一章 偏微分方程
根据上式 (1.12)，有
u(x
i
+ h, y
j
) − 2u(x
i
, y
j
) + u(x
i
− h, y
j
)
= h
2
∂
2
u
∂x
2
(x
i
, y
j
) +
h
4
24
∂
4
u
∂x
4
(ξ
+
, y
j
) +
∂
4
u
∂x
4
(ξ
−
, y
j
)
= h
2
∂
2
u
∂x
2
(x
i
, y
j
) +
h
4
12
∂
4
u
∂x
4
(ξ, y
j
) (1.13)
其中：x
i−1
⩽ ξ ⩽ x
i+1
。
在上式的最后一项中，我们利用了和的平均值定理，即如果有 g
i
⩾ 0 且有
J
i=1
g
i
= 1，那
么 ∃c，min x
i
⩽ c ⩽ max x
i
，使得
J
i=1
g
i
f(x
i
) = f (c) 成立。类似 x，我们对 y 也有
u(x
i
, y
j
+ k) −2u(x
i
, y
j
) + u(x
i
, y
j
− k)
= k
2
∂
2
u
∂y
2
(x
i
, y
j
) +
k
4
12
∂
4
u
∂y
4
(x
i
, η) (1.14)
其中：y
j−1
⩽ η ⩽ y
j+1
。
将式 (1.13) 除以 h
2
，式 (1.14) 除以 k
2
，然后将两式相加，可得
∆u(x
i
, y
j
) =
u(x
i
+ h, y
j
) − 2u(x
i
, y
j
) + u(x
i
− h, y
j
)
h
2
+
u(x
i
, y
j
+ k) −2u(x
i
, y
j
) + u(x
i
, y
j
− k)
k
2
−
h
4
12
∂
4
u
∂x
4
(ξ, y
j
) −
h
4
12
∂
4
u
∂y
4
(x
i
, η) (1.15)
上式忽略高阶项
∂
4
u
∂∗
4
，即可得到 ∆u(x
i
, y
j
) 的近似计算公式，然后将 ∆u(x
i
, y
j
) 带入 (1.11) 当中，
即可得到如下关系式
f(x
i
, y
j
) = −
u(x
i
+ h, y
j
) − 2u(x
i
, y
j
) + u(x
i
− h, y
j
)
h
2
−
u(x
i
, y
j
+ k) −2u(x
i
, y
j
) + u(x
i
, y
j
− k)
k
2
(1.16)
由此，我们得到了 u(x
i
, y
j
) 与 u(x
i+1
, y
j
)、u(x
i−1
, y
j
)、u(x
i
, y
j+1
) 和 u(x
i
, y
j−1
) 的等式关系
(1.16)，其示意图如图 (1.6) 所示，我们称此差分方程为泊松方程的五点差分方程，当然，选用不
同阶的泰勒展开公式会得到更多点的差分方程 (比如：九点差分方程)。
图 1.6: 五点位置示意图
综合所有的内点 (x
i
, y
j
) ∈ Ω
h
，我们会得到一个含有 (N − 1) × (M − 1) 个未知量 u
i,j
(i =
1, . . . , N − 1, j = 1, . . . , M − 1)(由于边界值已知，所以 i = 0, N ，y = 0, M 的值已知) 的
http://www.ma-xy.com 10 http://www.ma-xy.com
http://www.ma-xy.com
第一章 偏微分方程 1.3 偏微分方程数值方法
(N − 1) × (M − 1) 阶线性方程组。结合边界点方程
u(x
i
, y
j
) = g(x
i
, y
j
) (x
i
, y
j
) ∈ ∂Ω
h
由此得到 u(x, y) 的数值解。记由线性方程组得到的解为 ˆu(x
i
, y
j
)，它是 u(x
i
, y
j
) 的估计。为书
写简便，我们将 u(x
i
, y
j
) 记为 u
i,j
，同样有 u
i+1,j
, u
i−1,j
, u
i,j+1
, u
i,j−1
，以及它们对应的估计值
ˆu
i,j
。
问一：这种差分方法的数值计算误差哪里来？1. 步长 h, k 的大小；2. 被忽略的高阶项的大
小。问二：差分方程解的唯一性和收敛性如何？
下面，我们讨论 (1.16) 的截断误差。在计算 (1.16) 的时候，我们忽略了 (1.15) 的高阶项
∂
4
u
∂∗
4
从而带来误差，令
max
(x,y )∈Ω
∂
4
u
∂x
4
(x, y)
⩽ M
4
, max
(x,y )∈Ω
∂
4
u
∂y
4
(x, y)
⩽ M
4
则有如下截断误差项的估计
|R
i,j
(u)| =
h
4
12
∂
4
u
∂x
4
(ξ, y
j
) +
h
4
12
∂
4
u
∂y
4
(x
i
, η)
⩽
M
4
12
(h
2
+ k
2
)
所以，五点差分格式的数值计算误差的阶为 |R
i,j
(u)| = O(h
2
+ k
2
)。
再者，我们讨论差分方法的唯一性和收敛性。这里我们直接给出其收敛性和唯一性定理：
定理 (差分法的收敛性) 如果泊松方程第一边值问题的解 u(x, y) 在 Ω ∪∂Ω 上有四阶连续偏
导数 (即 u 是函数类中的一员)，则其五点差分的数值解一致收敛到精确解，解唯一并且有
max
Ω
h
|
ˆ
u
i,j
−
u
i,j
|
⩽
C
12
max(k, h)
2
其中：C 是与 h, k 无关的常数。
二维泊松方程第二边值条件
对于 Poisson 方程的 Dirichlet 边值问题，前一小节已经详细讨论了其差分格式的建立。可
以看出，只需要将边界点处的函数值直接带入差分格式中，就能获得 Poisson 方程 Dirichlet 边
值问题的有限差分格式。那么对于 Neumann 边值问题，差分格式的边界条件又如何处理呢？下
面具体讨论这一问题，不妨设 Poisson 方程的第二类边界条件 (即 Neumann 边界条件) 为
∂u
∂n
= α(x, y), (x, y) ∈ ∂Ω
为了简单起见，设 Ω = {(x, y) |0 < x < 1, 0 < y < 1}。考虑步长为 h 的正方形网格，显然，
在正方形的四个顶点上，外法线方向导数
∂u
∂n
没有定义，这表示 α(x, y) 在顶点处不连续。一般情
况下，取平均值来作为不连续点上的函数值。下面讨论 Neumann 边值问题的有限差分格式，即
对 (N + 1)
2
个网格点来构造离散的差分格式。
http://www.ma-xy.com 11 http://www.ma-xy.com
http://www.ma-xy.com
1.3 偏微分方程数值方法 第一章 偏微分方程
首先，对于所有的内节点 (x
i
, y
j
) ∈ Ω
h
，由式 (1.16)，可得离散的差分格式为
−
1
h
2
[u
i+1,j
+ u
i−1,j
+ u
i,j+1
+ u
i,j−1
− 4u
i,j
] = f
i,j
i = 1, 2, . . . , N − 1; j = 1, 2, . . . , N − 1 (1.17)
其次，对于边界点 (x
i
, y
j
) ∈ ∂Ω
h
，也需要对其进行离散，从而获得相应的差分方程。为了获
得差分方程，需要对网格区域 Ω
h
∪ ∂Ω
h
进行扩展，即对正方形网格区域四周各增加一排步长为
h 的节点，以左边界 x = 0 为例，在边界 x = 0 的左边引进一条直线 x = −h，将平行于 x 轴的
网格线 y = y
j
= jh 相应向左延长使其与直线 x = −h 相交，相应的一排交点即为新增加的网格
点，记作 (x
−1
, y
j
) = (−h, y
j
)，对其他三条边界作类似处理增加网格点。可以看出，新增加的节
点不属于网格区域 Ω
h
∪ ∂Ω
h
，当然也就不属于求解区域 Ω。不防将这些新增加的节点分别记作
(x
−1
, y
j
), (x
N+1
, y
j
), (j = −1, 0, . . . , N + 1)
(x
i
, y
−1
), (x
i
, y
N+1
), (i = 0, 1, . . . , N)
利用这些新增加的节点，根据第二边值条件，对于边界点，可以获得相应的离散格式。根据法向
导数与偏导数的关系，利用一阶中心差商代替导数，以 x = 0 为例，其上的网格点的方向导数
∂u
∂n
0,j
=
cos
π
∂u
∂x
0,j
+
sin
π
∂u
∂y
0,j
= −
∂u
∂x
0,j
≈
1
2h
(u
−1,j
− u
1,j
)
根据第二边值条件
∂u
∂n
0,j
= α
0,j
，即可得 x = 0 上导数边值条件的离散差分格式，而对于其他
三条边界，可作类似的离散化处理，综合可得在 x = 0, x = 1, y = 0, y = 1 上的导数边值条件的
离散差分格式分别为
1
2h
(u
−1,j
− u
1,j
) = α
0,j
j = 1, 2, . . . , N
1
2h
(u
N+1,j
− u
N−1,j
) = α
N,j
j = 1, 2, . . . , N
1
2h
(u
i,−1
− u
i,1
) = α
i,0
j = 1, 2, . . . , N
1
2h
(u
i,N+1
− u
i,N−1
) = α
i,N
j = 1, 2, . . . , N
(1.18)
此处，α
i,j
= α(x
i
, y
j
)。接下来，需要消去式 (1.18) 中新增加的节点处的函数值。我们以 x = 0
上的边值条件为例，消去式 (1.18) 中第一个离散格式中的函数值 u
−1,j
。这里需要借助前面构造
的五点差分格式来进行处理。不防在式 (1.17) 中取 i = 0，则有
u
−1,j
= −h
2
f
0,j
− u
1,j
− u
0,j+1
+ 4u
0,j
代入 (1.18) 第一式中，可得
4u
0,j
− 2u
1,j
− u
0,j−1
− u
0,j+1
= 2 hα
0,j
− h
2
f
0,j
j = 1, 2, . . . , N − 1
(1.19)
http://www.ma-xy.com 12 http://www.ma-xy.com
http://www.ma-xy.com
第一章 偏微分方程 1.3 偏微分方程数值方法
可以看出上式消去了新增加的网格点 u
−1,j
，故式 (1.19) 即为在 x = 0 上的差分格式。用同样的
方法，在 x = 1, y = 0, y = 1 上的网格点分别建立如下的离散格式：
4u
N,j
− 2u
N−1,j
− u
N,j−1
− u
N,j+1
= 2 hα
N,j
− h
2
f
N,j
j = 1, 2, . . . , N − 1
4u
i,0
− 2u
i,1
− u
i−1,0
− u
i+1,0
= 2 hα
i,0
− h
2
f
i,0
i = 1, 2, . . . , N − 1
4u
i,N
− 2u
i,N−1
− u
i−1,N
− u
i+1,N
= 2 hα
i,N
− h
2
f
i,N
i = 1, 2, . . . , N − 1
(1.20)
最后我们来讨论四个顶点 (x
0
, y
0
), (x
0
, y
N
), (x
N
, y
0
), (x
N
, y
N
) 上的离散格式。以顶点 (x
0
, y
0
) 为
例，在式 (1.19) 中令 j = 0，有
4u
0,0
− 2u
1,0
− u
0,−1
− u
0,1
= −h
2
f
0,0
+ 2hα
0,0
又
1
2h
(u
0,−1
− u
0,1
) = −α
0,0
消去 u
0,−1
可得
4u
0,0
− 2u
1,0
− 2u
0,1
= 4 hα
0,0
− h
2
f
0,0
(1.21)
同理可得其他三个顶点的离散格式为
4u
0,N
− 2u
1,N
− 2u
0,N−1
= 4 hα
0,N
− h
2
f
0,N
4u
N,0
− 2u
N,1
− 2u
N−1,0
= 4 hα
N,0
− h
2
f
N,0
4u
N,N
− 2u
N−1,N
− 2u
N,N−1
= 4 hα
N,N
− h
2
f
N,N
(1.22)
综合上面的分析，联立内节差分方程 (1.17)，边界点和顶点上的离散差分格式 (1.19)-(1.22)，
即得 Poisson 方程的 Neumann 边值问题的有限差分格式，很容易看出，在正方形区域 Ω =
{(x, y)|0 < x < 1, 0 < y < 1} 上的 (N + 1)
2
个节点上的差分方程是一个 (N + 1)
2
阶的线性方程
组。不难验证，这个线性方程组的系数矩阵是奇异的，且解不唯一。可以注意到，在 Poisson 方
程中，若取 f (x, y) = 0，Poisson 方程就变为 Laplace 方程，可以非常容易地写出 Laplace 方程的
Neumann 边值问题的有限差分方程。对于 Poisson 方程的混合边值条件，采用与分析 Neumann
边值问题相同的方法，可类似于得到 Poisson 方程的混合边值问题的离散的有限差分格式。值得
注意的是，混合边值问题的差分方程的解是唯一的。
1.3.2 变分原理
Sobolev 空间
继续考虑前面的弦振动问题 (1.1.2)，在 f(x) 的作用下，弦振动满足如下方程
∂
2
u
∂t
2
= a
2
∂
2
u
∂x
2
+ f (x, t)
http://www.ma-xy.com 13 http://www.ma-xy.com
http://www.ma-xy.com
1.3 偏微分方程数值方法 第一章 偏微分方程
我们假设 t 时刻弦处于平衡位置
∂
2
u
∂t
2
= 0 ，则平衡位置 ¯u(x) 满足
a
2
∂
2
u
∂x
2
= f(x)
u(0) = 0
u(l) = 0
(1.23)
其中：l 为弦长，u(l) = 0 表示端点为 0。这样，求弦的平衡位置 ¯u(x) 就变为常微分方程 (ODE)
中的两点边值问题。
另一方面，由“极小位能原理”，弦平衡位置 ¯u(x) 是一切位置 {u(x)} 中，使位能最小者。其
中，位能的定义为
J(u) = W
内
+ W
外
=
1
2
l
0
a
2
(u
′
x
)
2
dx −
l
0
f · udx
=
1
2
l
0
(a
2
u
′2
− 2uf )dx (1.24)
故，求平衡位置 ¯u(x) 即可解上式 (1.24) 的变分问题 J(u)。
可以想象，微分方程问题与变分问题具有某种联系，因为它们同样是在一个函数空间中寻找
最优函数 u，以达到某种目的。我们在谈论函数时，经常讨论集合、元素、映射 x ∈ R
n
f
−→ y ∈ R。
现在，我们有一个集合，这个集合是一个函数集合 {u(x)}。下面，就来看一下函数集合。
设 I = (a, b)，
¯
I = [a, b]，L
2
(I) 是定义在 I 上的平方可积的可测函数组成的集合 (空间)，定
义空间的内积和范数为
(f, g) =
b
a
f · gdx, f, g ∈ L
2
(I)
||f|| =
(f · f) =
b
a
|f|
2
dx
2
f ∈ L
2
(I)
L
2
(I) 关于加法和数乘是线性空间，关于 ( , ) 是完备内积空间，因此，L
2
(I) 是 Hilbert 空间。
所谓完备是指，Cauchy 收敛定理在 L
2
(I) 中成立，即 L
2
(I) 中任意一函数列 {f
n
}，如果函数列
关于度量 || ·|| 满足
||f
n
− f
m
|| → 0, (n, m → ∞)
则必有 f ∈ L
2
(I)，使得
||f
n
− f || → 0, (n, m → 0)
并记为 lim
n→∞
f
n
= f。
为了使 J 有意义，我们有必要要求 u ∈ L
2
(I), u
′
∈ L
2
(I)，即平方可积可测
b
a
|u|
2
dx < ∞
b
a
|u
′
|
2
dx < ∞
现在，我们可以定义一个空间 H(I)
H(I) = {f|f ∈ L
2
(I), f
′
∈ L
2
(I)}
http://www.ma-xy.com 14 http://www.ma-xy.com
http://www.ma-xy.com
第一章 偏微分方程 1.3 偏微分方程数值方法
但仅仅如此要求，是不够的，可能会遇到许多函数并不是处处可微的，所以有必要将导数的定义
f
′
扩大。我们引入下面的广义导数的概念:
定义 (广义导数) 我们用 C
∞
0
(I) 表示在 I 上无数次可微并且在端点 a, b 的某邻域内为 0 的
函数类 (函数集合)。设 f ∈ L
2
(I) 且为一次连续可微函数，若存在 g ∈ L
2
(I)，使得
b
a
g(x)φ(x)dx = −
b
a
f(x)φ
′
(x)dx
对 ∀φ ∈ C
∞
0
(I) 成立，则称 g(x) 为 f 在 I 的广义导数
¬
。
定义 (Sobolev 空间) 现在，我们可以来定义 Soboblev 空间了，定义
H
1
(I) = {f|f ∈ L
2
(I), f
′
∈ L
2
(I)}
为 (1 阶)Sobolev 空间，其中：f
′
为函数 f 的广义导数。
通俗地说，函数集合 H
1
(I) 中的函数具有这样的特点：函数 f 平方可积可测，并且它的广
义导数 f
′
平方可积可测。同理，我们可以根据 H
1
(I) 定义函数 f 的 m 阶广义导数，以及基于
m 阶广义导数的 m 阶 Sobolev 空间 H
m
(I)。
上面定义的 H
1
(I) 是在一维情况进行的 (I = [a, b])，下面，我们将域 I 推广到二维平面域
以及高维域中。设 Ω 为有界平面域，Γ 为分段光滑的简单闭曲线，
¯
Ω = Ω ∪ Γ 是 Ω 的闭包。用
C
∞
0
(Ω) 表示 Ω 上无穷次可微并且有紧支集的函数类，L
2
(Ω) 为 Ω 上平方可积可测的函数类，则
可以定义 Ω 上的广义导数和 Sobolev 空间 H
m
(Ω)。
下面，我们来讨论如下泛函 J(u) 的极值
J(u) =
b
a
f(x, u, u
′
)dx
其中：x ∈ R，u(x) : R → R。上述泛函取极值的必要条件为：u 满足如下欧拉方程
f
u
− f
u
′
x
− u
′
f
u
′
u
− u
′′
f
u
′
u
′
= 0
上述极值情况可以推广到多元函数的情况，如
J(u) =
Ω
∂u
∂x
2
+
∂u
∂y
2
+
∂u
∂z
2
+ 2udΩ
其中：u 为三元函数 u(x, y, z)。
1.3.3 Rizt 变分原理
考虑如下常微分方程 (ODE) 问题
Lu := −
d
2
u
dx
2
= f
a < x < b
u(a) = 0, u
′
(b) = 0
(1.25)
¬
注：δ 函数不属于 L
2
(I)；阶梯函数无广义导数。
http://www.ma-xy.com 15 http://www.ma-xy.com
http://www.ma-xy.com
1.3 偏微分方程数值方法 第一章 偏微分方程
其中：f ∈ C(a, b)。上述 ODE 问题与弦平衡方程 (1.23) 相同，我们对其作如下二次泛函
J(u) =
1
2
(Lu, u) −(f, u)
= −
1
2
b
a
u
′′
udx −
b
a
fudx (1.26)
其中：(f, u) 表示函数内积。
对上式右边的积分项 −
b
a
u
′′
udx 进行分部积分，并带入相应的边值条件，有
−
b
a
u
′′
udx = −
du
dx
u
b
a
+
b
a
d
2
u
dx
2
dx =
b
a
(u
′′
)
2
dx
注意：上式的中间部分采用了边值为 0 的条件。令
−a(u, v) =
b
a
u
′
v
′
dx
则 (1.26) 可以写为
J(u) =
1
2
a(u, u) −(f, u) (1.27)
于是，求两点边值问题就变为如下的变分问题
u
∗
= arg min
u∈H
1
E
J(u)
其中：H
1
E
为一维 Sobolev 空间 H
1
中满足 u(a) = 0 的函数组成的子空间，可写为交集的的形
式。
上述过程即为极小位能原理。我们记 C
k
(Ω) 为在域 Ω 上具有 k 阶连续偏导数的函数集，假
设 u ∈ C
2
(a, b) ∩C
1
(a, b) 满足方程 (1.27)，则称 u 为方程 (1.25) 的古典解。
注：我们对上面的 a(u, v) 是有要求的：
1、要求 a(u, v) 是对称形式的，即 ∀u, v ∈ H
1
(I)，有 a(u, v) = a(v, u)；
2、要求 a(u, v) 是正定的，即 ∀u ∈ H
1
E
，有 a(u, v) ⩾ r ||u||
2
1
；
3、a(u, v) 是连续的，即 ∀u, v ∈ H
1
，有 |a(u, v)| ⩽ M||u||
1
||v||
1
。
1.3.4 Galerkin 虚功原理
仍然考虑如下常微分方程 (ODE) 问题 (1.25)，用 v 乘以问题 (1.25) 方程两端，并在区间
[a, b] 进行积分，有
b
a
(Lu − f )vdx =
b
a
−
d
2
u
dx
2
vdx = 0
仍然利用分部积分法处理上式中间部分，并利用相应的边值条件，有
−
b
a
u
′′
vdx =
b
a
u
′
v
′
dx
http://www.ma-xy.com 16 http://www.ma-xy.com
http://www.ma-xy.com
第一章 偏微分方程 1.3 偏微分方程数值方法
令
−a(u, v) =
b
a
u
′
v
′
dx
则上式可以写为
a(u, v) − (f, v) = 0 (1.28)
上式 (1.28) 是常微分两点边值问题 (1.25) 的 Galerkin 变分问题。
对于 Galerkin 变分方程 (1.28) 的解和原两点边值问题 (1.25) 的解之间的关系，我们有
u ∈ C
2
是边值问题的解 ⇐⇒ u ∈ H
1
E
且是如下 Galerkin 变分方程的解
a(u, v) − (f, v) = 0, ∀v ∈ H
1
E
其中：H
1
E
是一个 I = [I
a
, I
b
] 域上的函数集合，是一维 Sobolev 空间 H
1
中满足 v(I
a
) = 0 的函
数组成的子空间，H
1
E
= {v|v ∈ L
2
(I), a(u, v) < ∞, v(I
a
) = 0}。
1.3.5 Ritz-Galerkin 方法
对常微的两点边值问题 (1.25)，Ritz 泛函形式为
J(u) =
1
2
a(u, v) − (f, u)
于是两点边值问题 Lu = f 等价于如下变分问题
u
∗
= arg min
u
J(u)
此为边值问题的极小位能原理。两点边值问题的另一个变分问题是 Galerkin 变分：求 u ∈ U，使
a(u, v) = (f, v) ∀v ∈ U
此为边值问题的虚功原理。极小位能原理对内积 a(u, v) 有要求 (1.3.3)，而虚功原理不要求 a(u, v)
对称和正定。
上述两种变分法的难点在于：如何在无穷维 (无数个可行个体) 空间 U 上求泛函的极小值
(函数 u
∗
)？Ritz-Galerkin 方法用有穷维空间来做近似替代，即不用找遍空间 U，而是找一个有
限的空间，有限空间中有有限可列个个体，因此问题变得可解。接下来的问题是如何选取有穷维
空间？
设 U
n
是 U 的 n 维子空间，φ
1
, φ
2
, . . . , φ
n
是 U
n
的一组基底，我们称其为基函数，则 U
n
中的任意一个函数 (个体)u
n
可表示为
u
n
=
n
i=1
c
i
φ
i
http://www.ma-xy.com 17 http://www.ma-xy.com
http://www.ma-xy.com
1.3 偏微分方程数值方法 第一章 偏微分方程
Ritz 方法 Ritz 法的目标是：求解系数 {c
i
}
n
i=1
，使下面的 J(u
n
) 取极小值。
J(u
n
) =
1
2
a(u
n
, v
n
) − (f, u
n
)
=
1
2
n
i,j=1
a(φ
i
, φ
j
)c
i
c
j
−
n
j=1
c
j
(f, φ
j
)
其中：J(u
n
) 是 c
1
, c
2
, . . . , c
n
的二次函数，并且由 a(u, v) 的对称性有 a(φ
i
, φ
j
) = a(φ
j
, φ
i
)。令
∂J(u
n
)
∂c
j
= 0
则 {c
i
} 满足下面的线性方程组。
n
i=1
a(φ
i
, φ
j
)c
i
= ( f, φ
j
) j = 1, 2, . . . , n
在求出 c
i
之后，我们就可以得出 u 的近似值 u
n
。上面是在有限维 n 维空间 U
n
上进行讨论
的，令 n → ∞，如果 lim
n→∞
n
i=1
c
i
φ
i
存在，我们就可得到函数
u
n
=
∞
i=1
c
i
φ
i
(x)
Ritz 方法得到一系列函数 u
1
(x), u
2
(x), . . . , u
n
(x)，我们设它们对应的变分值为 J
1
, J
2
, . . . , J
k
, . . . 。
在实际操作中，只要连续 m 次得到的值 J
i
, j
k
接近，那么得到的 u
k
(x) 是 u(x) 的较好的近似表
达。下面的问题是如何找基函数 {φ
i
}
n
i=1
？其实我们可以选取如下三角基和多项式基
{1, x, x
2
, . . . , x
n
, . . . }
{1, cos x, sin x, cos 2x, sin 2x, . . . }
Galerkin 方法 上面介绍的在有限维空间 U
n
求解函数 u 的近似是基于 Ritz 变分法的，下面，
我们介绍基于 Galerkin 变分法的有限维方法。
Step1. 建基。
选取一个相对完备的函数系 φ
1
, φ
2
, . . . , φ
k
，并且这些函数都可以满足常微的两点边值条件。
如果常微的边值条件为 u(a) = 0，则 φ
1
(a) = 0, φ
2
(a) = 0, . . . , φ
n
(a) = 0。
Step2. 构建近似解。
选取线性组合 u
n
(x) =
n
i=1
c
i
φ
i
为微分方程解 u 的近似解，显然近似解满足边值条件。其中，
系数 c
i
待定。
Step3. 构建 G 方程组。
用 φ
1
, φ
2
, . . . , φ
n
分别乘以 Lu
n
= f 的两边，然后再积分 (这就是 Galerkin 与 Ritz 的主要
区别)，有
Ω
(Lu
n
− f )φ
i
dx = 0, i = 1, 2, . . . , n
http://www.ma-xy.com 18 http://www.ma-xy.com
http://www.ma-xy.com
第一章 偏微分方程 1.3 偏微分方程数值方法
即
Ω
L(
n
i=1
c
i
φ
i
) − f )
φ
i
dx = 0, i = 1, 2, . . . , n (1.29)
Step4. 求解系数。
由积分后得到的关于 c
i
的 n 个方程组 (1.29) 可以解出系数 c
i
，即可得到 u 的估计 u
n
。
Step5. 终止条件。
如果连续多次得到的泛函值 J
i
的变化不大，则终止；否则，置 k = k + 1，返回 Step1。
注：U
n
为 U 上的 n 维有限子空间，φ
1
, φ
2
, . . . , φ
n
为基函数。U
n
可以由基 {φ
i
} 张成，我们记为
U
n
=
v
v =
n
i=1
c
i
φ
i
= span{φ
i
}
注意：无论是 Ritz 还是 Galerkin，我们都是通过逐渐扩大基 {φ
i
}
n
i=1
中基函数的个数 n 来实现
算法的。
上述 R-G 方法是在常微两点边值下引导出来的，对于偏微分方程，情况是相似的，具体情
况可以参考《微分方程数值解法 (第四版)》李荣华 P198。下面，我们仅作简单的介绍。考虑如
下二阶椭圆方程
− △u = −
∂
2
u
∂x
2
+
∂
2
u
∂y
2
= f(x, y)
u(x, y)
∂Ω
= g(x, y), (x, y) ∈ ∂Ω
(1.30)
其中：u 为待求函数；Ω ∈ R
2
为自变量 (x, y) 区域；Γ = ∂Ω ∈ R
2
为自变量取值边界。
仿照上面 R-G 方法，找方程对应的泛函
J(u) =
1
2
(−△u, u) −(f, u)
=
1
2
Ω
(−△u)u dxdy −
Ω
fu dxdy
其中：v 为 Ω 上的函数。由 Green 公式，我们有
Ω
(−△u)vdxdy =
Ω
(
∂u
∂x
∂v
∂x
+
∂u
∂y
∂v
∂y
)dxdy −
Γ
∂u
∂nvds
其中：n 表示边界 Γ 的单位外法向量。
若 u, v 满足边界条件 u|
Γ
= 0 ，则
Ω
(−△u)vdxdy =
Ω
(
∂u
∂x
∂v
∂x
+
∂u
∂y
∂v
∂y
)dxdy
令
a(u, v) =
Ω
(u
′
x
v
′
x
+ u
′
y
v
′
y
)dxdy
则可以将原问题 (1.30) 变为下列泛函极值问题
min
u∈U
J(u) =
1
2
a(u, v) − (f, u)
http://www.ma-xy.com 19 http://www.ma-xy.com
http://www.ma-xy.com
1.3 偏微分方程数值方法 第一章 偏微分方程
1.3.6 有限元方法
有限元方法的导出
有限元方法是 Ritz-Galerkin 方法的延伸，它用样条函数提供了一种选取“局部基函数”或
“分片多项式空间 (分段多项式函数的空间推广)”的技术，克服了 R-G 方法基函数选择的困难。
在古典 R-G 方法的基函数选择中，一般选取幂函数和三角函数等初等函数，而且往往要求基在
区域边界上满足边界条件。如果区域 Ω 是二维和三维的不规则区域，这样的基函数往往很难构
建出来。
Galerkin(虚功原理) 比 Ritz 应用更为广泛，因为它不要求 a(u, v) 对称和正定。下面，我们
就用 Galerkin 方法导出有限元方法。仍然讨论如下两点边值问题
Lu := −
d
2
u
dx
2
= f
a < x < b
u(a) = 0, u
′
(b) = 0
(1.31)
其中：f ∈ C(a, b)。上述方程对应的 G 泛函问题为：求 u，使得下式成立
a(u, v) = (f, v)
a(u, v) − (f, v) = 0, ∀v ∈ H
1
E
其中：H
1
E
是一个 I = [a, b] 域上的函数集合，是一维 Sobolev 空间 H
1
中满足 v(a) = 0 的函数
组成的子空间，H
1
E
= {v|v ∈ L
2
(I), a(u, v) < ∞, v(a) = 0}，注意：这里有两个 a，要注意区分
它们，一个是数值，一个是算子。
下面，我们就来寻找有限维函数空间 U
n
。对区域 Ω := I = [a, b] 进行离散化
a = x
0
< x
1
< ··· < x
n
= b
令 I
i
= [x
i−1
, x
i
] 为第 i 区间段，令 h
i
= x
i
− x
i−1
为第 i 区间段的长度。我们在每个小区间 I
i
上考虑线性函数，则在整个区间 I 上，函数为分段一次多项式的形式，记为 u
h
(x)，我们将所有
的分段一次多项式函数放在一起，形成函数空间 U
n
，显然 U
n
是有限维空间。
u
h
(x) =
x
i
− x
h
i
u
i−1
+
x − x
i−1
h
i
u
i
其中：{u
i
} 为待求量，x
i
∈ I
i
, i = 1, 2, . . . n。
不仅可以采用线性函数的形式，我们还可以在每个小区间 I
i
上考虑二次函数、三次函数等
其它的函数形式。由此，我们可以构建出有限维的函数空间 U
n
，但存留的问题是：该空间的基
函数是什么？如何求解参数 {u
i
}？
同一空间 U
n
可以有不用的基。我们对 U
n
考虑如图 (1.7) 不相关的基函数，左边的基函数
处理边值为 0 的情况，右边的基函数处理边值不为 0 的情况。
http://www.ma-xy.com 20 http://www.ma-xy.com
http://www.ma-xy.com
第一章 偏微分方程 1.3 偏微分方程数值方法
图 1.7: 分段线性的基函数示意图
下面，考虑基的一个具体形式 (边值不为 0)
φ
0
(x) =
1 +
x
1
− x
h
1
, x
0
⩽ x ⩽ x
1
0, 其它
φ
i
(x) =
1 +
x − x
i
h
i
, x
i−1
⩽ x ⩽ x
i
1 +
x − x
i
h
i+1
, x
i
⩽ x ⩽ x
i+1
0, 其它
φ
0
(x) =
1 +
x − x
n
h
n
, x
n−1
⩽ x ⩽ x
n
0, 其它
基于上面这种基 {φ
i
}，有限维空间 U
n
中的任意一个分段一次多项式 u
h
, (u
h
∈ U
n
) 都可以
用基函数表示
u
h
(x) =
n
i=1
φ
i
u
i
其中：我们可以看到函数 u
h
(x) 在点 x
i
的函数值为 u
i
= u
h
(x
i
)。
至此，我们已经构建了一个有限维空间 U
n
，以及它的基组 {φ
i
}，我们可以沿用 Golerkin 方
法得到最终的参数值 {u
i
}。
设 u 是 G 变分问题的解，v 为 U 中的任意函数，我们将它们用基函数估计 u, v，有
u
h
(x) =
n
j=1
u
j
φ
j
(x)
v
h
(x) =
n
j=1
v
i
φ
i
(x)
http://www.ma-xy.com 21 http://www.ma-xy.com
http://www.ma-xy.com
1.3 偏微分方程数值方法 第一章 偏微分方程
其中：u
h
, v
h
是 u, v 有限维空间中的估计。将上式带入到 G 变分问题当中，我们就可以形成如
下方程组，解之，即可得到 {u
i
}。
n
i=1
a(φ
i
, φ
i
)u
i
− (f, φ
j
) = 0, j = 1, 2, . . . , n
=⇒
n
i=1
a(φ
i
, φ
j
)u
i
=
b
a
fφ
j
dx, j = 1, 2, . . . , n
在实际处理过程中，我们并不把方程中系数矩阵的每个元素 a(φ
i
, φ
j
) 按 i 和 j 的顺序依次
计算，而是逐元计算其积分值，再累加起来的到系数矩阵。
有限元计算步骤
上面在常微边值问题的基础上介绍了有限元方法，下面用二维泊松方程来介绍有限元的计算
步骤，为了便于理解，我们对泊松方程进行简化，并采用三角形元来进行说明。
Lu := −△u = f
u
∂Ω
= g
(1.32)
为了简单，不妨令 f = x，g = 0，Ω 为 (0, −1), (0, 1), (1, 0) 构成的闭三角区域。
1. 区域 Ω 离散化 我们将 Ω 离散化为三角形元集合，如图 (1.8)，并设 N
e
为三角形元的个数，
e
k
, (k = 1, 2, . . . , N
e
) 为第 k 个三角形元 (单元 k)。e
k
有三个端点，分别记为 P
i
, P
j
, P
k
，i, j, k
为逆时针方向，则 e
k
可以表示为 e
k
(i, j, k)。我们记 N
p
为端点数。
图 1.8: 三角形元划分示意图
2. 选择近似函数 在小单元 e
k
(i, j, k) 上进行插值，这里，我们选取线性插值，即
u
h
(x, y) := u
e
(x, y) = ax + by + c
其中：a, b, c 为待定系数。
http://www.ma-xy.com 22 http://www.ma-xy.com
http://www.ma-xy.com
第一章 偏微分方程 1.3 偏微分方程数值方法
3. 求解单元参数 设 u
e
(x, y) 在小单元 e
k
在 i, j, k 三点上的函数值为 u
i
, u
j
, u
k
，则我们可以建
立下列方程组
u(x
i
, y
i
) = ax
i
+ by
i
+ c =: u
i
u(x
j
, y
j
) = ax
j
+ by
j
+ c =: u
j
u(x
k
, y
k
) = ax
k
+ by
k
+ c =: u
k
(1.33)
解上述方程组
(
1.33)
，有
a =
1
2∆e
k
(a
i
u
i
+ a
j
u
j
+ a
k
u
k
)
b =
1
2∆e
k
(b
i
u
i
+ b
j
u
j
+ b
k
u
k
)
c =
1
2∆e
k
(c
i
u
i
+ c
j
u
j
+ c
k
u
k
)
(1.34)
其中：∆e
k
为小单元 e
k
的面积
∆e =
1
2
x
i
y
i
1
x
j
y
j
1
x
k
y
k
1
a
i
=
y
j
1
y
k
1
, a
j
=
y
k
1
y
i
1
, a
k
=
y
i
1
y
j
1
b
i
= −
x
j
1
x
k
1
, b
j
= −
x
k
1
x
i
1
, b
k
= −
x
i
1
x
j
1
c
i
=
x
j
y
j
x
k
y
k
, c
j
=
x
k
y
k
x
i
y
i
, c
k
=
x
i
y
i
x
j
y
j
将 a
1
, a
2
, a
3
带回 u
e
(x, y)，有
u
e
(x, y) = N
i
(x, y)u
i
+ N
j
(x, y)u
j
+ N
k
(x, y)u
k
其中：N
i
(x, y), N
j
(x, y), N
k
(x, y) 为
N
i
(x, y) =
1
2∆e
(a
i
x + b
i
y + c
i
)
N
j
(x, y) =
1
2∆e
(a
j
x + b
j
y + c
j
)
N
k
(x, y) =
1
2∆e
(a
k
x + b
k
y + c
k
)
(1.35)
为了便于书写，我们引入向量
N = [N
i
(x, y) + N
j
(x, y) + N
k
(x, y)]
u
e
= [ u
i
, u
j
, u
k
]
T
http://www.ma-xy.com 23 http://www.ma-xy.com
http://www.ma-xy.com
1.3 偏微分方程数值方法 第一章 偏微分方程
则 u
e
(x, y) 可以写为
u
e
(x, y) = Nu
e
并且 u
e
(x, y) 的梯度可以表示为
∇u
e
=
∂u
e
∂x
∂u
e
∂y
= Bu
e
其中：
B =
∂N
i
∂x
∂N
j
∂x
∂N
k
∂x
∂N
i
∂y
∂N
j
∂y
∂N
k
∂y
=
1
2∆e
k
a
i
a
j
a
k
b
i
b
j
b
k
4. 构造有限维空间的基函数 前面，我们假设函数 u 在单元 e
k
上是线性的，从而构建了函数 u
的有限维空间 U
h
，U
h
是连续分片线性函数集。U
h
中的每一个连续分片线性函数都可以由其基
组 {φ
i
} 表示。这里，我们考虑高度为 1 的“角锥函数”为基函数 φ
i
(i = 1, 2, . . . , N
p
)，这样，U
h
就是一个 N
p
维的有限维空间，其中的函数 u
h
(x, y) ∈ U
h
表示为
u
h
(x, y) =
N
p
i=1
u
i
φ
i
(x, y)
其中，u
i
为待求系数。
5. 构建有限元方程 原方程 (1.32) 的 G 变分问题为：求 u(x, y) 使得
a(u, v) = (f, v), ∀v ∈ U (1.36)
现在，将函数空间 U 变为有限维空间 U
h
，并且还给出了它的基组，我们用 u
h
, v
h
来代替
u, v
u
h
(x, y) =
N
p
i=1
u
i
φ
i
(x, y)
v
h
(x, y) =
N
p
i=1
v
i
φ
i
(x, y)
http://www.ma-xy.com 24 http://www.ma-xy.com
http://www.ma-xy.com
第一章 偏微分方程 1.3 偏微分方程数值方法
则上述无穷维变分问题 (1.36) 变为有穷维变分问题：求 {u
i
}，使得 ∀v
h
∈ U
h
，有
a(u, v) = (f, v), ∀v ∈ U
=⇒ a(u
h
, v
h
) = (f, v
h
), ∀v
h
∈ U
h
=⇒
Ω
∇u
h
∇v
h
dxdy =
Ω
fv
h
dxdy, ∀v
h
∈ U
h
=⇒
N
e
k=1
e
k
∇u
h
· ∇v
h
dxdy =
N
e
k=1
e
k
fv
h
dxdy, ∀v
h
∈ U
h
(1.37)
上面的计算过程是先逐元计算积分，然后再叠加求和。下面，我们就来计算一下逐元积分，以及
它们的和。
6. 单元积分 不失为一般性，我们在单元 e
k
(i, j, k) 上进行分析。设 u, v 函数在三个定点 P
i
, P
j
, P
k
上的数值分别为
u
e
= [ u
i
, u
j
, u
k
]
T
v
e
= [ v
i
, v
j
, v
k
]
T
则在
e
k
上，有
u
e
(x, y) = Nu
e
v
e
(x, y) = Nv
e
∇u
e
= Bu
e
∇v
e
= Bv
e
于是单元 e
k
上的积分
e
k
∇u
h
∇v
h
dxdy 和
e
k
fv
h
dxdy 是可以计算的
e
k
∇u
h
· ∇v
h
dxdy =
e
k
(∇v
h
)
T
(∇u
h
)dxdy
=
e
k
(Bv
e
)
T
(Bu
e
)dxdy
= v
T
e
K
e
u
e
http://www.ma-xy.com 25 http://www.ma-xy.com
http://www.ma-xy.com
1.3 偏微分方程数值方法 第一章 偏微分方程
其中：
K
e
=
e
k
B
T
Bdxdy
= B
T
B
e
k
dxdy
= ∆ e
k
B
T
B
=
1
4∆e
k
a
i
b
i
a
j
b
j
a
k
b
k
a
i
a
j
a
k
b
i
b
j
b
k
∆e
k
为小单元 e
k
的面积。
我们将 K
e
记为
K
e
=
¯
k
ii
¯
k
ij
¯
k
ik
¯
k
ji
¯
k
jj
¯
k
jk
¯
k
ki
¯
k
kj
¯
k
kk
其中：
¯
k
st
=
1
4∆e
k
(a
s
a
t
+ b
s
b
t
) s, t = i, j, k
称 K
e
为单元刚度矩阵。
同理，我们来处理另一个单元积分
e
k
fv
h
dxdy =
e
k
(Nv
e
)
T
fdxdy
=
e
k
v
T
e
N
T
fdxdy
= v
T
e
F
e
其中：向量 N 的计算如公式 (1.35)，
F
e
=
e
k
N
T
fdxdy =
¯
F
i
¯
F
j
¯
F
k
¯
F
j
=
e
k
N
s
fdxdy s = i, j, k
N
i
(x, y) =
1
2∆e
(a
i
x + b
i
y + c
i
)
我们称 F
e
为单元载荷向量。
http://www.ma-xy.com 26 http://www.ma-xy.com
http://www.ma-xy.com
第一章 偏微分方程 1.3 偏微分方程数值方法
7. 单元积分求和 在计算完单元 e
k
上的积分
e
k
∇u
h
∇v
h
dxdy 和
e
k
fv
h
dxdy 之后，我们要
将所有单元的积分相加，即计算
N
e
k=1
e
k
(·)。
对某个小单元 e
k
来说，端点 i, j, k 是逆时针排列的，但是，这只是在一个小单元里的标记，
每个节点在的 N
p
个节点中有着 1 −N
p
的编号。令
v
e
= [ v
i
, v
j
, v
k
]
T
= C
e
v
u
e
= [ u
i
, u
j
, u
k
]
T
= C
e
u
其中：C
e
是一个 3 × N
p
的 0, 1 指示矩阵。C
e
在第一行对应 i 的列、第二行对应 j 的列、第三
行对应 k 的列上的矩阵元素为 1，其余矩阵元素为 0，v = [v
i
, v
j
, v
k
]。
这样，我们有
N
e
k=1
e
k
∇u
h
∇v
h
dxdy
=
N
e
k=1
v
T
e
K
e
u
e
=
N
e
k=1
v
T
C
T
e
K
e
C
e
u
= v
T
Ku (1.38)
其中：
K
=
N
e
k=1
C
T
e
K
e
C
e
。
我们称 K 为总刚度矩阵。同理，可以处理下一个单元累加
N
e
k=1
e
k
f∇v
h
dxdy
=
N
e
k=1
v
T
e
F
e
=
N
e
k=1
v
T
C
T
e
F
e
= v
T
F (1.39)
其中：F =
N
e
k=1
C
T
e
F
e
。我们称 F 为总载荷矩阵。
8. 求解有限元方程组 我们将 (1.38) 和 (1.39) 带入到有限元方程 (1.37)，有
N
e
k=1
e
k
∇u
h
· ∇v
h
dxdy =
N
e
k=1
e
k
fv
h
dxdy, ∀v
h
∈ U
h
=⇒ v
T
Ku = v
T
F, ∀v ∈ R
N
p
=⇒ v
T
(Ku − F) = 0, ∀v ∈ R
N
p
http://www.ma-xy.com 27 http://www.ma-xy.com
http://www.ma-xy.com
1.3 偏微分方程数值方法 第一章 偏微分方程
由此，我们得到 u 的方程组
Ku − F = 0
解上述方程组后即可得到 u，u 的分量 u
1
, u
2
, . . . , u
N
p
即为函数 u
h
(x, y) 在网格点 1 到 N
p
各点
的函数值，即 u(x, y) 的数值估计。之后，我们有函数 u
h
(x, y) =
N
p
i=1
u
i
φ
i
(x, y)，i = 1, 2, . . . , N
p
为定点。
注：在上述求解过程中，我们考虑的是二维泊松方程的第一边值条件，并且假设边值条件为 0，
即 g = 0。关于更一般的边值条件可参考其它书籍，如《偏微分方程数值解法 (第二版)》陆金甫
P228。
考虑：1、上面的讨论是用三角形元对区域 Ω 进行划分的，是否可以考虑其它元，比如：长
方形元。2、上面假设 u(x, y) 在小单元 e
k
上是一次函数，能否考虑高次函数，比如：二次函数。
1.3.7 谱分析方法
有限元法和谱方法的讨论
下面，我们在一维情况下讨论有限元方法和谱方法。考虑如下两个问题
问题 1：回归 (拟合问题)。根据测量数据 {x
i
y
i
}
m
i=1
，求回归函数 y = f(x), a < x < b。
问题 2：常微分方程问题。求函数 y = f(x)，满足如下常微分方程
y
′
= g(y, x)
y(a) = α
y(b) = β
上述两个问题都是在二维平面域 xOy 上求函数 y = f (x)，问题 1 要求 y(x) 距离真实数据
点 {x
i
y
i
}
m
i=1
合理，问题 2 要求 y(x) 满足方程。下面，不要考虑回归，也不要考虑 ODE，我们
来自己找一个函数 y(x)，到平面域上去找，可以有无数种可选择的函数 f
1
, f
2
, . . . , f
n
, . . . 。定义
在区间域 I = [a, b] 内所有函数的集合为 Φ
f
，Φ
f
中包含了各种各样的函数：处处可微的、光滑
的、分段的、处处不可导的、无穷大的。显然，我们不能在无穷维函数空间 Φ
f
中一个一个函数
的找，那太恐怖了，也是不现实的。我们要在 Φ
f
的子空间 (子集 or 子域) 中寻找 f。显然，应该
把那些值无穷的函数类去掉，但仅仅如此还是不够的，我们不妨要求 f 具有特征：光滑或者可积
或者平方可积或者分段或者其它，记这样的特征函数集合为 H = {f|特征}。先不讨论 H，我们
继续看，如果 f 是解，那么我们可以称 H 为特征可行域，但在 H 中找 f 仍然不太现实，例如
在 I 上的平方可积的函数有无穷多个，H 仍然是一个无穷维空间，(f
1
, f
2
, . . . , f
n
, . . . ) = {f }
∞
，
既然如此，我们不妨将空间 H 进一步缩小，比如我们考虑如下的 n 次多项式函数空间 H
n
H
n
=
f
1
(x) = ax + b
f
2
(x) = ax
2
+ bx + c
.
.
.
f
n
(x) = ax
n
+ bx
n−1
+ . . .
(1.40)
http://www.ma-xy.com 28 http://www.ma-xy.com
http://www.ma-xy.com
第一章 偏微分方程 1.3 偏微分方程数值方法
我们就在 H
n
中找 f，看哪一个 f 相对而言最好 (满足一定的指标)。现在好了，H
n
是空间 H
的一个有限维子空间，是可以搜索的，我们称 H
n
为函数 f 的探索空间 (前面记为 U
n
)。对一个
小的 H
n
而言，我们可以将其中的元素 f 列举出来，但对于一个大空间，就无能为力了。考虑到
三维空间中的任一向量可以由基向量表示，那么，我们能不能找到 H
n
中的基函数 {φ
i
}，使 H
n
中的所有函数 f
i
∈ H
n
都可以用基函数表示？可以的，比如上面的 n 次多项式函数集 (1.40)H
n
中的任意一个函数 f
i
都可以用 {x
0
, x
1
, x
2
, . . . , x
n
} 表示
f =
n
i=0
c
i
x
i
其中：c
i
为待求系数。
我们将 {x
0
, x
1
, x
2
, . . . , x
n
} 作为基组，{φ
i
}
n
i=1
= {x
0
, x
1
, x
2
, . . . , x
n
}，并称 H
n
是基组 {φ
i
}
n
i=1
张成的探索空间，记为
H
n
= span{φ
1
, φ
2
, . . . , φ
n
}
其中：n 为基组中基函数的个数。
同一空间中的基组有不同的形式 (即基组不唯一)，并且，观察到三维空间 Oxyz 的基向
量 e
1
, e
2
, e
3
两两之间相互垂直，为正交基向量 (正交基底)，我们也可以在 H
n
中选择正交基组
{φ
i
}
n
i=1
：∀i, j ∈ n，i = j，φ
i
, φ
j
∈ { φ
i
}
n
i=1
，有
b
a
φ
i
φ
j
dx = 0
由此，我们有 H
n
在 I = [a, b] 上的正交基组。
仔细看 f =
n
i=0
c
i
x
i
，这很像傅里叶级数，我们来回顾一下傅里叶级数：将 f(x) 在区间
[−π, π] 上展开，有
ˆ
f(x) =
a
0
2
+
∞
k=1
(a
k
cos kx + b
k
sin kx)
其中：
a
0
=
1
π
π
−π
f(x)dx
a
n
=
1
π
π
−π
f(x) cos nxdx
b
n
=
1
π
π
−π
f(x) sin nxdx
ˆ
f(x) 为函数 f(x) 的逼近 (估计)，f (x) ∼
ˆ
f(x)，当然我们要判断
ˆ
f(x) 的收敛性，这里我们
不做讨论。可以看到 f(x) 在 [−π, π] 上被展成了 {1, cos x, sin x, cos 2x, sin 2x, . . . , cos nx, sin nx}
的线性组合的形式，所以 {1, cos x, sin x, cos 2x, sin 2x, . . . , cos nx, sin nx} 是 H
n
的一个基组，并
http://www.ma-xy.com 29 http://www.ma-xy.com
http://www.ma-xy.com
1.3 偏微分方程数值方法 第一章 偏微分方程
且该基组是正交基组，因为 ∀i, j ∈ n，i = j，有
π
−π
cos(ix) sin(jx)dx = 0
π
−π
cos(ix) cos(jx)dx = 0
π
−π
sin(ix) sin(jx)dx = 0
π
−π
cos(ix)dx = 0
π
−π
sin(ix)dx = 0
即基组中的任意两个不同的基函数的区间卷积为 0。
在这一章后面，我们会讨论一些其它形式的正交基，如：Legendre 多项式、Chebushev 多项
式、Hermite 多项式、Jacobi 多项式、Laguerre 多项式等。
在选取基组 {φ
i
} 时，我们自然希望由 {φ
i
} 张成的空间 H
n
尽可能大。现在，回到我们的问
题当中，设要求解的函数为 u(一维空间为 u(x)，二维空间为 u(x, y))，设其自变量区域为 Ω，区
域边界 ∂Ω，探索空间为 U
n
(前面写为 H
n
)，U
n
的基组为 {φ
i
}
n
i=1
，我们的目标是在 U
n
中寻找
u(x) 的替代 (估计)u
n
(x)。由于 u
n
(x) 可以由基组表示
u
n
(x) =
n
i=1
c
i
φ
i
(x)
那么，寻找 u
n
(x) 的过程就变为寻找系数 c
i
的过程，这样，我们就把泛函极值问题变为了参数
极值问题。
如果对 u(x) 的要求不同 (如 u(x) 不必处处可微)，那么探索空间 U
n
就会不同，从而基组
{φ
i
} 就不同，前面说到，我们希望 U
n
尽可能大。¬一种想法是假设 U
n
(x) 为分段多项式形式，
u
n
(x) ∈ U
n
在 I = [a, b] 上可以是分段一次多项式，也可以是二次和更多次的，分段多项式次
数不同，基组也就不同。不妨设 u
n
(x) 在 I 上为分段一次多项式，将 I 离散化为 n 段，并且
有 n 个基函数 φ
i
(x)。如果基函数 φ
x
(x) 在 x
i
处的函数值为 1，其它基函数在 x
i
处的函数值为
0。那么，基函数的系数 c
i
就可以写为函数值 u
i
，这样，最终求得的 {u
i
} 即为 u(x) 的数值解，
这正是有限元的思想。另一种想法是假设 U
n
(x) 在 Ω 内光滑或者无穷次可微，设探索空间为
U
n
= {u
n
(x)|特征}，基组为 {φ
i
}，则函数 u
n
(x) 可以表示为
u
n
(x) =
n
i=1
c
i
φ
i
(x)
于是求 u
n
(x) 就变为求 c
i
，取不同基组 {φ
i
}，可以得到不同结果，这是谱分析的思想。比较有限
元方法和谱方法，有限元方法以分段函数 (集) 作为探索空间，以分段函数作为基组；谱方法以整
体无限光滑的函数集作为探索空间，以无限光滑函数作为基组 (如 Legendre 多项式、Chebushev
多项式、Hermite 多项式、Jacobi 多项式、Laguerre 多项式等)。
http://www.ma-xy.com 30 http://www.ma-xy.com
http://www.ma-xy.com
第一章 偏微分方程 1.3 偏微分方程数值方法
周期边值问题的谱方法
上面介绍了有限元方法和谱方法的思想，下面，我们来具体介绍谱方法。谱方法和有限元方
法都是基于 R-G 方法的，谱方法有 Galerkin 谱方法、Tau 方法和配点法 (拟谱方法) 之分，在后
面的嫦娥 3 号登月的建模中，我们介绍了配点法，这里就不再重述了。其实，谱方法的思想来源
已久，早在 1820 年，Navier 就用谱方法的思想求解了弹性薄板问题，但是很长一段时间内，谱
方法由于计算量大，在边界处要满足条件等劣势而没有被广泛应用。1965 年，快速傅里叶变换
FFT 的出现使谱方法得以重获新生。近 30 年来，谱方法更是快速发展起来，值得注意的是，谱
方法的基函数是一个定义在 n 维长方体上的正交多项式，因此，谱方法要求求解的区域 Ω 满足
一定的条件。将谱方法用于更广泛的区域 Ω 形式是一项值得研究的工作，
前面，我们提到过，谱方法的探索函数 u
n
(x) 被取为无穷可微函数，Galerkin 谱方法的检验
函数 v
n
(x) 与试探函数 u
n
(x) 属于同一空间 U
n
，并在边界处满足条件，Tau 方法则不要求 v
n
(x)
满足边界条件。用 Galerkin 谱方法求解如下周期边值问题
Lu := −u
′′
= f(x)
u(0) = u(2π) = 0
x ∈ [0, 2π]
(1.41)
用 v 乘上式 (1.41) 两边，并在区间 [0, 2π] 上积分，有
2π
0
Lu vdx =
2π
0
fxdx
然后，对上式左边的积分用分部积分法并结合边值条件，有
2π
0
Lu vdx =
2π
0
u
′
v
′
dx − u
′
v
2π
0
+
2π
0
uvdx
=
2π
0
u
′
v
′
dx
令
a(u, v) =
2π
0
u
′
v
′
dx
(f, v) =
2π
0
fvdx
于是原问题 (1.41) 的 G 变分问题为：求 u(x) ∈ U ，使得对 ∀v ∈ U, v(0) = v(2π) = 0，有
a(u, v) = (f, v) (1.42)
取无穷维空间 U 的有限维探索空间 U
n
= {u
n
(x)|特征}，探索空间的基组为 {φ
k
}
N
k=−N
，基
函数的具体形式为
φ
k
(x) = e
−ikx
− 1 k = ±1, ±2, . . . , ±N
http://www.ma-xy.com 31 http://www.ma-xy.com
http://www.ma-xy.com
1.3 偏微分方程数值方法 第一章 偏微分方程
U
n
由基组张成，U
n
= span{φ
k
}
N
k=−N
。于是，
u
n
(
x
)
∈
U
n
可以表示成如下的形式
u
n
(x) =
N
k=−N
c
k
φ
k
=
N
k=−N
c
k
(e
−ikx
− 1)
于是，原变分问题变为求参数 {c
k
}
N
k=−N
。我们将 u, v 的估计 u
n
(x), v
n
(x)
u
n
(x) =
N
k=−N
c
k
φ
k
=
N
k=−N
c
k
(e
−ikx
− 1)
v
n
(x) =
N
k=−N
c
j
φ
j
=
N
j=−N
c
k
(e
−ijx
− 1)
带入 G 变分问题 (1.42) 中来进行整理，有
a(u, v) = (f, v) ∀v(x) ∈ U
=⇒ a(u
n
, v
n
) = (f, v
n
) ∀v
n
(x) ∈ U
n
(1.43)
由于 v
n
(x) 的任意性 (∀v
n
(x) ∈ U
n
)，不妨令 v
n
(x) 就为空间 U
n
的基函数 φ
j
(x)，于是有下面的
参数方程
a(u, v) = (f, v) ∀v(x) ∈ U
=⇒ a(u
n
, φ
j
) = (f, φ
j
) ∀φ
j
(x) ∈ U
n
, j = ±1, ±2, . . . , ±N
=⇒ a
N
k=−N
c
k
(e
−ikx
− 1), φ
j
= ( f, φ
j
) ∀j = ±1, ±2, . . . , ±N
=⇒
N
k=−N
c
k
a
e
−
ikx
− 1, φ
j
= ( f, φ
j
)
=⇒
N
k=−N
c
k
a
e
−ikx
, φ
j
= ( f, φ
j
)
=⇒
N
k=−N
c
k
a
e
−ikx
, e
−ijx
− 1
= ( f, φ
j
)
=⇒
N
k=−N
c
k
a
e
ikx
, e
ijx
−
N
k=−N
c
k
a
e
−ikx
, 1
= ( f, φ
j
)
=⇒
N
k=−N
c
k
a
e
ikx
, e
ijx
= ( f, φ
j
) ∀j = ±1, ±2, . . . , ±N (1.44)
注意，还有
N
k=−N
c
k
= 0
http://www.ma-xy.com 32 http://www.ma-xy.com
http://www.ma-xy.com
第一章 偏微分方程 1.3 偏微分方程数值方法
我们下面的工作就是求 a
e
ikx
, e
ijx
和 (f, φ
j
)。
a
e
ikx
, e
ijx
=
2π
0
e
ikx
e
ijx
dx
=
2π
0
−(ik)(ij)e
ikx
e
−ijx
dx
=
2πj
2
, k = j
0, k = j
(1.45)
特别地，当 k = j = 0 时，有 a(1, 1) = 0。将上式 (1.45) 带入到参数方程 (1.44)，对 ∀j =
±1, ±2, . . . , ±N，有
a(u
n
, φ
j
) =
N
k=−N
c
k
a
e
ikx
, e
ijx
=: ( f, φ
j
)
=
N
k=−N
c
k
2πj
2
, k = j
N
k=−N
c
k
0, k = j
=: ( f, φ
j
)
= c
j
2π(j
2
) =: (f, φ
j
)
由于 v
n
的任意性 (∀j = ±1, ±2, . . . , ±N )，方程中共 2N + 1 个参数 {c
i
}，所以我们要将 2N
个基函数 φ
j
都带入方程，以形成如下 2N + 1 个方程组
c
j
2π(j
2
) = (f, φ
j
) j = ±1, ±2, . . . , ±N
N
k=−N
c
k
= 0
共有 2N + 1 个待求参数 {c
k
}
N
k=−N
，有 2N + 1 个方程，解之，可得 {c
k
}
N
k=−N
。
正交多项式
上面的谱方法是基于常微周期边界的，对于非周期边界，我们用 Legendre 多项式、Chebyshev
多项式等作为基组来进行求解，为此，我们给出 Legendre、Chebyshev、Hermite、Jacobi 和 Laguerre
等多项式的定义以及部分性质。在介绍多项式之前，先来定义一下基组中基函数的正交性。
正交性 设 f, g ∈ C(I)，ρ 是 I = [a, b] 上的权函数，若
(f, g)
ρ
=
b
a
ρ(x)fgdx = 0
则称 f, g 在 I 上带权 ρ 正交。
http://www.ma-xy.com 33 http://www.ma-xy.com
http://www.ma-xy.com
1.3 偏微分方程数值方法 第一章 偏微分方程
Legendre 多项式 Legendre 多项式是由下列 Legendre 方程导出的
(1 − x
2
)y
′′
− 2xy
′
+ n(n + 1)y = 0
上述方程的通解为
y(x) = AP
n
(x) + BQ
n
(x)
其中：A, B 为常数，P
n
, Q
n
分别为第一、第二类 Legendre 函数。
Legendre 函数的一般形式为
P
n
(x) =
[n/2]
m=0
(−1)
m
(2n − 2m)!
2
n
m!(n − m)!(n − 2m)!
x
n−2m
(1.46)
Q
n
(x) = P
n
(x)
1
(1 − x
2
)P
2
n
dx
其中：[·] 表示取整。
Legendre 函数的递推公式为
P
0
(x) = 1
P
1
(x) = x
P
2
(x) =
1
2
(−1 + 3x
2
)
P
3
(x) =
1
2
(−3x + 5x
2
)
P
4
(x) =
1
8
(3 − 30x
2
+ 34x
4
)
P
5
(x) =
1
8
(15x − 70x
3
+ 63x
5
)
. . .
(n + 1)P
n+1
(x) − (2n + 1)xP
n
(x) + nP
n−1
(x) = 0
n 从 0 到 5 的函数图像如图 (1.9)
所示
https://en.wikipedia.org/wiki/Legendre_polynomials
http://www.ma-xy.com 34 http://www.ma-xy.com
http://www.ma-xy.com
第一章 偏微分方程 1.3 偏微分方程数值方法
图 1.9: Legendre 多项式 MATLAB 图像
Legendre 函数的正交性
1
−1
P
n
(x)P
m
(x)dx =
2
2n + 1
, n = m
0, n = m
=
2
2n + 1
δ
mn
其中：δ
mn
为示性函数，当 m = n 时才起作用。
Chebyshev 多项式 Chebyshev 是下列方程的解
(1 − x
2
)y
′′
− xy
′
+ n
2
y = 0
Chebyshev 函数的递推公式为
T
0
(x) = 1
T
1
(x) = x
T
2
(x) = 2x
2
− 1
T
3
(x) = 4x
3
− 3x
T
4
(x) = 8x
4
− 8x
2
+ 1
. . .
T
n+1
(x) = 2xT
n
(x) − T
n−1
(x) (1.47)
n 取值 0 到 4 时的函数图像如图 (1.10)
®
所示
®
https://en.wikipedia.org/wiki/Chebyshev_polynomials
http://www.ma-xy.com 35 http://www.ma-xy.com
http://www.ma-xy.com
1.3 偏微分方程数值方法 第一章 偏微分方程
图 1.10: Chebyshev 多项式 MATLAB 图像
注意到三角关系式
cos((n + 1)θ) + cos((n −1)θ) = 2 cos θ cos(nθ)
将上式与递推公式 (1.47) 比较，当 θ = arccos x 时，发现 cos(x arccos x) 满足递推关系式，所以
有
T
n
(x) = cos nθ
θ = arccos x
n ⩾ 0
x ∈ I
Chebyshev 的正交性
π
0
cos(mθ) cos(nθ)dx = 0 m = n
或者
(T
m
, T
n
) =
1
−1
1
√
1 − x
2
T
m
(x)T
n
(x)dx = 0 n = m
所以，Cheby 多项式在 [−1, 1] 上带权 ρ(x) =
1
√
1−x
2
正交。
Jacobi 多项式 Chebyshev 多项式和 Legendre 多项式都是 Jacobi 多项式的特殊情况，Jacobi
是下面奇异 Sturm-liourille 问题的特征函数
d
dx
(w
αβ
(x)
d
dx
u(x)) + λ
αβ
n
w
αβ
(x)u(x) = 0
其中：特征值为
λ
αβ
n
= n(n + α + β = 1) n ⩾ 0, α, β ⩾ −1
http://www.ma-xy.com 36 http://www.ma-xy.com
http://www.ma-xy.com
第一章 偏微分方程 1.3 偏微分方程数值方法
w
αβ
(x) 为权重系数。
Jacobi 的表达式为
J
α,β
n
(x) =
Γ(α + n + 1)
n!Γ(α + β + n + 1)
n
m=0
n
m
Γ(α + β + m + n + 1)
Γ(α + m + 1)
x − 1
2
m
其中：Γ(x) 是 Gamma 函数。
Jacobi 的递推关系式为
J
α,β
0
(x) = 1
J
α,β
1
(x) =
1
2
(α + β + 2)x +
1
2
(α − β)
J
α,β
n+1
(x) = (a
α,β
n
x − b
α,β
n
)J
α,β
n
(x) − c
α,β
n
J
α,β
n−1
(x)
其中：系数为
a
α,β
n
=
(2n + α + β + 1)(2n + α + β + 2)
2(n + 1)(n + α + β + 1)
b
α,β
n
=
(β
2
− α
2
)(2n + α + β + 1)
2(n + 1)(n + α + β + 1)(2n + α + β)
c
α,β
n
=
(n + α)(n + β)(2n + α + β + 2))
(n + 1)(n + α + β + 1)(2n + α + β)
Jcacobi 的正交性如下
1
−1
(1 − x)
α
(1 + x)
β
J
α,β
m
(x)J
α,β
n
(x)dx
=
2
α+β+1
2n + α + β + 1
Γ(n + α + 1)Γ(n + β + 1)
Γ(n + α + β + 1)n!
δ
mn
其中：αβ > −1。
Hermite 多项式 Hermite 多项式有两种计算公式，分别为物理学和统计学中的公式，下面介
绍的是物理学公式
H
n
(x) =
[n/2]
m=0
(−1)
m
n!
m!(n − 2m)!
(2x)
n−2m
(1.48)
Hermite 的递推公式为
H
0
(x) = 1
H
1
(x) = 2x
H
2
(x) = 4x
2
− 2
H
3
(x) = 8x
3
− 12x
H
4
(x) = 16x
4
− 48x
2
+ 12
H
5
(x) = 32x
5
− 160x
3
+ 120x
http://www.ma-xy.com 37 http://www.ma-xy.com
http://www.ma-xy.com
1.3 偏微分方程数值方法 第一章 偏微分方程
n 取值从 0 到 5 时的函数图像如图 (1.11)
¯
所示
图 1.11: Hermite 多项式 MATLAB 图像
Hermite
多项式的正交性
∞
−∞
e
−x
2
H
n
(x)H
m
(x)dx = 2
n
n!
√
πδ
mn
即 H
n
(x) 在 [−∞, ∞] 上带权 ρ(x) = e
−x
2
正交。
Laguerre 多项式 Lagurre 是下面方程的解
xy
′′
+ (1 − x)y
′
+ ny = 0
上述方程的广义形式为
xy
′′
+ (α + 1 − x)y
′
+ ny = 0
Lagurre 的表达式为
L
n
(x) =
n
k=0
n
k
(−1)
k
k!
x
k
(1.49)
广义表达式为
L
α
n
(x) =
n
k=0
(−1)
k
n + α
n − k
x
k
k!
¯
https://en.wikipedia.org/wiki/Hermite_polynomials
http://www.ma-xy.com 38 http://www.ma-xy.com
http://www.ma-xy.com
第一章 偏微分方程 1.3 偏微分方程数值方法
Lagurre 的递推公式为
L
0
(x) = 1
L
1
(x) = −x + 1
L
2
(x) =
1
2
(x
2
− 4x + 2)
L
3
(x) =
1
6
(−x
3
+ 9x
2
− 18x + 6)
L
4
(x) =
1
24
(x
4
− 16x
3
+ 72x
2
− 96x + 24)
.
.
.
L
n+1
(x) =
(2n + 1 − x)L
n
(x) − nL
n−1
(x)
n + 1
广义形式的递推公式为
L
α
n
(x) =
α + 1 − x
n
L
α+1
n−1
(x) −
x
n
L
α+2
n−2
(x)
Lagurre 的正交性为
∞
0
x
α
e
−x
L
α
n
(x)L
α
m
(x)dx =
Γ(n + α + 1)
n!
δ
mn
非周期边值问题的谱方法
上面介绍了一些正交多项式，下面，我们以 Legendre 多项式为基组，来求解如下泊松方程
Lu := △u = f(x, y) (x, y) ∈ Ω
u(±1, y) = u(x, ±1) = u
∂Ω
= 0 (x, y) ∈ ∂Ω
其中：Ω = {(x, y)|0 ⩽ x, y ⩽ 1}。
上述问题的 G 变分问题为：求 u(x, y) ∈ U，使得 ∀v(x, y) ∈ U，v(±1, y) = v(x, ±1) = 0，有
a(u, v) = (f, v)
我们取 U 的有限维空间 U
n
= {u
n
(x, y)|特征} 为搜索空间，其基组为 {φ
i
(x, y)}。这里，我们用
Legendre 多项式 L
n
(x), L
m
(y) 作为基函数，则基组为 {L
i
(x), L
j
(y)}
n
i,j=0
，于是，u
n
(x, y) 可以
用基组表示为
u
n
(x, y) =
n
k=0
c
k
φ
k
(x, y) =
n
i=0
n
j=0
a
ij
L
i
(x)L
j
(y) (1.50)
注意：L
i
(x), L
j
(y) 不必满足边界条件，这点与 Galerkin 谱方法不同。
同理，我们可以用基组表示出 v(x, y) 的估计 v
n
(x, y)
v
n
(x, y) =
n
k=0
c
k
φ
k
(x, y) =
n
i=0
n
j=0
b
ij
L
i
(x)L
j
(y) (1.51)
http://www.ma-xy.com 39 http://www.ma-xy.com
http://www.ma-xy.com
1.4 符号和函数空间注记 第一章 偏微分方程
由于 v
n
(x, y) 具有任意性，那么，不妨令
v
ij
= φ
ij
(x, y) = p
i
(x)p
j
(y) =
2i + 1
2
L
i
(x)
2j + 1
2
L
j
(x)
φ
ij
(x, y) 可以取遍所有基组 (i, j = 0, 1, . . . , n −1)，并且注意，只有当 i = j 时，φ
ij
才有值，否
则为 0。我们记边界条件为
B
1
(u) = u(x, −1)
B
2
(u) = u(x, 1)
B
3
(u) = u(−1, y)
B
4
(u) = u(1, y)
那么，a(u, v) = (f, v) 可以表示为
a(u, v) = (f, v) ∀v ∈ U
=⇒ a(u
n
, v
ij
) = (f, v
ij
) ∀v
ij
∈ U
n
=⇒
Ω
∂
2
u
n
∂x
2
+
∂
2
u
n
∂y
2
φ
ij
(x, y)dxdy =
Ω
f(x, y)φ
ij
(x, y)dxdy
上述方程中共有 n + 1 个参数 a
ij
，所以我们要构建方程组。当 i, j = 0, 1, . . . , n − 1 时，可
形成如下方程组
a(u
n
, v
ij
) = (f, v
ij
) i, j = 0, 1, . . . , n − 1
1
−1
B
1
(
u
n
)
p
i
(
x
)
d
x
= 0
i
= 0
,
1
, . . . , n
1
−1
B
2
(u
n
)p
i
(x)dx = 0 i = 0, 1, . . . , n
1
−1
B
3
(u
n
)p
j
(x)dy = 0 j = 0, 1, . . . , n
1
−1
B
4
(u
n
)p
j
(x)dy = 0 j = 0, 1, . . . , n
利用 Legendre 多项式 L
n
(x) 的正交性，解上述方程可以得到解 {a
ij
}
n
i,j=0
。
1.4 符号和函数空间注记
1.4.1 符号注记
PDE 篇常用的数学符号解释
N 表示整数集合
R 表示实数集合
R
n
表示 n 维欧式空间
http://www.ma-xy.com 40 http://www.ma-xy.com
http://www.ma-xy.com
第一章 偏微分方程 1.4 符号和函数空间注记
Ω 表示 R
n
的子集 (域)，是连同开集
∂Ω 表示 Ω 的边界
△ 表示 Laplace 算子，△u =
n
i=1
∂
2
u
∂x
2
i
∇ 表示梯度算子 (哈密顿算子)
C
m
0
(Ω) 表示 Ω 上有紧支集的函数集合
C
∞
0
(Ω) 表示 Ω 上有紧支集的无穷次可微函数的集合
C
∞
(Ω) 表示 Ω 上无穷次可微函数的集合
L
1
loc
(Ω) 表示 Ω 上局部可积函数集，对任意有界可测集 E ⊂ Ω，f 在 E 上 Lebesgue 可积
L
p
(Ω) 表示 Ω 上 p 次 Lebesgue 可积函数集
U 表示函数 u 的求解空间，是 Sobolev 空间
U
n
表示无穷维空间 U 的 n 维有穷维空间
u
n
表示函数 u 的估计，u
n
∈ U
n
H
1
E
是函数集 H
n
≡ U
n
的子集 (子空间)
1.4.2 函数空间注记
距离空间 E 是一个非空集合 (不妨设为函数集合)，如果 ∀x, y ∈ E，∃d(x, y) ∈ R，d 满足：1、
非负 2、对称 3、三角不等式，则称 d 为 x, y 的距离，E 为距离空间。
线性空间 H 是一个非空集合，如果 H 中的元素 x, y 对加法和数量乘法封闭，则称 H 为线性
空间。
线性赋范空间 H 是一个线性空间，H 中的每一元素 x，按某一规则对应一非负实数 ||x||，且
满足
1、||x|| ⩾ 0。如果 ||x|| = 0，则 x = 0，这里的 0 是空间 H 的 0 元素。
2、||αx|| = |α| ||x||，∀α ∈ R, ∀x ∈ H。
3、||x + y|| ⩽ ||x|| + ||y||，∀x, y ∈ H。
则称 H 为赋范线性空间 (线性并且范数存在，线性赋范)，称 ||x|| 为 x 的范数。
Banach 空间 设 H 是一个距离空间，M, E ⊂ H。如果对于 H 中的任意元素 x 的任意邻域均
含有 M 中的点，则称 M 在 E 中稠密；如果 H 中的任一基本序列都有机选，且极限点属于 H，
则称 H 为完备空间。完备的线性赋范空间为 Banach 空间。
L
p
(Ω) 空间 记 Ω 上 Lebesgue 可测且 p 次可积
Ω
|f(x)|
p
dx < +∞ 的函数全体为 L
p
(Ω)。其
中：1 < p < ∞。
若令
||f||
L
p
(Ω)
=
Ω
|f(x)|
p
dx
1
p
为范数，则 L
p
(Ω) = {f | ||f||
L
p
(Ω)
< ∞}。
http://www.ma-xy.com 41 http://www.ma-xy.com
http://www.ma-xy.com
1.4 符号和函数空间注记 第一章 偏微分方程
L
p
(Ω) 是一个 Banch 空间，L
p
(Ω) 空间中有下列重要不等式
1、Minibowski 不等式。如果 1 ⩽ p ⩽ +∞，则 ∀f, g ∈ L
p
(Ω)，有
||f + g||
L
p
(Ω)
⩽ ||f|| + ||g||
2、Holder 不等式。如果 1 ⩽ p, q ⩽ +∞，
1
p
+
1
q
= 1 ，则 ∀f ∈ L
p
(Ω), g ∈ L
q
(Ω)，有 fg ∈∈ L
1
(Ω)
且
||fg||
L
1
(Ω)
⩽ ||f||
L
p
(Ω)
||g||
L
q
(Ω)
3、Schwarz 不等式。如果 p = q = 2，则 ∀f, g ∈ L
2
(Ω)，有 f g ∈ L
1
(Ω) 且
||fg||
L
1
(Ω)
⩽ ||f||
L
2
(Ω)
||g||
L
2
(Ω)
内积空间 设 H 为线性空间，若 ∀x, y ∈ H，∃
x, y
∈ R 与之对应，且满足 1、
x, y
⩾ 0，
∀x, y ∈ H
2、
x + y, z
=
x, z
+
y, z
，∀x, y, z ∈ H
3、
αx, y
= α
x, y
，∀α ∈ R, ∀x, y ∈ H
4、
x, y
=
y, x
，∀x, y ∈ H
则称 < x, y > 为 x, y 的内积，H 为内积空间。在实内积空间中，有
||x|| =
x, y
θ = arccos
x, y
||x|| ||y||
Hilbert 空间 完备的内积空间为 Hilbert 空间，Hilbert 是 Banch 空间。西伯尔特空间为基于
任意正交系上的多项式表示的傅里叶级数和傅里叶变换提供了有效的表达方式。
广义导数 设 u ∈ L
1
loc
(Ω)(即 u 在 Ω 上局部可积)，而且 ∃v ∈ L
1
loc
(Ω)，有
Ω
vφdx = −
Ω
u
∂φ
∂x
i
dx ∀φ ∈ C
∞
0
(Ω)
其中：x
i
为自变量 x 的分量，i = 1, 2, . . . , n。称 v 是 u 对 x
i
的一阶广义导数，记为 v =
∂φ
∂x
i
。
如果 u ∈ L
1
loc
(Ω)，∃v ∈ L
1
loc
(Ω) ，使得
Ω
vφdx = (−1)
|α|
Ω
u∂
α
φdx ∀φ ∈ C
∞
0
(Ω)
则称 v 是 u 的 α 阶广义导数，记为 v = ∂
α
u。
Sobolev 空间 设 m 为非负正数，1 ⩽ p ⩽ +∞，记
S
mp
(Ω) = {u ∈ L
p
(Ω)|∂
α
u ∈ L
p
(Ω), |α| ≤ m}
http://www.ma-xy.com 42 http://www.ma-xy.com
http://www.ma-xy.com
第一章 偏微分方程 1.5 MATLAB 解偏微分方程
并且赋与范数
||u||
m,p
=
|α|≤m
||∂
α
u||
p
L
p
(Ω)
1
p
||u||
m,∞
= max
|α|≤m
||∂
α
u||
0,∞
则称线性赋范空间 S
mp
(Ω) 为 Sobolev 空间。
1.5 MATLAB 解偏微分方程
下面偏微分方程的写法我们用 Laplace 算子和哈密顿算子表示。Laplace 算子和哈密顿算子
在前面介绍过了，这里不再介绍。我们可以将二维泊松方程写为如下三种形式
−
∂
u
∂x
2
+
∂
u
∂y
2
= f(x, y)
− △u = f(x, y)
− ∇ · ∇u = f(x, y)
给出方程的表达式之后，我们需要继续讨论变量的范围以及方程的边界条件。我们仍然沿用
前面的符号，用 Ω 表示 R
n
上的有界区域，用 ∂Ω 表示 Ω 的边界，有了这些之后，我们就可以
求解具体的偏微分方程了。下面，将展示如何运用 MATLAB 来求解 PDE 问题。值得一提的是，
由于二维图像更容易绘制，所以我们往往在 Ω ∈ R
2
上求解函数 u(x, y)，将一维情况 u(x) 扩展
到二维，将三维情况 u(x, y, t) 压缩到二维，下面的示例都是求解 u(x, y) 或者 u(x, t) 的情形。关
于 PDE 的数值方法，前面我们主要介绍了有限差分法、有限元方法和谱方法，下面，我们给出
其 MATLAB 示例。
1.5.1 有限元解法
MATLAB 自带的工具箱 PDEtoolbox 是运用有限元方法求解偏微分问题的，它们是使用最
广泛的方法和工具箱。我们先来看一下 PDEtoolbox 能求解的方程类型和边界条件。
方程类型 PDEtoolbox 能求解的方程类型有如下几种：
1. 椭圆型方程
−∇
(
c
∇
u
) +
au
=
f,
Ω
其中：Ω ∈ R
2
，c, a, f 为 Ω 上的实函数，u 为 Ω 上的待求函数。
2. 抛物型方程
d
∂u
∂t
− ∇(c∇u) + au = f, Ω
3. 双曲型方程
d
∂
2
u
∂t
2
− ∇(c∇u) + au = f, Ω
http://www.ma-xy.com 43 http://www.ma-xy.com
http://www.ma-xy.com
1.5 MATLAB 解偏微分方程 第一章 偏微分方程
4. 特征方程
−∇(c∇u) + au = λdu, Ω
其中：Ω ∈ R
2
，c, a, f 为标量复函数形式系数，可以是时间 t 的函数，λ 为待求本征值，u 为 Ω
上的待求函数。
5. 非线性椭圆方程
−∇(c(u)∇u) + a(u)u = f(u), Ω
其中：c, a, f 可以是 u 的函数。
6. 偏微分方程组
−∇(c
11
∇u
1
) − ∇(c
12
∇u
2
) + a
11
u
1
+ a
12
u
2
= f
1
−∇(c
21
∇u
1
) − ∇(c
22
∇u
2
) + a
21
u
1
+ a
22
u
2
= f
1
边值条件类型 PDEtoolbox 能求解的方程边值条件有如下几种类型：
1. Dirichlet 条件
hu = r, ∂Ω
2. Neumann 条件
n(c∇u) + qu = g, ∂Ω
其中：n 为 ∂Ω 上的单位外法向矢量；g, r, q, h, c 为 ∂Ω 上的函数。
对于特征值问题，要求 g = r = 0；对于非线性问题，g, q, r, h, c 可以是 u 的函数；对于抛物
型和双曲型方程，g, r, h, c, q 可以是时间 t 的函数；对于方程组而言，Dirichlet 条件为
h
11
u
1
+ h
12
u
2
= r
1
h
21
u
1
+ h
22
u
2
= r
2
而 Neumann 条件为
n(c
11
∇u
1
) + n(c
12
∇u
2
) + q
11
u
1
+ q
12
u
2
= g
1
n(c
21
∇u
1
) + n(c
22
∇u
2
) + q
21
u
1
+ q
22
u
2
= g
2
混合边界条件为
h
11
u
1
+ h
12
u
2
= r
n(c
11
∇u
1
) + n(c
12
∇u
2
) + q
11
u
1
+ q
12
u
2
= g
1
+ h
11
µ
n(c
21
∇u
1
) + n(c
22
∇u
2
) + q
11
u
1
+ q
12
u
2
= g
1
+ h
12
µ
其中：µ 的计算要使得 Dirichlet 条件满足。
http://www.ma-xy.com 44 http://www.ma-xy.com
http://www.ma-xy.com
第一章 偏微分方程 1.5 MATLAB 解偏微分方程
有限元步骤 下面，我们用椭圆型方程快速概述一下有限元法的步骤。考虑基本的椭圆方程
− ∇(c∇u) + au = f Ω
hu = r ∂Ω
n(c∇u) + qu = g ∂Ω
¬上述问题的 G 变分问题为：求 u ∈ U，使得 ∀v ∈ U，有
Ω
[−(∇(c∇u)v + auv)]dxdy =
Ω
fvdxdy
对上面的 G 变分问题使用 Green 公式，有
Ω
[(c∇u)∇v + auv]dxdy −
Ω
(−qu + g)vds =
Ω
fvdxdy
设 u(x, y), v(x, y) 可以近似表示为分片多项式的形式。{φ
i
(x, y)} 为空间 U
n
的基组，有
u
n
(x, y) =
n
i=1
φ
i
c
i
v
n
(x, y) =
n
i=1
φ
i
d
i
其中：c
i
为待求系数。
®区域 Ω 的划分。用三角形元将 Ω 划分为 N
e
个单元，N
p
个顶点。一般而言，在使用三角
形元时，要注意：1. 单元顶点不能是相邻单元边上的内点；2. 避免出现大的钝角和大的边；3. 在
u(x, y) 的梯度变化较剧烈的地方，网格要加密；4. 单元编号可以任意，但节点编号应尽量使所有
两个相邻节点编号只差的绝对值中最大值越小越好
min max
|i−j|=1
|i, j|
s.t.
i = 1, 2, . . . , N
p
j = 1, 2, . . . , N
p
¯单元 e 上的多项式。在单元 e 上假设 u
e
(x, y) 为分片一次多项式，即
u
e
(x, y) = ax + by + c
设在 e 的各定点 (i, j, k) 处的函数值为 u
i
, u
j
, , u
k
，记 u
e
= [ u
i
, u
j
, u
k
]
T
，可以计算
u
e
(x, y) = Nu
e
及梯度
∇u
e
(x, y) = Bu
e
°计算单元刚度矩阵，单元载荷向量。
N
e
k=1
e
k
[(c∇u
e
)∇v
e
+ au
e
v
e
]de −
∂Ω
(−qu
e
+ g)v
e
ds =
N
e
k=1
e
k
fv
e
de
http://www.ma-xy.com 45 http://www.ma-xy.com
http://www.ma-xy.com
1.5 MATLAB 解偏微分方程 第一章 偏微分方程
这里，三角单元的梯度 ∇u
e
和面积 S
e
=
e
k
·de 是通过 MATLAB 的 pdetrg 实现的，得到单元
刚度矩阵 K
e
，和单元载荷向量 F
e
。
±总刚度矩阵和总载荷向量。首先，对 F
e
, K
e
进行扩充，MATLAB 通过 eassemple 函数实
现的；然后，再求解总刚度矩阵 K 和总载荷向量 F ，MATLAB 通过 assempde 函数实现。
²对于 Neumann 条件，∂Ω 上不需要满足任何约束条件，由此可以求得
KU = F
1.5.2 MATLAB 应用实例
示例 1：单位圆上的泊松问题 (椭圆型)
− △u = 1
Ω = {(x, y)|x
2
+ y
2
< 1}
u
∂Ω
= 0
此问题的精确解为
u(x, y) =
1 − x
2
− y
2
4
此问题的求解程序如下
1 %%%%%%%%%%%%%%% 示 例 1 ： 单 位 圆 上 的 Poisson 方 程 %%%%%%%%%%%%%%%%
2 %% 旧 方 法
3 %−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
4 %一 般 方 程 形 式 为 : −div ( c*grad (u) ) + a*u = f ，
5 %我 们 要 解 的 方 程 为 ： c = 1 ; a = 0; f = 1 ;
6 %−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
7 clc , c l e a r
8 g = @ c i r c l e g ; b = @ circ leb 1 ;
9 c = 1; a = 0 ; f = 1;
10 [ p , e , t ] = i nitmesh (g , ’hmax ’ , 1) ; %初 始 化 网 格 。 x = p ( 1 , : ) ; y = p ( 2 , : )
11 % pdemesh (p , e , t );%绘 制 网 格 , p 的 第 一 行 和 第 二 行 为 x , y坐 标 ， e 是 边 界 矩 阵
12 e r ro r = [ ] ;
13 e r r = 1 ;
14 while e r r > 0. 0 01 ;
15 [ p , e , t ] = ref inemesh (g , p , e , t ) ; %加 密 网 格
16 u = assempde (b , p , e , t , c , a , f ) ; %求 解u。 组 成 偏 微 分 方 程 刚 度 矩 阵 及 其 右 边
17 exact = (1−p ( 1 , : ) . ^ 2 − p ( 2 , : ) .^ 2 ) ’/ 4 ; %精 确 解
18 e r r = norm(u − exact , ’ i n f ’ ) ; %计 算 误 差
19 er ro r = [ er r or , e r r ] ;
20 end
21 f i g u r e
22 pdemesh (p , e , t ) ;%绘 制 网 格
23 f i g u r e
24 pdesur f (p , t , u) ;%绘 制 解 的 曲 面
25 f i g u r e
http://www.ma-xy.com 46 http://www.ma-xy.com
http://www.ma-xy.com
第一章 偏微分方程 1.5 MATLAB 解偏微分方程
26 pdesur f (p , t , u − exact ) ;%绘 制 误 差 图
27 %% 新 方 法
28 %−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
29 %一 般 方 程 形 式 为 ：m*Utt + d*Ut − div ( c*grad (u) ) + a*u = f ，
30 %我 们 要 解 的 方 程 为：− div ( c*grad (u ) ) = 1 .
31 %−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
32 clc , c l e a r
33 % 确 定 方 程 （ 模 型 系 数 及 边 界 形 状 ）
34 m = 0 ; d = 0 ; c = 1 ; a = 0 ; f = 1 ;
35 g = @ c i r c l e g ;%边 界 形 状 g ， @lshapeg
36 f i g u r e
37 pdegplot (g , ’ EdgeLabels ’ , ’ on ’ ) ;%查 看 边 界 形 状g ， 并 打 开 边 界 标 号
38 ax i s equal
39 % 创 建 一 个PDE模 型−−pdemodel
40 numberOfPDE = 1;%偏 微 分 方 程 （ 组 ） 的 个 数
41 pdemodel = createp de (numberOfPDE) ;%返 回 一 个n 个 方 程 的 偏 微 分 方 程 模 型 容 器 。
42 geometryFromEdges ( pdemodel , g ) ;%创 建 一 个 几 何 实 体 ， oumiga
43 s p e c i f y C o e f f i c i e n t s ( pdemodel , ’m’ ,m, ’ d ’ ,d , ’ c ’ , c , ’ a ’ , a , ’ f ’ , f ) ;%指 定 模 型 的 系 数
44 applyBoundaryCondition ( pdemodel , ’Edge ’ , ( 1 : 4 ) , ’ u ’ ,0 ) ;%设 置pdem的 边 界 条 件
45 hmax = 1 ;
46 generateMesh (pdemodel , ’Hmax ’ ,hmax) ;%形 成 单 元 e
47 %
48 f i g u r e
49 pdemesh ( pdemodel ) ;
50 ax i s equal
51 er = I nf ;
52 while er > 0.001
53 hmax = hmax/ 2 ;
54 generateMesh ( pdemodel , ’Hmax’ ,hmax) ;
55 r e s u l t = s olvepde ( pdemodel ) ; %解pdem
56 msh = pdemodel . Mesh ;% 第 一 行 为x 坐 标 ， 第 二 行 为 y坐 标
57 u = r e s u l t . NodalSolution ;
58 exact = (1 − msh. Nodes ( 1 , : ) . ^2 − msh. Nodes ( 2 , : ) .^ 2 ) ’/ 4 ; %解 析 解
59 er = norm(u−exact , ’ i n f ’ ) ;
60 f p r i n t f ( ’ Error : %e . Number of nodes : %d\n ’ , er , s i z e (msh. Nodes , 2 ) ) ;
61 end
62 f i g u r e
63 pdemesh ( pdemodel ) ;
64 ax i s equal
65 f i g u r e
66 pdeplot (pdemodel , ’XYData ’ ,u , ’ZData ’ , u , ’ ColorBar ’ , ’ o f f ’ )
67
新旧方法求解结果如图 (1.12) 所示
http://www.ma-xy.com 47 http://www.ma-xy.com
http://www.ma-xy.com
1.5 MATLAB 解偏微分方程 第一章 偏微分方程
(a) 旧方法 (b) 新方法
图 1.12: 单位圆上的泊松问题
示例 2：单位正方形上的波动方程 (双曲型)
∂
2
u
∂t
2
− ∇ · ∇u = 0
Ω = {(x, y)|0 ⩽ x, y ⩽ 1}
初值条件：
u(x, 0) = arctan
cos
πx
2

∂u
∂t
t=0
= 3 sin(πx)e
sin
(
πy
2
)
边界条件：
u
x=±1
= 0 2,4 边为 Dirichlet 条件
∂u
∂n
y =±1
= 0 1,3 边为 Neumann 条件
求解程序如下
1 %%%%%%%%%%%%%%%%%%%%%%%%% 示 例 2 : 波 动 方 程 %%%%%%%%%%%%%%%%%%%
2 %−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
3 %波 动 方 程 的 一 般 形 式 为 ： d*Utt − div ( c*giadU ) + a*U = f
4 %我 们 要 解 决 的 二 维 波 动 方 程 为 ： Utt − div ( c *giadU ) = 0；
5 %−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
http://www.ma-xy.com 48 http://www.ma-xy.com
http://www.ma-xy.com
第一章 偏微分方程 1.5 MATLAB 解偏微分方程
6 %% 旧 方 法
7 clc , c l e a r
8 c = 1; a = 0 ; f = 0; d = 1 ;
9 g = @squareg ;%边 界 区 域
10 b = @squareb3 ;
11 [ p , e , t ] = i nitmesh (g ) ;%初 始 化 网 格
12 x = p( 1 , : ) ’ ;
13 y = p( 2 , : ) ’ ;
14 u0 = atan ( co s ( p i /2*x ) ) ;%初 始 条 件
15 ut0 = 3* s i n ( pi *x ) . * exp ( s i n ( pi /2*y) ) ;%初 始 条 件
16 n = 31 ;
17 t l i s t = l i n s p a c e ( 0 , 5 , n) ;%离 散 时 间
18 uu = hy p e rboli c ( u0 , ut0 , t l i s t , b , p , e , t , c , a , f , d) ;%求 解 双 曲 线 型PDE方 程 ， uu 是 一
个 31 列 的 矩 阵 ， 每 个 时 间 t 为 一 列
19 del t a = −1: 0 . 1 :1;
20 [ uxy , tn , a2 , a3 ] = t r i 2 g r i d (p , t , uu ( : , 1) , delt a , delta ) ;%将 三 角 形 网 格 转 化 为 矩 形 网
格
21 gp = [ tn ; a2 ; a3 ] ;
22 umax = max(max(uu) ) ;
23 umin = min(min(uu) ) ;
24 newplot
25 M = moviein (n ) ;
26 f o r i = 1 : n
27 pdeplot ( p , e , t , ’ xydata ’ , uu ( : , i ) , ’ zdata ’ , uu ( : , i ) , . . .
28 ’mesh ’ , ’ o f f ’ , ’ xygrid ’ , ’ on ’ , ’ gridparam ’ , gp , . . .
29 ’ colorbar ’ , ’ o f f ’ , ’ z s t y l e ’ , ’ continuous ’ ) ;
30 a x is ([ −1 , 1 , −1, 1 , umin , umax] ) ;
31 c axi s ( [ umin , umax ] ) ;
32 M( : , i ) = getframe ;
33 end
34 movie (M, 10) ;
35 %% 新 方 法
36 %−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
37 %方 程 的 一 般 形 式 为 ：m*Utt + d*Ut − div ( c*giadU ) + a*U = f
38 %我 们 要 解 决 的 二 维 波 动 方 程 为 ： Utt − div ( c *giadU ) = 0；
39 %−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
40 clc , c l e a r
41 % 1 、 确 定 方 程 （ 模 型 系 数 ， 边 界 形 状 ）
42 m = 1 ; d = 0 ; c = 1 ; a = 0 ; f = 0 ;
43 g = @squareg ;
44 numberOfPDE = 1;
45 model = createp de (numberOfPDE) ;
46 geometryFromEdges ( model , g ) ;
47 % pdegplot (model , ’ edgeLabels ’ , ’ on ’ ) ;
48 % ylim ([ −1.1 1 . 1 ] ) ;
49 % axi s equal , t i t l e ’Geometry With Edge Labels Displayed ’ , x l a b el x , y la b e l y
50 %2 、 确 定 边 界 条 件
51 s p e c i f y C o e f f i c i e n t s ( model , ’m’ ,m, ’d ’ , d , ’ c ’ , c , ’ a ’ , a , ’ f ’ , f ) ;
52 applyBoundaryCondition ( model , ’Edge ’ , [ 2 , 4 ] , ’u ’ , 0) ;%d i r i c h l e t 边 界 条 件
53 applyBoundaryCondition ( model , ’Edge ’ , [ 1 3] , ’ g ’ , 0) ;%neumann边 界 条 件
54 %3 、 形 成 最 终 模 型
55 generateMesh (model ) ;%形 成 单 元 e
56 u0 = @( l o c a t i o n ) atan ( cos ( pi /2* lo c a t i o n . x ) ) ;%初 始 条 件
http://www.ma-xy.com 49 http://www.ma-xy.com
http://www.ma-xy.com
1.5 MATLAB 解偏微分方程 第一章 偏微分方程
57 ut0 = @( lo c a t i o n ) 3* si n ( pi * l o c a t i o n . x ) . * exp ( s i n ( pi /2* l o c a t i on . y ) ) ;%初 始 条 件
58 s e t I n i t i a l C o n d i t i o n s ( model , u0 , ut0 ) ;%设 置 初 始 条 件
59 %4 、 求 解 模 型
60 n = 31 ;
61 t l i s t = l i n s p a c e ( 0 , 5 ,n) ;
62 model . SolverOptions . R e po r t S t a t i s t i c s = ’ on ’ ;
63 r e s u l t = s olvepde (model , t l i s t ) ;%解PDE模 型
64 u = r e s u l t . NodalSolution ;%u是 一 个 31 列 的 矩 阵 ， 每 个 时 间 t 为 一 列
65 f i g u r e
66 umax = max(max(u) ) ;
67 umin = min(min(u ) ) ;
68 f o r i = 1 : n
69 pdeplot ( model , ’ xydata ’ ,u ( : , i ) , ’ zdata ’ ,u ( : , i ) , ’ z s t y l e ’ , ’ continuous ’ , . . .
70 ’mesh ’ , ’ o f f ’ , ’ xygrid ’ , ’ on ’ , ’ colorbar ’ , ’ o f f ’ ) ;
71 a x is ([−1 1 −1 1 umin umax] ) ;
72 c axi s ( [ umin umax ] ) ;
73 x la b e l x , yl ab e l y , z l a b e l u
74 M( i ) = getframe ;
75 end
76 movie (M, 1) ;
77
用新方法求解上述问题，最终停止时的结果如图 (1.13) 所示
图 1.13: 单位正方形上的波动方程
示例 3：薄板中的非线性传热 (抛物型方程)
一般形式的热传导方程如下
ρc
∂T
∂t
− ∇(k∇T ) = Q + h(T
Q
− T )
其中：ρ 为密度，c 为比热容，k 为导热系数，Q 为热源，h 为对流热的传输系数，T
Q
为环境温
度，h(T
Q
− T ) 表示从环境向区域内部传输的热量。
如果仅考虑稳态情况
∂T
∂t
= 0 ，有
−∇(k∇T ) = Q + h(T
1
− T )
http://www.ma-xy.com 50 http://www.ma-xy.com
http://www.ma-xy.com
第一章 偏微分方程 1.5 MATLAB 解偏微分方程
令
Q
c
= h
c
(T − T
a
)
其中：T
a
是环境温度，h
c
是热对流系数。
从单位面积辐射的热转移量为
Q
r
= ϵσ(T
4
− T
4
a
)
其中：ϵ 为薄板发射率，σ 是波尔兹曼常数。
ρct
2
∂T
∂t
− kt
2
∇
2
T = −2Q
c
− 2Q
r
其中：ρ 为密度，c 为比热容，t
2
为薄板厚度。
ρct
2
∂T
∂t
− kt
2
∇
2
T + 2h
c
T + 2ϵσT
4
= 2 h
c
T
a
+ 2ϵσT
4
a
上述问题的求解程序如下
1 %%%%%%%%%%%%%%% 示 例 3 ： 金 属 板 的 热 传 导 问 题 %%%%%%%%%%%%%%%%
2 %% 旧 方 法
3 %−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
4 %抛 物 线 方 程 的 一 般 形 式 ： d*Ut − div ( c *gradU ) + a*U = f
5 %带 求 解 的 抛 物 线 方 程 为 ： d*Ut − div ( c *gradU ) = 0
6 %−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
7 clc , c l e a r
8 % 设 置 方 程 系 数
9 c = 1; a = 0 ; f = 1; d = 1 ;
10 %设 置 区 域
11 R1 = [3; 4; −1;1 ;1; −1; −1; −1;1; 1];
12 C1 = [ 1 ; 0 ; 0 ; 0 . 4 ] ;
13 C1 = [ C1; z e r o s ( l eng th (R1) − len gth (C1) ,1 ) ] ;
14 gd = [ R1, C1 ] ;
15 s f = ’R1+C1 ’ ;
16 ns = char ( ’R1 ’ , ’C1 ’ ) ’ ;
17 g = decsg (gd , sf , ns ) ;
18 %
19 numberOfPDE = 1;
20 model = createp de (numberOfPDE) ;%创 建 一 个PDEmodel
21 geometryFromEdges ( model , g ) ;%创 建 二 维 几 何
22 applyBoundaryCondition ( model , ’ Edge ’ , ( 1 : 4 ) , ’u ’ ,0 ) ;%边 界 条 件
23 s e t I n i t i a l C o n d i t i o n s ( model , 0 ) ;%初 始 条 件 ， 圆 ( Face2 ) 内 为1 ， 圆 外 为 0 。
24 s e t I n i t i a l C o n d i t i o n s ( model , 1 , ’ Face ’ ,2 ) ;%初 始 条 件
25 msh = generateMesh (model) ;%有 限 元 划 分
26 %
27 nframes = 20;
28 t l i s t = l i n s p a c e ( 0 ,0 . 1 , nframes ) ;%设 置 离 散 时 间
29 model . SolverOptions . R e po r t S t a t i s t i c s = ’ on ’ ;
30 u = parabolic ( 0 , t l i s t , model , c , a , f , d) ;%用 pa r a b o l i c 求 解 抛 物 型 方 程
31 %
32 f i g u r e
33 umax = max(max(u) ) ;
34 umin = min(min(u ) ) ;
http://www.ma-xy.com 51 http://www.ma-xy.com
http://www.ma-xy.com
1.5 MATLAB 解偏微分方程 第一章 偏微分方程
35 f o r j = 1 : nframes ,
36 pdeplot ( model , ’XYData ’ ,u ( : , j ) , ’ ZData ’ ,u ( : , j ) ) ;
37 c axi s ( [ umin umax ] ) ;
38 a x is ([−1 1 −1 1 0 1 ] ) ;
39 M( j ) = getframe ;
40 end
41 %% 新 方 法
42 %−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
43 %方 程 的 一 般 形 式 ：m*Utt + d*Ut − div ( c*gradU ) + a*U = f
44 %带 求 解 的 抛 物 线 方 程 为 ： d*Ut − div ( c *gradU ) = 0
45 %−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
46 clc , c l e a r
47 m = 0 ; d = 1 ; c = 1 ; a = 0 ; f = 1 ;
48 %
49 R1 = [3; 4; −1;1 ;1; −1; −1; −1;1; 1];
50 C1 = [ 1 ; 0 ; 0 ; 0 . 4 ] ;
51 C1 = [ C1; z e r o s ( l eng th (R1) − len gth (C1) ,1 ) ] ;
52 gd = [ R1, C1 ] ;
53 s f = ’R1+C1 ’ ;
54 ns = char ( ’R1 ’ , ’C1 ’ ) ’ ;
55 g = decsg (gd , sf , ns ) ;
56 %
57 numberOfPDE = 1;
58 model = createp de (numberOfPDE) ;
59 geometryFromEdges ( model , g ) ;
60 applyBoundaryCondition ( model , ’ Edge ’ , ( 1 : 4 ) , ’u ’ ,0 ) ;
61 s p e c i f y C o e f f i c i e n t s ( model , ’m’ ,m, . . .
62 ’d ’ , d , . . .
63 ’ c ’ ,c , . . .
64 ’ a ’ , a , . . .
65 ’ f ’ , f ) ;
66 s e t I n i t i a l C o n d i t i o n s ( model , 0 ) ;
67 s e t I n i t i a l C o n d i t i o n s ( model , 1 , ’ Face ’ ,2 ) ;
68 msh = generateMesh (model) ;
69 nframes = 20;
70 t l i s t = l i n s p a c e ( 0 ,0 . 1 , nframes ) ;
71 model . SolverOptions . R e po r t S t a t i s t i c s =’ on ’ ;
72 r e s u l t = s olvepde (model , t l i s t ) ;
73 u1 = r e s u l t . NodalSolution ;
74 %
75 f i g u r e
76 umax = max(max(u1) ) ;
77 umin = min(min( u1) ) ;
78 f o r j = 1 : nframes ,
79 pdeplot ( model , ’XYData ’ ,u1 ( : , j ) , ’ ZData ’ ,u1 ( : , j ) ) ;
80 c axi s ( [ umin umax ] ) ;
81 a x is ([−1 1 −1 1 0 1 ] ) ;
82 Mv( j ) = getframe ;
83 end
84
最终停止时的结果如图 (1.14) 所示
http://www.ma-xy.com 52 http://www.ma-xy.com
http://www.ma-xy.com
第一章 偏微分方程 1.5 MATLAB 解偏微分方程
图 1.14: 薄板中的非线性传热
http://www.ma-xy.com 53 http://www.ma-xy.com