GSoC 2017：PDE 算子的高效离散化

2017 年 9 月 6 日 | Shivin Srivastava, Christopher Rackauckas

本项目旨在为 JuliaDiffEq 使用 有限差分法(FDM) 方法构建一个 PDE 求解器。我们选择 FDM 方法而不是 FEM 和 FVM，因为 FEM 和 FVM 已经有很多现成的工具箱，而 FDM 则不然。此外，在很多情况下，问题的几何形状非常简单，可以使用 FDM 方法进行求解，这些方法由于能够避免矩阵乘法的瓶颈步骤而快得多，而是使用线性变换来模拟矩阵乘法的效果。由于矩阵乘法本质上将一个向量元素转换为相邻元素的加权和，因此这可以很容易地通过一个特殊函数来实现，该函数以最优的 O(n)\mathcal{O}(n) 时间对向量进行操作。

结果是名为 DiffEqOperators.jl 的新软件包，它创建了偏微分算子的高效离散化，从而将 PDE 转换为 ODE，可以通过现有的 ODE 求解器进行高效地求解。DerivativeOperator 基于中心差分方案，用于近似某个点的导数，而 UpwindOperators 基于单边差分方案，其中解通常是一个在特定方向上传播的波。该软件包还支持各种边界条件，如 Dirichlet、Neumann、周期 和 Robin 边界条件。

  1. 动机
  2. 使用 DiffEqOperators 求解热方程
  3. 中心导数在何处失效？
  4. 未来工作
  5. 致谢

动机

有限差分法的基本思想是生成与微分算子相对应的有限差分权重，以允许一定程度的近似。时间和空间变量被划分成一个网格，其中 h=Δx=1N+1h = \Delta x = \frac{1}{N+1} 和 xi=ihx_i = ih 对于 i=0,1,...,N+1i = 0, 1,...,N+1 和 k=Δt=TM+1k = \Delta t = \frac{T}{M+1} 和 tj=jkt_j = jk 对于 j=0,1,...,M+1j = 0, 1,...,M+1.
离散未知数是标量 uiju_i^j 对于上面给出的i和j的值，希望 uiju_i^j 是 u(xi,tj)u(x_i, t_j) 的近似值。该方程的右侧通过设置 fij=f(xi,tj)f_i^j = f(x_i, t_j) 进行离散化。我们也使用符号

Uj=(u1ju2j⋮uNj) U^{j} = \begin{pmatrix} u_1^j\\ u_2^j \\ \vdots \\ u_N^j \end{pmatrix}

来表示时间 tjt_j 在空间网格上的近似值向量。

∂u∂t(xi,tj)≈u(xi,tj+1)−u(xi,tj)k\frac{\partial u}{\partial t}(x_i,t_j) \approx \frac{u(x_i,t_{j+1}) - u(x_i,t_j)}{k}

对于二阶空间导数，通过结合前向差分商和后向差分商，我们得到了中心逼近

∂2u∂x2(xi,tj)≈u(xi+1,tj)−u(xi,tj)h−u(xi,tj)−u(xi−1,tj)hh=u(xi−1,tj)−2u(xi,tj)+u(xi+1,tj)h2 \frac{\partial^2u }{\partial x^2}(x_i,t_j) \approx \frac{\frac{u(x_{i+1},t_j) - u(x_i,t_j)}{h} - \frac{u(x_i,t_j) - u(x_{i-1},t_j)}{h}}{h} = \frac{u(x_{i-1},t_j) - 2u(x_i,t_j) + u(x_{i+1}, t_j)} {h^2}

因此，对应二阶偏导数的权重为 [1,−2,1][1, -2, 1]。有限差分法通过用离散未知数代替网格点上解的精确值来模仿这些逼近。

在这个特定情况下，我们最终得到了以下方案

{ujj+1−uijk−ui−1j−2uij+ui+1jh2=fij for i=1,...,N,j=1,...,Mui0=u0(xi) for i=1,...,Nu0j=uN+1j=0 for j=1,...,M+1 \left\{\begin{matrix} \frac{u_j^{j+1} - u_i^j}{k} - \frac{u_{i-1}^j - 2u_i^j + u_{i+1}^j}{h^2} = f_i^j \quad\textit{ for }\quad i = 1,...,N,j = 1,...,M \\ u_i^0 = u_0(x_i) \quad\textit{ for }\quad i = 1,...,N\\ u_0^j = u_{N+1}^j = 0 \quad\textit{ for }\quad j = 1,...,M+1\\ \end{matrix}\right.

这是将模板重写为递归式。用向量形式写出来，我们得到：- Uj+1−Ujk+AhUj=Fj for j=1,...,M\frac{U^{j+1} - U^j}{k} + A_h U^j = F^j \quad\textit{ for }\quad j = 1,...,M

当我们要将此运算符应用于向量 UU 时，权重向量就变成了一个称为变换矩阵的矩阵 AhA_h

Ah=(2−10⋯0−12−1⋯0⋮⋱⋱⋱⋮0⋯−12−10⋯0−12)such that ∂2Uj∂x2≈Ah∗Uj A_h = \begin{pmatrix} 2 & -1 & 0 & \cdots & 0\\ -1 & 2 & -1 & \cdots & 0\\ \vdots & \ddots & \ddots & \ddots & \vdots \\ 0 & \cdots & -1 & 2 & -1 \\ 0 & \cdots & 0 & -1 & 2 \\ \end{pmatrix} \quad\textit{such that }\quad \frac{\partial^2U^j}{\partial x^2} \approx A_h*U^j

但矩阵乘法成本很高，因此最好使用二阶偏导数的线性算子，而不是变换矩阵。 它看起来像这样：-

function double_partial(x,dx) for i in 2:length(dx)-1 dx[i] = -1*x[i-1] + 2*x[i] + -1*x[i+1] end dx[1] = 2*x[1] + -1*x[2] dx[end] = -1*x[end-1] + 2*x[end] end

此函数对向量进行操作，具有最佳的 O(n)\mathcal{O}(n) 时间复杂度，与矩阵乘法的低效 O(n2)\mathcal{O}(n^2) 时间复杂度相比，同时避免了稀疏矩阵的开销。

为了将偏微分方程转换为常微分方程，我们在空间上进行离散化，但在时间上不进行离散化。然后，这个常微分方程可以通过现有的求解器有效地求解。我们的半线性热方程，也称为反应扩散方程，转换为以下常微分方程：ui′=Ahui+f(t,ui)u_i' = A_{h}u_i + f(t,u_i) 其中 AA 是线性算子，而不是变换矩阵。因此，我们需要使 **DifferentialEquations.jl** 的常微分方程求解器也兼容线性算子。

由于手动计算泰勒系数很繁琐，Fornberg 给出了一种 算法 ，可以高效地计算任意导数和逼近阶的泰勒系数。这些模板可以通过对相邻点的适当加权求和来有效地计算任意点的导数。例如，[−1,2,−1][-1, 2, -1] 是用于计算某一点的二阶导数的二阶模板。

在 **DiffEqOperators.jl** 中，我们可以轻松地从算子中提取任意导数和逼近阶的模板。例如。

# Define A as a DerivativeOperator of 4th order and of 2nd order of accuracy

julia> A = DerivativeOperator{Float64}(4,2,1.0,10,:Dirichlet0,:Dirichlet0)
julia> A.stencil_coefs
7-element SVector{7,Float64}:
  -0.166667
   2.0
  -6.5
   9.33333
  -6.5
   2.0
  -0.166667

如果我们想将算子作为矩阵乘法（稀疏或稠密）应用，我们可以通过提取线性算子的变换矩阵来轻松实现，它看起来像这样：-

julia> full(A)
10×10 Array{Float64,2}:
 9.33333   -6.5        2.0       …   0.0        0.0        0.0
-6.5        9.33333   -6.5           0.0        0.0        0.0
 2.0       -6.5        9.33333       0.0        0.0        0.0
-0.166667   2.0       -6.5           0.0        0.0        0.0
 0.0       -0.166667   2.0          -0.166667   0.0        0.0
 0.0        0.0       -0.166667  …   2.0       -0.166667   0.0
 0.0        0.0        0.0          -6.5        2.0       -0.166667
 0.0        0.0        0.0           9.33333   -6.5        2.0
 0.0        0.0        0.0          -6.5        9.33333   -6.5
 0.0        0.0        0.0           2.0       -6.5        9.33333

julia> sparse(A)
 10×10 SparseMatrixCSC{Float64,Int64} with 58 stored entries:
 [1 ,  1]  =  9.33333
 [2 ,  1]  =  -6.5
 [3 ,  1]  =  2.0
 [4 ,  1]  =  -0.166667
 [1 ,  2]  =  -6.5
 [2 ,  2]  =  9.33333
 [3 ,  2]  =  -6.5
 ⋮
 [7 ,  9]  =  2.0
 [8 ,  9]  =  -6.5
 [9 ,  9]  =  9.33333
 [10,  9]  =  -6.5
 [7 , 10]  =  -0.166667
 [8 , 10]  =  2.0
 [9 , 10]  =  -6.5
 [10, 10]  =  9.33333

模板乘法是极度并行的，**DiffEqOperators.jl** 已经考虑到了这一点。

使用 DiffEqOperators 求解热方程

现在让我们使用二维 space x time 网格上的显式离散化来求解著名的热方程。热方程为：-

∂u∂t−∂2u∂x2=0\frac{\partial u}{\partial t} - \frac{\partial^2 u}{\partial x^2} = 0

对于这个例子，我们考虑一个狄利克雷边界条件，初始分布为抛物线。由于我们在边界上固定了值（在本例中相等），经过很长时间后，我们预计一维杆将以线性方式被加热。

julia> using DiffEqOperators, DifferentialEquations, Plots
julia> x = -pi : 2pi/511 : pi;
julia> u0 = -(x - 0.5).^2 + 1/12;
julia> A = DerivativeOperator{Float64}(2,2,2pi/511,512,:Dirichlet,:Dirichlet;BC=(u0[1],u0[end]));

这是设置问题的代码。首先，我们定义域，它只是一条直线，被分成 512 段。然后，我们定义初始条件，它是一个关于 x 坐标的抛物线函数。

最后，我们初始化二阶导数和二阶逼近阶的 DerivativeOperator。我们告诉网格步长值、域的总长度和两端的边界条件。请注意，由于我们在这里应用狄利克雷边界条件，我们需要告诉边界上的值，它以元组的形式作为最后一个参数给出。

现在将方程作为常微分方程求解，我们有：-

julia> prob1 = ODEProblem(A, u0, (0.,10.));
julia> sol1 = solve(prob1, dense=false, tstops=0:0.01:10);
# try to plot the solution at different time points using
julia> plot(x, [sol1(i) for i in 0:1:10])

Heat Equation

请注意热分布如何随着时间的推移“变平”，正如预期的那样，最终趋于从左端到右端线性增加。

中心导数在何处失效？

并非所有偏微分方程都可以用中心差分法求解，例如KdV 波方程。经过常微分方程求解器的几次迭代后，波开始分裂，即它很快变得不稳定。

KdV using central derivatives

对于这些非常特殊的情况，已经设计出了 迎风格式。它表示一类用于求解双曲型偏微分方程的数值离散化方法。它们试图通过使用根据特征速度符号确定的方向进行偏差的差分来离散化双曲型偏微分方程。例如一维线性对流方程

∂u∂t+a∂u∂x=0 \frac{\partial u}{\partial t} + a\frac{\partial u}{\partial x}=0

描述了一个沿 x 轴以速度 aa 传播的波。如果 aa 为正，则上述方程的传播波解向右传播，左侧称为迎风侧，右侧称为背风侧。如果空间导数 ∂u∂x\frac{\partial u}{\partial x} 的有限差分格式包含迎风侧的更多点，则该格式称为 **迎风格式**。考虑二阶迎风格式的情况，定义

a+=max(a,0) a^{+} = max(a,0) a−=min(a,0)a^{-} = min(a,0) ux−=3uin−4ui−1n+ui−2n2Δxu_x^- = \frac{3u_i^n-4u_{i-1}^n+u_{i-2}^n}{2\Delta x}