偏微分方程数值解程序设计与实现—

常用算子符号

梯度算子
Rd\mathbb{R}^dRd空间中标量函数u(x)u(\bf{x})u(x)，其梯度算子定义如下：
gradu(x)=∇u(x)=[∂u∂x0∂u∂x1⋮∂u∂xd−1]grad u(\mathbf{x})=\nabla u(\mathbf{x})= \begin{bmatrix} \frac{\partial u}{\partial x_0} \\ \frac{\partial u}{\partial x_1} \\ \vdots \\ \frac{\partial u}{\partial x_{d-1}} \end{bmatrix} gradu(x)=∇u(x)=⎣⎢⎢⎢⎡∂x0∂u∂x1∂u⋮∂xd−1∂u⎦⎥⎥⎥⎤
梯度的几何意义：标量函数在一点的梯度向量，指向该点处函数值增加最快的方向，且长度是该点处沿这个方向的函数变化率。
散度算子
Rd\mathbb{R}^dRd空间中标量函数u(x)=[u0(x),u1(x),⋯,ud−1(x)]T\mathbf{u}(\mathbf{x})=[u_0(\mathbf{x}),u_1(\mathbf{x}), \cdots,u_{d-1}(\mathbf{x})]^Tu(x)=[u0(x),u1(x),⋯,ud−1(x)]T，其散度算子定义如下：
divu(x)=∇⋅u(x)=∂u0∂x0+∂u1∂x1+⋯+∂ud−1∂xd−1div \mathbf{u}(\mathbf{x})=\nabla \cdot \mathbf{u}(\mathbf{x})= \frac{\partial u_0}{\partial x_0} + \frac{\partial u_1}{\partial x_1}+\cdots +\frac{\partial u_{d-1}}{\partial x_{d-1}} divu(x)=∇⋅u(x)=∂x0∂u0+∂x1∂u1+⋯+∂xd−1∂ud−1
散度的物理意义：假设u(x)\mathbf{u}(\bf{x})u(x)代表一个稳定流动的不可压缩流体（密度为1）的速度场，他在x\bf{x}x处的散度，就是穿出单位体积边界的通量，也叫通量密度，而通量是流场单位时间内沿指定侧通过的一个曲面的量。给定一个位于封闭曲面S内点x\bf{x}x,S围成的区域记为Ω\OmegaΩ,
∇⋅u(x)=limS→x1∣Ω∣∫Su⋅nds\nabla \cdot \mathbf{u}(\mathbf{x})=lim_{S\rightarrow\mathbf{x}}\frac{1}{|\Omega|}\int_S \mathbf{u}\cdot \mathbf{n} d\mathbf{s} ∇⋅u(x)=limS→x∣Ω∣1∫Su⋅nds
旋度算子
R3\mathbb{R}^3R3空间中标量函数u(x)=[u0(x),u1(x),u2(x)]T\mathbf{u}(\mathbf{x})=[u_0(\mathbf{x}),u_1(\mathbf{x}), u_{2}(\mathbf{x})]^Tu(x)=[u0(x),u1(x),u2(x)]T，其旋度算子定义如下：
curlu(x)=∇×u(x)=[∂u2∂y−∂u1∂z∂u0∂z−∂u2∂x∂u1∂x−∂u0∂y]curl u(\mathbf{x})=\nabla \times u(\mathbf{x})= \begin{bmatrix} \frac{\partial u_2}{\partial y}-\frac{\partial u_1}{\partial z} \\ \frac{\partial u_0}{\partial z}-\frac{\partial u_2}{\partial x} \\ \frac{\partial u_1}{\partial x}-\frac{\partial u_0}{\partial y} \end{bmatrix} curlu(x)=∇×u(x)=⎣⎢⎡∂y∂u2−∂z∂u1∂z∂u0−∂x∂u2∂x∂u1−∂y∂u0⎦⎥⎤
散度的物理意义：向量场u\mathbf{u}u在点x\bf{x}x处的旋度向量，用来描述流体以x\bf{x}x为中心的漩涡强度和方向，它的指向与漩涡旋转最快的方向满足右手法则，长度为旋转最快的方向的角速度的两倍。
Laplace算子
Rd\mathbb{R}^dRd空间中标量函数u(x)u(\bf{x})u(x)的梯度的散度给出的微分算子称为Laplace算子，通常写成
Δu(x)=∇2u(x)=∇⋅∇u(x)=∑i=0d−1∂2u∂xi2\Delta u(\mathbf{x})=\nabla^2u(\mathbf{x})=\nabla\cdot \nabla u(\mathbf{x})=\sum_{i=0}^{d-1}\frac{\partial^2u}{\partial x^2_i} Δu(x)=∇2u(x)=∇⋅∇u(x)=i=0∑d−1∂xi2∂2u

散度定理

给定定义在Ω∈Rd\Omega \in\mathbb{R}^dΩ∈Rd向量函数F(x)\bf{F}(\bf{x})F(x)，记n\bf{n}n为边界∂Ω\partial \Omega∂Ω上的单位外法线向量，则有
∫Ω∇⋅Fdx=∫∂ΩF⋅nds\int_{\Omega}\nabla \cdot \mathbf{F}d\mathbf{x}=\int_{\partial\Omega} \mathbf{F}\cdot \mathbf{n}d\mathbf{s} ∫Ω∇⋅Fdx=∫∂ΩF⋅nds
这就是散度定理。

散度定理的应用

1）取F=v∇u\mathbf{F}=v\nabla uF=v∇u，可得：
∫Ω∇⋅(v∇udx)=∫∂Ωv∇u⋅nds\int_{\Omega}\nabla \cdot (v\nabla ud\mathbf{x})=\int_{\partial\Omega} v\nabla u\cdot \mathbf{n}d\mathbf{s} ∫Ω∇⋅(v∇udx)=∫∂Ωv∇u⋅nds

∫ΩvΔudx+∫Ω∇u⋅∇vdx=∫∂Ωv∇u⋅nds\int_{\Omega}v\Delta ud\mathbf{x}+\int_{\Omega}\nabla u \cdot \nabla vd\mathbf{x}=\int_{\partial\Omega} v\nabla u\cdot \mathbf{n}d\mathbf{s} ∫ΩvΔudx+∫Ω∇u⋅∇vdx=∫∂Ωv∇u⋅nds
2）取F=[u,0,0]T\mathbf{F}=[u,0,0]^TF=[u,0,0]T，可得：
∫Ω∇⋅[u00]dx=∫Ωux=∫∂Ωunxds\int_{\Omega}\nabla \cdot \begin{bmatrix} u\\0\\0 \end{bmatrix}d\mathbf{x}=\int_{\Omega}u_x=\int_{\partial\Omega} u \mathbf{n}_xd\mathbf{s} ∫Ω∇⋅⎣⎡u00⎦⎤dx=∫Ωux=∫∂Ωunxds
3）

函数空间

这里所用到的函数空间，简单来讲就是具有某些共同性质，且对于线性运算封闭的函数组成集合，它里面通常有无穷多个函数。

线性运算封闭性：对于函数空间V中任意两个函数vvv,www，及实数空间R\mathbb{R}R中任意两个常数c0c_0c0和c1c_1c1，满足
c0v+c1w∈Vc_0v+c_1w\in V c0v+c1w∈V
函数空间的基：设S=ϕiS={\phi_i}S=ϕi是VVV的一个子集，如果SSS中的元素线性无关，且VVV中的任意一个函数vvv都可以由它们线性表出，即存在一组和SSS中函数一样多的常数集合ci{c_i}ci，使得
v=∑i=0ciϕiv=\sum_{i=0}c_i\phi_i v=i=0∑ciϕi
则称ϕi{\phi_i}ϕi为VVV的一组基。
(1)函数空间由基唯一决定，即给定一组基，就唯一张成一个空间V=span{ϕi}V=span\{\phi_i\}V=span{ϕi}
(2)函数空间的基不是唯一的，问题是编程实现的时候你要选择哪一组？
(3)给定空间的一组基，就可以建立起空间VVV到空间Rd\mathbb{R}^dRd的一一映射，其中ddd是向量空间的维数，可以为无穷
v∈V↔[c0,c1,⋯]∈Rdv\in V \leftrightarrow [c_0,c_1,\cdots]\in\mathbb{R}^d v∈V↔[c0,c1,⋯]∈Rd
函数空间的维数：函数空间基函数的个数就成为函数空间的维数。包括有限维空间和无限维空间。编程实现只能处理有限位的函数空间！
常见的函数空间：
记Ω\OmegaΩ为Rd\mathbb{R}^dRd空间中的任意子区域
1）L2(Ω)L^2(\Omega)L2(Ω)：对于任意的$V \in L2(Ω)L^2(\Omega)L2(Ω)，∫Ωv2dx<∞\int_{\Omega}v^2d\mathbf{x}<\infty∫Ωv2dx<∞。
2）H1(Ω)H^1(\Omega)H1(Ω)：函数及其导数都属于L2(Ω)L^2(\Omega)L2(Ω)空间。
3）H(div,Ω)\mathbf{H}(div,\Omega)H(div,Ω)：向量函数空间，向量函数的每个分量及其散度都属于L2(Ω)L^2(\Omega)L2(Ω)空间。
4）H(curl,Ω)\mathbf{H}(curl,\Omega)H(curl,Ω)：向量函数空间，向量函数的每个分量及其旋度都属于L2(Ω)L^2(\Omega)L2(Ω)空间。
5）Pk(Ω)\mathbb{P}_k(\Omega)Pk(Ω)：不大于kkk次的多项式函数组成的空间。
6）Pk(Ω；Rd)\mathbb{P}_k(\Omega；\mathbb{R}^d)Pk(Ω；Rd)：每个分量都是不大于kkk次的多项式函数的向量函数空间。
多项式函数空间是编程实现有限元算法的基础空间。

偏微分方程数值解：从无限到有限

“Hello World"模型：Possion方程

给定区域Ω⊂Rd\Omega\subset\mathbb{R}^dΩ⊂Rd，其边界记为∂Ω\partial \Omega∂Ω，求标量函数u(x)u(\mathbf{x})u(x)，满足
−Δu(x)=f(x),∀x∈Ωu(x)=0,∀x∈∂Ω-\Delta u(\mathbf{x})=f(\mathbf{x}), \forall \mathbf{x}\in \Omega\\ u(\mathbf{x})=0,\forall \mathbf{x}\in \partial\Omega −Δu(x)=f(x),∀x∈Ωu(x)=0,∀x∈∂Ω

Laplace算子Δ\DeltaΔ的定义
- d=1d=1d=1：Δ(x)=Δu(x)=ux\Delta (\mathbf{x})=\Delta u(x)=u_xΔ(x)=Δu(x)=ux
- d=2d=2d=2：Δ(x)=Δu(x,y)=uxx+uyy\Delta (\mathbf{x})=\Delta u(x,y)=u_{xx}+u_{yy}Δ(x)=Δu(x,y)=uxx+uyy
- d=3d=3d=3：Δ(x)=Δu(x,y,z)=uxx+uyy+uzz\Delta (\mathbf{x})=\Delta u(x,y,z)=u_{xx}+u_{yy}+u_{zz}Δ(x)=Δu(x,y,z)=uxx+uyy+uzz
未知量是什么？需要求无穷多个点处对应的函数值。
该问题求解的困难之处：算子和未知量的无限性。

从微分到积分：连续的弱形式

从无限到有限：离散的弱形式

从无限到有限：数值积分

总结

上面通过弱形式（变分）把偏微分方程转化为代数方程的方法统称为Galerkin方法。
偏微分数值计算方法就是为解决偏微分方程模型的无限性和计算资源的有限性这对本质矛盾而出现，并不断前进发展的。
与Galerkin方法不同的另一种方法——有限差分方法，是通过数值微分来解决微分算子定义的无限性问题。

全部内容来自湘潭大学魏华祎老师暑期学校课程https://www.bilibili.com/video/BV1np4y1S7LK/