PDE的数值解法（有限元，有限差分法）综合介绍

以下内容均可参考本人知乎文章添加链接描述和添加链接描述

有限差分法

有限差分法finite difference(FD)是求解微分方程的最为容易理解的方法，下面将针对几类常见的PDE来做一些具体的介绍。由于本人知识有限，关于误差分析和收敛性证明都不会介绍.

一维例子

我们以一个一维PDE的求解来介绍有限差分算法，这里给出精确解$\hat{u}=x(1-x)\exp (x)$.

{ − u x x = f , x ∈ ( 0 , 1 ) , f = exp ⁡ ( x ) ( − x 2 − x + 1 ) , u ( 0 ) = 0 , u ( 1 ) = 0. \left\{

\right. ⎩ ⎨ ⎧​−uxx​=f,x∈(0,1),f=exp(x)(−x2−x+1),u(0)=0,u(1)=0.​

对区域$[0,1]$剖分，均匀分成$N$份，即$x_0=0<x_i<x_N=1$，其中$h_x=\frac{1}{N}$，现在要求$u_i=u(x_i),1\le i\le N-1$的取值.利用$\frac{u_{i+1}-2u_i+u_{i-1}}{h_x^2},1\le i\le N-1$可以近似代替$u_{xx}(x_i)$，那么得到一个线性系统：

{ − ( u i + 1 − 2 u i + u i − 1 ) = f i h x 2 , 1 ≤ i ≤ N − 1 , u 0 = 0 , u N = 0 , A u = b . \left\{

\right. ⎩ ⎨ ⎧​−(ui+1​−2ui​+ui−1​)=fi​hx2​,1≤i≤N−1,u0​=0,uN​=0,Au=b.​

这里的$b,u$是一个$N+1$维向量，其中$b_0=0,b_N=0$,其他部分$b_i=f_i{h}_x^2,1\le i\le N-1$，这里的A是一个$(N+1)\times (N+1)$的带状矩阵，满足$A_{0,0}=A_{N,N}=1$，其余对角线元素$A_{i,i}=2$，就是下面这个样子.

A = ( 1 0 0 ⋯ 0 − 1 2 − 1 ⋯ 0 ⋮ ⋮ ⋱ ⋮ 0 ⋯ − 1 2 − 1 0 0 ⋯ 0 1 ) A=\left(

\right) A= ​1−1⋮00​02⋮⋯0​0−1⋱−1⋯​⋯⋯⋮20​00−11​ ​

FD求解泊松方程

Dirichlet边界条件

给定泊松方程如下，其中$\Omega =[0,1{]}^2,\mathbf{x}=(x,y)$：

{ − Δ u = f , x ∈ Ω , u = u 0 , x ∈ ∂ Ω . \left\{

\right. {−Δu=f,x∈Ω,u=u0​,x∈∂Ω.​

引入剖分$M,N,h_x=\frac{1}{M},h_y=\frac{1}{N}$，将区域离散化以后得到网格点$(x_i,y_j),0\le i\le M,0\le j\le N$,记$u_{i,j}=u(x_i,y_j)$，利用中心差分公式???" role="presentation" style="position: relative;">???$\text{???}$得到离散以后的系统，网格剖分如下图所示

{ u x x ( x i , y j ) = u i + 1 , j − 2 u i , j + u i − 1 , j h x 2 + O ( h x 2 ) , 1 ≤ i ≤ M − 1 , u y y ( x i , y j ) = u i , j + 1 − 2 u i , j + u i , j − 1 h y 2 + O ( h y 2 ) , 1 ≤ j ≤ N − 1. \left\{

\right. {uxx​(xi​,yj​)=hx2​ui+1,j​−2ui,j​+ui−1,j​​+O(hx2​),1≤i≤M−1,uyy​(xi​,yj​)=hy2​ui,j+1​−2ui,j​+ui,j−1​​+O(hy2​),1≤j≤N−1.​

{ − u i + 1 , j − 2 u i , j + u i − 1 , j h x 2 − u i , j + 1 − 2 u i , j + u i , j − 1 h y 2 = f i , j , 1 ≤ i ≤ M − 1 , 1 ≤ j ≤ N − 1 , u i , j = u 0 ( x i , y j ) , ( x i , y j ) ∈ ∂ Ω . \left\{

\right. {−hx2​ui+1,j​−2ui,j​+ui−1,j​​−hy2​ui,j+1​−2ui,j​+ui,j−1​​=fi,j​,1≤i≤M−1,1≤j≤N−1,ui,j​=u0​(xi​,yj​),(xi​,yj​)∈∂Ω.​

通过上面的公式可以形成一个线性方程组$Ax=b,A\in R^{(M+1)(N+1)\times (M+1)(N+1)}$，不过由于在实际操作过程中$M,N$很大，会导致$h_x^2,h_y^2$太大，为了避免矩阵$A$病态，我们可以在???" role="presentation" style="position: relative;">???$\text{???}$第一个公式两边同时乘一个$h_x{h}_y$，经过整理以后得到下面这个公式
{ − h y h x u i − 1 , j + ( 2 h x h y + 2 h x h y ) u i , j − h y h x u i + 1 , j − h x h y u i , j − 1 − h x h y u i , j + 1 = f i , j h x h y , u i , j = u 0 ( x i , y j ) , ( x i , y j ) ∈ ∂ Ω . \left\{

\right. {−hx​hy​​ui−1,j​+(2hy​hx​​+2hy​hx​​)ui,j​−hx​hy​​ui+1,j​−hy​hx​​ui,j−1​−hy​hx​​ui,j+1​=fi,j​hx​hy​,ui,j​=u0​(xi​,yj​),(xi​,yj​)∈∂Ω.​

有了上述公式以后，下面我们开始在python中实现$A,b$的构造，这里说明$x=(u_{0,0},u_{0,1},\dots ,u_{0,N},\dots ,u_{M,0},u_{M,1},\dots ,u_{M,N})$.

图片展示了五点差分法构造矩阵的全过程，$u$在边界点的取值已知，所以当循环遇到边界点的索引，此时矩阵的对角元素设定为1.事实上，如果想缩小矩阵的规模，可以把已知量挪到右边，只针对位于网格点内部的函数值$u_{i,j}$求解，形成一个$(M-1)(N-1)\times (M-1)(N-1)$的矩阵，读者可以自己尝试.引入精确解$u=\log (10(x+y{)}^2+(x-y{)}^2+0.5)$，然后形成对应的右端项和边界条件测试的结果如下所示:

代码参考添加链接描述

Neumann边界条件

给定泊松方程如下，其中$\Omega =[0,1{]}^2,\mathbf{x}=(x,y)$：

{ − Δ u + u = f , x ∈ Ω , ∂ u ∂ n = u 0 , x ∈ ∂ Ω . \left\{

\right. {−Δu+u=f,x∈Ω,∂n∂u​=u0​,x∈∂Ω.​

针对Neumann边界条件，泊松方程的设定和上面的Dirichlet???" role="presentation" style="position: relative;">???$\text{???}$不太一样，在控制方程中多了一项关于$u$的线性项，这个设定可以保证该方程有唯一解.如果没有这一项，假设$u^{\ast }$是精确解，那么会发现$u^{\ast }+c$也是精确解，其中$c$是任意常数.在控制方程的处理和上面Dirichlet几乎一样，离散化以后会得到下面这个公式：

{ − h y h x u i − 1 , j + ( 2 h x h y + 2 h x h y + h x h y ) u i , j − h y h x u i + 1 , j − h x h y u i , j − 1 − h x h y u i , j + 1 = f i , j h x h y , ∂ u ∂ n ( x i , y j ) = u 0 ( x i , y j ) , ( x i , y j ) ∈ ∂ Ω . \left\{

\right. {−hx​hy​​ui−1,j​+(2hy​hx​​+2hy​hx​​+hx​hy​)ui,j​−hx​hy​​ui+1,j​−hy​hx​​ui,j−1​−hy​hx​​ui,j+1​=fi,j​hx​hy​,∂n∂u​(xi​,yj​)=u0​(xi​,yj​),(xi​,yj​)∈∂Ω.​

重点是在边界上的处理，方向导数放在四个边界上是有两个$u_x,u_y$，因此这里需要逼近$u_x,u_y$，可以使用一阶差分或者二阶差分和虚拟点来逼近，下面重点介绍虚拟点方法.

这里我们选取边界点$(x_i,y_0),1\le i\le M-1$举例子，注意在4个角需要做特殊处理.则会有

$\frac{\partial u}{\partial n}(x_i,y_0)=\frac{u_{i,1}-u_{i,-1}}{2h_y}=u_0(x_i,y_0),1\le i\le M-1.$

然后代入方程中替换$u_{i,-1}=u_{i,1}-2u_0(x_i,y_0)h_y$即可，对应需要修改右端项$b$的定义.代码参考添加链接描述.

FD求解含时间项的泊松方程

引入方程定义：

{ u ( x , t ) − Δ u ( x , t ) = f ( x , t ) , ( x , t ) ∈ Ω × [ 0 , 1 ] , u ( x , t ) = g ( x , t ) , ( x , t ) ∈ ∂ Ω × [ 0 , 1 ] , u ( x , 0 ) = h ( x ) , x ∈ Ω . \left\{

\right. ⎩ ⎨ ⎧​u(x,t)−Δu(x,t)=f(x,t),(x,t)∈Ω×[0,1],u(x,t)=g(x,t),(x,t)∈∂Ω×[0,1],u(x,0)=h(x),x∈Ω.​

针对方程???" role="presentation" style="position: relative;">???$\text{???}$引入4种差分格式如下，为了方便记忆，定义:\
${\delta }_x{u}_{i,j}^n=\frac{u_{i+1,j}^n-2u_{i,j}^n+u_{i-1,j}^n}{h_x^2}$,
${\delta }_y{u}_{i,j}^n=\frac{u_{i,j+1}^n-2u_{i,j}^n+u_{i,j-1}^n}{h_y^2}$ \
Forward-Euler:

$\frac{u_{i,j}^{n+1}-u_{i,j}^n}{\tau }-{\delta }_x{u}_{i,j}^n-{\delta }_y{u}_{i,j}^n=f_{i,j}^{n+1}.$

Backward-Euler:

$\frac{u_{i,j}^{n+1}-u_{i,j}^n}{\tau }-{\delta }_x{u}_{i,j}^{n+1}-{\delta }_y{u}_{i,j}^{n+1}=f_{i,j}^{n+1}.$

BDF:

$\frac{3u_{i,j}^{n+1}-4u_{i,j}^n+u_{i,j}^{n-1}}{2\tau }-{\delta }_x{u}_{i,j}^{n+1}-{\delta }_y{u}_{i,j}^{n+1}=f_{i,j}^{n+1}.$

CN:

$\frac{u_{i,j}^{n+1}-u_{i,j}^n}{\tau }-\frac{1}{2}({\delta }_x{u}_{i,j}^{n+1}+{\delta }_y{u}_{i,j}^{n+1}+{\delta }_x{u}_{i,j}^n+{\delta }_y{u}_{i,j}^n)=f_{i,j}^{n+1}.$

有了这些格式以后，我们引入$U^n=(u_{0,0}^n,\dots ,u_{0,N}^n,\dots ,u_{M,0}^n,\dots ,u_{M,N}^n)$，则上述4个格式可以理解为4个迭代的线性系统.
Forward-Euler:

$U^{n+1}=A{U}^n+b.$

Backward-Euler:

$A{U}^{n+1}=B{U}^n+b.$

BDF:

$A{U}^{n+1}=B{U}^n+C{U}^{n-1}+b.$

CN:

$A{U}^{n+1}=B{U}^n+b.$

在代码的实现过程中，Forward-Euler不需要求解线性系统，但是这里需要对时间步长$\tau$和空间步长$h_x,h_y$做一个限制来满足稳定性条件，另外三个格式对于任意的时间空间步长都是稳定的，而上面提到的4个格式里面的矩阵需要自己根据离散格式来确定.这里仅仅针对BDF格式做一个具体的说明，其他的类似.

BDF格式

对于时间层，我们每隔$\tau$时间做了一次计算，$U^n$表示当时间层到了第n层的时候那一层的解.将BDF格式整理如下：

( 3 + 4 τ ( 1 h x 2 + 1 h y 2 ) ) u i , j n + 1 − 2 τ h y 2 ( u i , j + 1 n + u i , j − 1 n ) − 2 τ h x 2 ( u i + 1 , j n + u i − 1 , j n ) = 4 u i , j n − u i , j n − 1 + 2 f i , j n + 1 τ .

​(3+4τ(hx2​1​+hy2​1​))ui,jn+1​−2hy2​τ​(ui,j+1n​+ui,j−1n​)−2hx2​τ​(ui+1,jn​+ui−1,jn​)=4ui,jn​−ui,jn−1​+2fi,jn+1​τ.​

具体代码参考下面图片

在BDF代码实现中，计算$U^2$的时候需要$U^1,U^0$，其中$U^0$可以直接通过初始条件获得，$U^1$可以通过CN格式或者Backward-Euler格式获得.
这里我们取$u=(x^2+y^2)\exp (t-t^2)$来做算例测试BDF格式,这里只计算最后一层$u(x,y,1)$的误差.

区域分解算法

上面的FD把整个区域进行剖分，当剖分很细的时候，得到的线性系统维度很大，为了降低存储成本和计算成本，可以将区域划分为若干区域，在每个区域上求解，具体过程将以数学公式来说明.这种算法也称之为cauchy-schwarz算法.

首先需要把区域进行划分，先假设划分为两个区域$\Omega ={\Omega }_0\bigcup {\Omega }_1$，其中${\Omega }_0=[0,0.5]\times [0,1],{\omega }_1=[0.4,1]\times [0,1]$，注意两个区域必须要有重叠overlap.做剖分的时候注意，为了避免麻烦，应该让边界线 { 0.4 } × [ 0 , 1 ] , { 0.5 } × [ 0 , 1 ] \{0.4\}\times [0,1],\{0.5\}\times [0,1] {0.4}×[0,1],{0.5}×[0,1]都位于网格线上,针对泊松方程引入两个子问题.

{ − Δ u 1 n + 1 = f , x ∈ Ω 0 , u 1 n + 1 = u 0 , x ∈ ∂ Ω 0 ∩ ∂ Ω , u 1 n + 1 = u 2 n , x ∈ ∂ Ω 0 ∖ ∂ Ω . \left\{

\right. ⎩ ⎨ ⎧​−Δu1n+1​=f,x∈Ω0​,u1n+1​=u0​,x∈∂Ω0​∩∂Ω,u1n+1​=u2n​,x∈∂Ω0​∖∂Ω.​

{ − Δ u 2 n = f , x ∈ Ω 1 , u 2 n = u 0 , x ∈ ∂ Ω 1 ∩ ∂ Ω , u 2 n = u 1 n + 1 , x ∈ ∂ Ω 1 ∖ ∂ Ω . \left\{

\right. ⎩ ⎨ ⎧​−Δu2n​=f,x∈Ω1​,u2n​=u0​,x∈∂Ω1​∩∂Ω,u2n​=u1n+1​,x∈∂Ω1​∖∂Ω.​

则cauchy-schwarz算法的定义是:

上面有两个子问题，两个线性系统的区别在于右端项，因此可以提前构建好矩阵$A$，然后根据不同的右端项定义两个求解函数solve，具体过程可以参考下列图片：

如果将区域分成3个区域，比如说${\Omega }_0=[0,0.5]\times [0,1],{\omega }_1=[0.4,0.8]\times [0,1],{\omega }_1=[0.7,1]\times [0,1]$，仍然要保证两两之间有overlap，过程与上面算法一样，代码可以参考添加链接描述

1：在上述方法中，最终的问题都被转化成求解一个线性系统$Ax=b$，这里本人调用的是python库numpy自带的线性求解器$linalg.solve(A,b)$，如何设计G-S迭代方法求解，如何证明G-S迭代法对这类线性系统的收敛性.\par
2：线性系统???" role="presentation" style="position: relative;">???$\text{???}$是稀疏的，是否可以设计一种针对这类问题的带状问题进行求解，并且减少矩阵存储量.
3：在介绍FD求解含时泊松方程时，向前欧拉格式不稳定，如何确定合适的$\tau ,h_x,h_y$使得该格式稳定.
4：在区域分解算法中，本人的划分都是针对一个维度，如果对区域做棋盘划分，这个算法该如何写，参考提示添加链接描述

有限元法

FD方法的特点是只能求解离散点的取值，如果想要知道不在网格点的函数值，就需要重新做剖分或者做插值得到，基于此有限元方法finite element(FE)发展起来.FE可以被理解为一种插值方法，还是以二维区域为例，将区域$\Omega$进行剖分得到$M\times N$个小区域，在每个小区域上定义一个分片函数，然后根据网格点的取值进行插值，下面直接以数学公式来进行说明.

首先引入单元函数$\psi (x,y)$，对区域进行剖分以后得到$(x_i,y_j)$，则会有${\phi }_{i,j}=\psi (\frac{x-x_i}{h_x},\frac{y-y_j}{h_y})$，那么这样定义的基函数${\phi }_{i,j}$就满足在第$(i,j)$个网格点上取值为1，在其他网格点上取值为0.
ψ ⁡ ( x , y ) = { ( 1 − x ) ∗ ( 1 − y ) , 0 ≤ x ≤ 1 , 0 ≤ y ≤ 1 , ( 1 + x ) ∗ ( 1 − y ) , − 1 ≤ x ≤ 0 , 0 ≤ y ≤ 1 , ( 1 + x ) ∗ ( 1 + y ) , − 1 ≤ x ≤ 0 , − 1 ≤ y ≤ 0 , ( 1 − x ) ∗ ( 1 + y ) , 0 ≤ x ≤ 1 , − 1 ≤ y ≤ 0. \operatorname{\psi}(x,y)=\left\{

\right. ψ(x,y)=⎩ ⎨ ⎧​(1−x)∗(1−y),0≤x≤1,0≤y≤1,(1+x)∗(1−y),−1≤x≤0,0≤y≤1,(1+x)∗(1+y),−1≤x≤0,−1≤y≤0,(1−x)∗(1+y),0≤x≤1,−1≤y≤0.​
于是我们的有限元近似解定义为，接下来的重点就是确定$u_{i,j}$的取值.有了上面FD的介绍，其实可以先通过FD求出对应的$u_{i,j}$，然后根据???" role="presentation" style="position: relative;">???$\text{???}$求出数值解$u_h$，这个方法很容易实现.不过FE并非如此.
$u_h={\sum }_{i=0}^{i=M}{\sum }_{j=0}^{j=N}{u}_{i,j}{\phi }_{i,j}.$
FE首先需要引入弱形式???" role="presentation" style="position: relative;">???$\text{???}$，其中 v ∈ V h = { v ∣ v ∣ ∂ Ω = 0 } v \in V_h=\{v| v|_{\partial \Omega} = 0\} v∈Vh​={v∣v∣∂Ω​=0}.

{ a ( u , v ) = ∮ Ω ∇ u ⋅ ∇ v d x d y , f ( v ) = ∮ Ω f v d x d y . \left\{

\right. {a(u,v)=∮Ω​∇u⋅∇vdxdy,f(v)=∮Ω​fvdxdy.​

根据方程???" role="presentation" style="position: relative;">???$\text{???}$，我们知道这个公式需要求解$(M+1)(N+1)$个未知数，在边界上的取值已知，还需要求出在区域内部的取值，因此取$v_{i,j}={\phi }_{i,j},1\le i\le M-1,1\le j\le N-1$带入方程???" role="presentation" style="position: relative;">???$\text{???}$就可以得到$(M+1)(N+1)$个方程，这里面有$(M-1)(N-1)$个积分需要计算.但是可以发现由于${\phi }_{i,j}$在大部分区域都是0，因此可以采取扫描单元的做法来简化计算.\
重点：代入$u_h,v={\phi }_{i,j}$得到线性系统得到公式(???)" role="presentation" style="position: relative;">(???)$\text{(???)}$.\
{ a ( u h , φ i , j ) = f ( φ i , j ) , 1 ≤ i ≤ M − 1 , 1 ≤ j ≤ N − 1 , u h ( x i , y j ) = u 0 ( x i , y j ) , ( x i , y j ) ∈ ∂ Ω . \left\{

\right. {a(uh​,φi,j​)=f(φi,j​),1≤i≤M−1,1≤j≤N−1,uh​(xi​,yj​)=u0​(xi​,yj​),(xi​,yj​)∈∂Ω.​
由于$u_h$和${\phi }_{i,j}$的特殊性，区域$\Omega$上的积分只有一小部分非0.
{ ∫ Ω ( ∑ j = 0 j = N u i , j ∇ φ i , j ) ∇ φ i ^ , j ^ d Ω = ∫ Ω f φ i ^ , j ^ d Ω , ∑ m = 0 M N − 1 ∫ Ω m ( ∑ i = 0 i = M ∑ j = 0 j = N u i , j ∇ φ i , j ) ∇ φ i ^ , j ^ d Ω = ∑ m = 0 M N − 1 ∫ Ω m f φ i ^ , j ^ d Ω . \left\{

∫Ω(∑j=0j=Nui,j∇φi,j)∇φi^,j^dΩ=∫Ωfφi^,j^dΩ,∑m=0MN−1∫Ωm(∑i=0i=M∑j=0j=Nui,j∇φi,j)∇φi^,j^dΩ=∑m=0MN−1∫Ωmfφi^,j^dΩ." role="presentation">∫Ω(∑j=Nj=0ui,j∇φi,j)∇φi^,j^dΩ=∫Ωfφi^,j^dΩ,∑MN−1m=0∫Ωm(∑i=Mi=0∑j=Nj=0ui,j∇φi,j)∇φi^,j^dΩ=∑MN−1m=0∫Ωmfφi^,j^dΩ.$$\begin{array}{l}{\int }_{\mathrm{\Omega }}(\sum _{j=0}^{j=N}{u}_{i,j}\mathrm{\nabla }{\phi }_{i,j})\mathrm{\nabla }{\phi }_{\hat{i},\hat{j}}d\mathrm{\Omega }={\int }_{\mathrm{\Omega }}f{\phi }_{\hat{i},\hat{j}}d\mathrm{\Omega },\\ \sum _{m=0}^{MN-1}{\int }_{{\mathrm{\Omega }}_m}(\sum _{i=0}^{i=M}\sum _{j=0}^{j=N}{u}_{i,j}\mathrm{\nabla }{\phi }_{i,j})\mathrm{\nabla }{\phi }_{\hat{i},\hat{j}}d\mathrm{\Omega }=\sum _{m=0}^{MN-1}{\int }_{{\mathrm{\Omega }}_m}f{\phi }_{\hat{i},\hat{j}}d\mathrm{\Omega }.\end{array}$$

\right. {∫Ω​(∑j=0j=N​ui,j​∇φi,j​)∇φi^,j^​​dΩ=∫Ω​fφi^,j^​​dΩ,∑m=0MN−1​∫Ωm​​(∑i=0i=M​∑j=0j=N​ui,j​∇φi,j​)∇φi^,j^​​dΩ=∑m=0MN−1​∫Ωm​​fφi^,j^​​dΩ.​
当$(x_i,y_j),(x_{\hat{i}},y_{\hat{j}})$同时位于${\Omega }_m$的顶点，才有${\int }_{{\Omega }_m}\nabla {\phi }_{i,j}\nabla {\phi }_{\hat{i},\hat{j}}\ne 0$.
${\int }_{{\Omega }_m}({\sum }_{i=0}^{i=M}{\sum }_{j=0}^{j=N}{u}_{i,j}\nabla {\phi }_{i,j})\nabla {\phi }_{\hat{i},\hat{j}}d\Omega ={\int }_{{\Omega }_m}{\sum }_{|i-\hat{i}|\le 1,|j-\hat{j}|\le 1}\nabla {\phi }_{i,j}\nabla {\phi }_{\hat{i},\hat{j}}d\Omega .$
因此才有下面的公式成立.
$a_{i,j}={\sum }_{Nodes[i],Nodes[j]\in T_m}{\oint }_{T_m}\nabla {\varphi }_i\nabla {\varphi }_{j}dxdy+{\sum }_{Nodes[i],Nodes[j]\in \partial {\Omega }_m}{\oint }_{\partial {\Omega }_m}\beta {\varphi }_i{\varphi }_{j}ds$
其中Nodes$[i] $表示在整体排序中的第$i$个网格点,$T_m$表示第m个单元，区域一共有$M\times N$个单元.
类似地，我们可以给出右端项，也就是所谓的载荷向量的b的定义
$b_i={\sum }_{Nodes[i]\in T_m}{\oint }_{T_m}f{\varphi }_{i}dxdy+{\sum }_{self.Nodes[i]\in \partial {\Omega }_m}{\oint }_{\partial {\Omega }_m}\beta g{\varphi }_{i}ds$

在给出这些矩阵元素和右端项元素的定义计算以后，我们就必须要解决积分问题，我们需要对区域进行离散化处理，选取有限个点作为积分节点，根据数值逼近的原理，我们知道在区域$[0,1{]}^2$中，我们只要选取以下4个点$A_1=(\frac{1-\sqrt{3}/3}{2},\frac{1-\sqrt{3}/3}{2}),A_2=(\frac{1+\sqrt{3}/3}{2},\frac{1-\sqrt{3}/3}{2}),A_3=(\frac{1+\sqrt{3}/3}{2},\frac{1+\sqrt{3}/3}{2}),A_4=(\frac{1-\sqrt{3}/3}{2},\frac{1+\sqrt{3}/3}{2})$作为积分点，
${\oint }_{[0,1{]}^2}fdxdy=0.25\ast {\sum }_{i=1}^4(f(A_i))$
那么上述积分近似就有2阶代数精度，即对于$1,x,y,x^2,y^2,xy$这些函数的线性组合，上述积分等式恒成立。本人代码里面FENET里面定义了$self.gp\_pos=[\frac{1-\sqrt{3}/3}{2},\frac{1-\sqrt{3}/3}{2}]$以便于在每个积分单元采取这样4个高斯积分节点。在这些定义完以后，本人定义了$Int\_basic\_basic$函数表示在某积分单元上两个基函数之间梯度的积分，定义了$Int\_F\_basic$表示在某积分单元上f和某基函数的积分，定义$Bou\_basic\_basic$表示在某Neumann边界线段单元上两个基函数的线积分，定义$Bou\_Neu\_basic$表示在Neumann边界上函数g和基函数在该线段单元上的线积分，代码参考添加链接描述

深度学习求解PDE

传统方法FD,FE最大的问题在于当坐标是高维向量的时候，无法对高维空间进行剖分，这个会造成维度灾难curse of dimensionality.而这些年深度学习的兴起给我们提供了一个新的方法来求解PDE，下图???" role="presentation" style="position: relative;">???$\text{???}$就是一个典型的全连接网络，该网络接受4个维度，输出一个标量，可以这么理解$FNN(\mathbf{x})$定义了一个从$R^4$到$R^1$的函数，这个函数是通过激活函数和线性组合得到的，使用深度学习求解PDE的核心就是选择适当的参数$W,b$使得$FNN(\mathbf{x})$尽可能逼近精确解.

一个通用的全连接网络NN表达式如下???" role="presentation" style="position: relative;">???$\text{???}$所示，在图片???" role="presentation" style="position: relative;">???$\text{???}$中有两个隐藏层，因此$n=2$，网络搭建的细节可以参考文献\cite{2nd}.
{ f i ( s ) = ϕ ( W i , 2 ⋅ ϕ ( W i , 1 s + b i , 1 ) + b i , 2 ) , z = f n ∘ f n − 1 ∘ … ∘ f 0 ( x ) , y ( x ) = W z + b . \left\{

\right. ⎩ ⎨ ⎧​fi​(s)=ϕ(Wi,2​⋅ϕ(Wi,1​s+bi,1​)+bi,2​),z=fn​∘fn−1​∘…∘f0​(x),y(x)=Wz+b.​
正如上面提到的，利用深度学习求解PDE的核心在于选择合适的权重和偏置，我们以$\theta$来表示神经网络NN中的权重和偏置，这些也被称为NN的参数，我们需要根据PDE来选择等价的优化问题，然后利用优化中的迭代算法来迭代更新参数$\theta$，下面我们引入几个基本的算法来做介绍，测试函数都是$u=\log (10(x+y{)}^2+(x-y{)}^2+0.5)$.

Deep Ritz

对于泊松方程，利用极小位能原理可以把PDE转化为等价的变分问题
{ min ⁡ u ∫ Ω ∣ ∇ u ∣ 2 − 2 f u d Ω ,  subject to  u = u 0 , x ∈ ∂ Ω . \left\{

\right. ⎩ ⎨ ⎧​​umin​∫Ω​∣∇u∣2−2fudΩ, subject to u=u0​,x∈∂Ω.​
我们需要根据变分问题???" role="presentation" style="position: relative;">???$\text{???}$来更新参数，首先就需要把积分离散化，我们将在区域$\Omega$采样$n_I$个点 { x i } i = 0 n I − 1 \{\mathbf{x}_i\}_{i=0}^{n_I - 1} {xi​}i=0nI​−1​作为积分点.下面我们将详细介绍利用python如何实现deep ritz方法：

第一：我们需要对区域$\Omega$进行离散化，需要在内部采集样本点 X = { x i } i = 0 n I − 1 X = \{\mathbf{x}_i\}_{i=0}^{n_I - 1} X={xi​}i=0nI​−1​，从上面公式中可以发现还需要准备好函数$f(X)$数据集，采样参考以下图片.

第二：我们需要在区域边界$\partial \Omega$进行采样，得到边界上的样本点 X = { x i } i = 0 n B − 1 X=\{\mathbf{x}_i\}_{i=0}^{n_B - 1} X={xi​}i=0nB​−1​，以及$u_0(X)$.

第三：我们需要搭建全连接网络，pytorch自带的torch.nn.Linear来作为每一层的初始参数，搭建好一个输入节点为2，输出节点为1的Net.

第四：我们需要根据网络和区域内部边界样本点和???" role="presentation" style="position: relative;">???$\text{???}$定义一个损失函数，对于这种约束优化问题，我们可以采取惩罚函数方法来定义一个损失函数???" role="presentation" style="position: relative;">???$\text{???}$，其中$u(\theta )$表示FNN的预测解，$\nabla u(\theta )$通过调用pytorch中的自动求导函数torch.autograd.grad函数实现，损失函数以及一部分的训练过程参考以下图片.
{ L = L I + β L B , L I = 1 n I ∑ i = 0 n I − 1 ∣ ∇ u ( θ ) i ∣ 2 − 2 f i u ( θ ) i , L B = 1 n B ∑ i = 0 n B − 1 ( u ( θ ) − u 0 ) i 2 . \left\{

\right. ⎩ ⎨ ⎧​​L=LI​+βLB​,LI​=nI​1​i=0∑nI​−1​∣∇u(θ)i​∣2−2fi​u(θ)i​,LB​=nB​1​i=0∑nB​−1​(u(θ)−u0​)i2​.​

第五：有了损失函数以后，我们需要对网络Net的参数做迭代优化，比如最简单的梯度下降法：
${\theta }^{k+1}={\theta }^k-c\nabla L({\theta }^k)$
对于大部分人来说，手动求出$\nabla L({\theta }^k)$极为困难，幸好pytorch自带这个功能，只需要调用L.backwad()函数就可以得到，并且调用pytorch自带的优化器比如说torch.optim.Adam(Net.parameters(),lr)就可以自动实现参数的更新,lr是步长，也称learning rate，参考论文，代码参考添加链接描述.

考虑该变分问题，可以发现这种方法有几个特点，首先这个问题只需要利用一阶梯度信息，其次这个问题不需要额外增加测试函数，正是这两点，使得Ritz方法非常适合求解高维PDE。但是另一方面，如果使用惩罚函数方法求解该优化问题，损失函数的最优值不见得是0，此时优化过程难以控制。因此如果可以寻找到一种方法，使得近似解自动满足边界条件，就可以避免惩罚项的出现。

Deep Nitsche

为了进一步提高精度，有的论文修改了针对泊松方程???" role="presentation" style="position: relative;">???$\text{???}$的变分问题，得到了一个Deep Nitsche Method，和Deep Ritz唯一的区别在于变分问题不一样，其他的区域采样或者训练如出一辙，细节参考论文，代码参考添加链接描述
Galerkin
Galerkin的特点，只需要利用pytorch对模型求一阶梯度（可以大幅度降低计算时间和存储成本，相比于PINN而言，PINN需要额外存储几个梯度数组），但是必须要提前准备好测试函数数据，这个对数学功底要求很高。以及，Galerkin仍然不能用来求解高维PDE，原因在于高维空间中的测试函数 $v$ 是定义在高维超立方体的函数，这个定义一定需要利用for循环定义，而且高维超立方体的采样不容易。
Galerkin原理是从另一个角度来对泊松方程进行降阶处理，具体参考公式
{ ∫ Ω ∇ u ⋅ ∇ v − f v d Ω = 0 ,  subject to  u = u 0 , x ∈ ∂ Ω , v ∈ V h = { v ∣ v ∣ ∂ Ω = 0 } . \left\{

\right. ⎩ ⎨ ⎧​​∫Ω​∇u⋅∇v−fvdΩ=0, subject to u=u0​,x∈∂Ω,v∈Vh​={v∣v∣∂Ω​=0}.​
在代码的实现过程中，除了采集样本点$X$和准备右端项$f(X)$以外，这里还需要额外准备好$v(X),\nabla v(X)$，在本人的代码实现中，$v$取得是前面有限元提到的单元函数${\phi }_{i,j}$，值得注意的是，此时我们需要选取尽可能多个测试函数$v$带入方程???" role="presentation" style="position: relative;">???$\text{???}$才能得到好的效果，针对测试函数的选取，又可以有多种策略.
特别注意，Galerkin对数学功底要求很大的另一个原因在于推导弱形式，上面这个二维泊松方程的弱形式只需要利用二维积分的格林公式即可，但是我特意没用含时间的热方程举例子，原来就在于即使是那种简单例子，弱形式的推导就已经有难度了。下面罗列一个NS方程的例子.

上面这个NS方程的弱形式就连很多数学系本科的同学也很难推导（此时反而利用PINN会更加方便）。PINN和Galerkin都有一个特点，就是损失函数直接描述了PDE的求解情况，根据损失函数的定义我们知道，神经网络的训练就是尽可能让损失函数接近于0，和0越接近，说明模型训练的效果越好。

PFNN

在论文中，以区域内点为中心构造了${\phi }_{i,j},1\le M-1,1\le j\le N-1$的$(M-1)(N-1)$个测试函数，然后代入方程???" role="presentation" style="position: relative;">???$\text{???}$来构造损失函数，同时使用length factor和另一个网络将区域内部和边界分开处理取得了不错的效果，损失函数参考???" role="presentation" style="position: relative;">???$\text{???}$，其中由于${\phi }_{i,j}$的特殊性质，$L_{i,j}$的积分计算可以限制在以$(x_i,y_j)$为中心的周围4个单元上，这个可以简化积分的计算.
{ L = ∑ i = 1 M − 1 ∑ j = 1 N − 1 L i , j + β L B , L i = ( ∫ Ω ∇ u ⋅ ∇ φ i , j − f φ i , j d Ω ) 2 , L B = 1 n B ∑ i = 0 n B − 1 ( u ( θ ) − u 0 ) i 2 . \left\{

\right. ⎩ ⎨ ⎧​​L=i=1∑M−1​j=1∑N−1​Li,j​+βLB​,Li​=(∫Ω​∇u⋅∇φi,j​−fφi,j​dΩ)2,LB​=nB​1​i=0∑nB​−1​(u(θ)−u0​)i2​.​
代码参考添加链接描述
\subsubsection{WAN}
对抗神经网络求解泊松方程的核心是利用两个网络，$N_u$模拟精确解$u$，$N_v$模拟测试函数$v$，在训练过程中需要不断调整网络$N_v$的参数使得$v$尽量填充整个$V_h$空间.假设$\theta$表示网络$N_u$的参数，$\eta$表示网络$N_v$的参数，那么WAN的核心思想可以理解为求解下列优化问题???" role="presentation" style="position: relative;">???$\text{???}$，其中$L$的定义类似于:
${\max }_{\eta }{\min }_{\theta }L.$
这里需要注意，为了使得$v\in V_h$，我们需要在网络$N_v$的基础上乘一个length factor，与此同时这两个网络需要交替训练，具体细节参考论文和代码添加链接描述

PINN

physics informed neural network最早出现于论文，核心在于下面这个二阶损失函数???" role="presentation" style="position: relative;">???$\text{???}$，里面的梯度也是调用pytorch的自动求导功能完成的，这个损失函数非常直观，事实上，根据本人的数值实验，在求解二维泊松方程这类问题上，PINN的效果远比上面提到的方法要好，PINN的代码实现和上面的如出一辙，这里不赘述，代码参考添加链接描述.
{ L = L I + β L B , L I = 1 n I ∑ i = 0 n I − 1 ( − Δ u − f ) i 2 , L B = 1 n B ∑ i = 0 n B − 1 ( u ( θ ) − u 0 ) i 2 . \left\{

\right. ⎩ ⎨ ⎧​​L=LI​+βLB​,LI​=nI​1​i=0∑nI​−1​(−Δu−f)i2​,LB​=nB​1​i=0∑nB​−1​(u(θ)−u0​)i2​.​

PINN的优势在于损失函数的形成非常直观，如果使用高阶优化器可以得到非常好的效果。PINN的缺点在于不能求解高维PDE，原来在于自动求导的时候需要求二阶导数，想象一下，如果问题是一个100维的，求解一阶导数的时候只需要利用下面这一行代码就可以解决。

u_x, = torch.autograd.grad(u, inset.X, create_graph=True, retain_graph=True,
                                   grad_outputs=torch.ones(inset.size,1))
1
2

但是，如果求解二阶导数，就必须要使用for循环，在形成损失函数的时候出现这个for循环将会导致优化过程特别缓慢，没有特别意义。

u_xx = torch.zeros(100,N)
for i in range(100):
u_xx[i:i + 1,:], = torch.autograd.grad(u_x[i:i + 1,:], inset.X, create_graph=True, retain_graph=True,
                                   grad_outputs=torch.ones(inset.size,1))
1
2
3
4

parametric PDE

仍然考虑泊松方程，不过这里的右端项$f=f(\mathbf{x},\mu ),\mu \in \mathcal{P}$，现在需要求解对于不同参数$\mu$的精确解，我们以
$u(x,y,\mu )=\exp (-\frac{(x-{\mu }_1{)}^2+(y-{\mu }_2{)}^2}{2{\mu }_3^2})$
为例子，其中$\mu =({\mu }_1,{\mu }_2,{\mu }_3)\in [0,1]\times [0,1]\times [0.1,1]$.对于这个问题，其实做法和上面的泊松方程处理类似，可以把$u$理解为在一个5维空间的函数，也就是说在求解过程中，只需要在采样过程将原来的二维空间采样修改成现在的5维空间采样，把空间坐标和参数都当成输入传进网络就可以得到一个含参问题的解.引入网络Net接收$(x,y,{\mu }_1,{\mu }_2,{\mu }_3)$为输入，输出对应的解即可，网络搭建参考???" role="presentation" style="position: relative;">???$\text{???}$，这里为了简单放了一个二维参数空间的FNN示意图，代码参考添加链接描述，求解结果如下，参考论文

代码说明

为了更好理解，这里我们将针对一个典型的PINN求解含参PDE问题代码做一个细致的说明，代码参考添加链接描述

第一步：导入相关的库，torch是用于神经网络搭建和训练的库，安装参考官方说明添加链接描述，matplotlib用于可视化，用法比较复杂，本人也需要百度才能知道如何针对需求使用命令，time用于计时，numpy是python最主要的数学库，注意最后那个bfgs是本人课题组自己编写的一个优化器，读者可以注释掉，最后优化的时候用torch自带的LBFGS优化器代替.

第二步：定义代码的运行设备device，这里本人选择GPU也就是cuda，读者需要自己确定自己电脑或者服务器是否有GPU，其次定义精确解UU和右端项FF，UU里面的order表示精确解的阶数，比如说$[0,1]$表示$u_y$，bounds表示求解区域和参数空间.

第三步：定义INSET，BDSET,TESET这三个类，分别表示内点集合，边界点集合，测试集集合，其中还会封装一些跟样本点有关的数据进去，比如说右端项.

第四步：搭建网络，这个Net类接受两个参数，layer表示网络的层数和每层的神经元个数，dtype表示网络的数据类型，我后面取的是double，根据我的经验，数据精度高得到的效果也会更好一些，当然了，这个也更消耗内存.

第五步：修改数据类型，上面的INSET那三个类里面的数据默认都是float，为了和Net数据类型匹配，这里需要把INSET那几个类里面的数据类型进行修改，比如说$x=x.type(dtype)$就可以把$x$从float变成dtype一样的类型，然后为了后面训练的时候加速，我们一开始就提出要把代码送到GPU上训练，但是数据一开始默认在CPU上，因此我们还需要写一个函数把所有的数据搬运到GPU，为了训练完以后数据可以回到CPU，这里相应地也要定义一个函数把数据从GPU搬运回CPU.

第六步：定义损失函数和训练函数，损失函数之前讲过，这里不赘述，主要介绍一下训练函数的定义.训练需要优化器optim，我一般都取LBFGS或者我们课题组自己写的BFGS，如果使用Adam优化器，训练过程的核心就是里面注释掉的代码，但是BFGS或者LBFGS的核心就是closure这个函数.这里做一些解释，optim.zero$\_$grad表示清除优化器里面存储的梯度，防止每一次的梯度进行叠加，loss.backward表示反向传播，这里会自动计算损失函数关于网络参数$W,b$的梯度，optim.step是负责根据优化器来对参数$W,b$做更新，三者缺一不可，等到这三个结束以后还需要重新计算一下loss，这里本人不是特别理解，但是本人的实验指出这一步很重要，如果注释掉，结果不对.

最后就是对上面定义的函数，网络，类进行实例化和调用，然后一切就结束了.

总结

上述提到的几种方法各有特点，这里对这些方法做一个总结.

神经网络方法求解PDE的困难和弊端
1：神经网络方法求解PDE缺乏收敛性保证。通过上面的介绍，我们可以这么理解，神经网络求解PDE是把原始问题转化为一个优化问题，然后不断通过更新神经网络的参数来达到min loss function的做法，比如说PINN和Galerkin思想，神经网络形成的损失函数只要不断接近0，那么说明神经网络学习的越好。但是现在有一个问题，会不会出现一种情况，形成了一个损失函数，不管怎么更新网络的参数，损失函数都维持在一个较高的水平呢？或者另一种情况，神经网络需要花大量时间训练才能把损失函数降到一个较低的水平。下面我们看一个障碍流问题的例子。

上面这个就是计算区域，这个问题描述了流体从入口进入管道，管道内部有一个障碍，然后现在想预测流体在出口的流速。这个问题可以使用PINN的思想，先采集样本点，然后搭建网络，优化损失函数，下面先给出这个问题的一个参考解，这个解是通过OpenFoam求出来的。

本人训练过程就发现，这个问题的损失函数特别难下降，具体可以参考下面这个图片，虽然原理在那摆着，但是这个问题却难以求解。

本人总结了一句话：PINN或者是Galerkin这类神经网络方法，都只能求解简单的PDE，对于简单的PDE，往往可以得到优于传统数值方法的效果，但是训练时间绝对更长。对于复杂的PDE，比如说NS方程，尤其是带有复杂边界的问题，这类简单的NN方法无法得到多好的效果，但是另一方面，传统数值方法求解复杂的PDE也同样需要消耗大量的计算精力。

个人以为NN求解PDE的优势应该在于求解高维的简单的PDE，比如说高维的Dirichlet边界泊松方程，这种方程由于计算区域维度更高，所以传统数值方法很容易发生维度灾难，但是神经网络方法却可以通过蒙特卡洛等方法对计算区域进行采样，轻而易举求解。

NN求解PDE目前更多还是应用于科研领域，对于实际应用目前难以推广，一个非常大的原因就是NN求解PDE没有收敛性的理论保证，对于PDE问题，NN就像一个炼丹炉，谁也不知道这个问题到底会求解成什么样子。反而是传统数值方法（FD,FE），虽然求解PDE很消耗时间，计算复杂度也高，但是这类方法的理论收敛经过几十年的研究论证，已经可以做到在求解问题之前就推断出计算结果大概可以到什么程度，这一点对于应用非常重要。