基于通量限制器的时间分数阶Burgers方程数值解法

Numerical Solution of Time-Fractional Burgers Equation Based on Flux Limiter

摘要: 对于求解可能具有激波等不连续点的时间分数阶守恒律，当分数阶γ接近于0时，目前还没有有效的方法求解时间分数阶非线性离散系统，本文以时间分数阶Burgers方程为例，运用多重网格迭代方法进行求解，对于对流项，采用通量限制器，使得新的数值通量在光滑区域变为高阶通量而在间断附近变为低阶通量，从而使问题的解达到更高阶精度，并在不同的γ取值以及不同的初边值条件下进行了有效的数值实验。

Abstract: For solving the fractional time conservation law that may have discontinuities such as shock waves, when the fractional order γ is close to 0, there is no effective method to solve the fractional nonline-ar discrete system. Taking the time fractional Burgers equation as an example, this paper applies the multi-grid iterative method to solve the problem. For the convection term, a flux limiter is adopted to make the new numerical flux become a higher order flux in the smooth area and a lower order flux near the discontinuity, thus achieving a higher order precision. The effective numerical experiments are carried out with different values of γ and different initial boundary values.

1. 引言

时间分数阶标量守恒律方程

$\frac{{\partial }^{\gamma }u\left(x,t\right)}{\partial t^{\gamma }}+\frac{\partial \left(f\left(u\right)\right)}{\partial x}=0,\left(x,t\right)\in \left[a,b\right]\times \left(0,T\right]$ (1)

当 $f\left(u\right)=u^2$ 时，我们称其为时间分数阶Burgers方程。

本文考虑下述时间分数阶Burgers方程：

$\frac{{\partial }^{\gamma }u\left(x,t\right)}{\partial t^{\gamma }}+\frac{\partial {\left(u\left(x,t\right)\right)}^2}{\partial x}=0$ (2)

其中 $a\le x\le b$， $0<t\le T$。

满足初值条件

$u\left(x,0\right)=g\left(x\right),a\le x\le b$

以及合适的边界条件，其中 $\gamma \in \left(0,1\right)$， $g\left(x\right)$ 为给定函数，

$\frac{{\partial }^{\gamma }u\left(x,t\right)}{\partial t^{\gamma }}=\frac{1}{\Gamma \left(1-\gamma \right)}{\displaystyle {\int }_0^t{\left(t-\tau \right)}^{-\gamma }\frac{\partial u\left(x,\tau \right)}{\partial \tau }\text{d}\tau }$

为 $\gamma$ 阶Caputo分数阶导数， $\Gamma \left(x\right)$ 为Gamma函数。

近十年来，许多学者运用不同的数值方法对含黏时间分数阶Burgers方程进行了求解，例如，Duangpan运用了有限积分法结合移位Chebyshev多项式 [1]，Akram等提出了有限差分格式 [2]，Esen和Tasbozan提出了三次b样条有限元搭配方法 [3] 等等。

2. 离散

定义空间步长 $h=\left(b-a\right)/\left(N+1\right)$，空间网格节点 $x_i=a+ih$，其中 $i=0,1,\cdots ,N$ ；定义空间网格半节点 $x_{i+\frac{1}{2}}=a+\left(i+1/2\right)h$，得到 $N+1$ 个网格单元 $\left[x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}\right]$，其中 $i=-1,0,1,\cdots ,N$。定义时间步长 $\Delta t=T/M$，时间网格节点 $t_k=k\Delta t$，其中 $k=0,1,\cdots ,M$。

对于函数 $\mathbb{F}$，使用如下简写形式：

${\mathbb{F}}_{i+\frac{1}{2},n}=\mathbb{F}\left(x_i,t_n\right)$ ${\mathbb{F}}_{i+\frac{1}{2}}=\mathbb{F}\left(x_{i+\frac{1}{2}}\right)$

在 $t=t_n$ 处计算式(1.1)并在单元 $\left[x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}\right]$ 上对方程进行积分，得到

$h\frac{{\partial }^{\gamma }{\bar{u}}_i\left(t_n\right)}{\partial t^{\gamma }}+\left(f{\left(u\right)}_{i+\frac{1}{2},n}-f{\left(u\right)}_{i-\frac{1}{2},n}\right)=0,i=0,1,\cdots ,N$ (3)

其中 ${\bar{u}}_i\left(t\right)={\displaystyle {\int }_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}u\left(x,t\right)\text{d}x/h}$ 表示 $u\left(x,t\right)$ 在区间 $\left(x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}\right)$ 上的平均值。

定义 $U_i^n={\bar{u}}_i\left(t_n\right)$，对式(3)的左端第一项运用 $L_1$ 公式逼近，取 $a_k={\left(k+1\right)}^{1-\alpha }-k^{1-\alpha }$ [4]：

${h\frac{{\partial }^{\gamma }{\bar{u}}_i\left(t_n\right)}{\partial t^{\gamma }}|}_{t_n}\approx \frac{h\Delta t^{-\gamma }}{\Gamma \left(2-\gamma \right)}\left(U_i^n-{\displaystyle \underset{k=1}{\overset{n-1}{\sum }}\left(a_{n-k-1}-a_{n-k}\right)U_i^k}-a_{n-1}{U}_i^0\right).$ (4)

通量近似为

$F_{j+\frac{1}{2},n}^n=f{\left(u\right)}_{i+\frac{1}{2},n}\approx F\left(U_j^n,U_{j+1}^n\right)$ (5)

$F_{j-\frac{1}{2},n}^n=f{\left(u\right)}_{i-\frac{1}{2},n}\approx F\left(U_{j-1}^n,U_j^n\right)$ (6)

结合式(3)、(4)、(5)及(6)，得到最终的离散形式为

$\frac{h\Delta t^{-\gamma }}{\Gamma \left(2-\gamma \right)}{\left(U_i^n-{\displaystyle \underset{k=1}{\overset{n-1}{\sum }}\left(a_{n-k-1}-a_{n-k}\right)U_i^k}-a_{n-1}U\right)}_i^0+F_{i+\frac{1}{2}}^n-F_{i-\frac{1}{2}}^n=0$ (7)

3. 通量限制器

为了使问题的解达到更高阶精度，我们首先选择一个在光滑区域表现很好的高阶通量和一个在间断附近表现很好的低阶通量，然后将数值通量重新定义为低阶通量和高阶通量的线性组合，使得新的数值通量在光滑区域变为高阶通量而在间断附近变为低阶通量：

$F_{j\pm \frac{1}{2}}^n=F_{j\pm \frac{1}{2}}^{n,L}+{\varphi }_{j\pm \frac{1}{2}}\left(F_{j\pm \frac{1}{2}}^{n,H}-F_{j\pm \frac{1}{2}}^{n,L}\right)$

式中 $F_{j\pm \frac{1}{2}}^{n,L}$ 为低阶通量， $F_{j\pm \frac{1}{2}}^{n,H}$ 为高阶通量； ${\varphi }_{j\pm \frac{1}{2}}$ 为通量限制器，是变差比 $r_{j\pm \frac{1}{2}}$ 的函数，即

${\varphi }_{j\pm \frac{1}{2}}={\varphi }_{j\pm \frac{1}{2}}\left(r_{j\pm \frac{1}{2}}\right)$

$r_{j\pm \frac{1}{2}}$ 的取值与流动方向相关，对于均匀网格，定义

$r_{j+\frac{1}{2}}=\{\begin{array}{l}\frac{U_j^n-U_{j-1}^n}{U_{j+1}^n-U_j^n},U_j^n\ge 0,\\ \frac{U_{j+2}^n-U_{j+1}^n}{U_{j+1}^n-U_j^n},U_j^n\le 0,\end{array}$

$r_{j-\frac{1}{2}}=\{\begin{array}{l}\frac{U_{j-1}^n-U_{j-2}^n}{U_j^n-U_{j-1}^n},U_{j-1}^n\ge 0,\\ \frac{U_{j+1}^n-U_j^n}{U_j^n-U_{j-1}^n},U_{j-1}^n\le 0,\end{array}$

可以看到，当 ${\varphi }_{j\pm \frac{1}{2}}=0$ 时，通量为低阶格式；当 ${\varphi }_{j\pm \frac{1}{2}}=1$ 时，通量为高阶格式。

为保证格式稳定性， $\varphi \left(1\right)=1$ 是必要的；同时，格式需要满足二阶精度，即对于限制器 $\varphi$ 有如下约束条件(若 $r<0$，则令 $\varphi \left(r\right)=0$ )：

$\{\begin{array}{l}r\le \varphi \left(r\right)\le 1,0<r\le 1,\\ 1\le \varphi \left(r\right)\le r,r\ge 1.\end{array}$

我们把Lax-Friedrichs通量作为低阶通量：

$F_{j-\frac{1}{2}}^{n,L}=\frac{f\left(U_{j-1}^n\right)+f\left(U_j^n\right)}{2}-\frac{{\alpha }_{j-\frac{1}{2}}}{2}\left(U_j^n-U_{j-1}^n\right)$

$F_{j+\frac{1}{2}}^{n,L}=\frac{f\left(U_j^n\right)+f\left(U_{j+1}^n\right)}{2}-\frac{{\alpha }_{j+\frac{1}{2}}}{2}\left(U_{j+1}^n-U_j^n\right)$

如果 ${\alpha }_{j\pm \frac{1}{2}}$ 取常数 $\alpha$，令 $\alpha =\underset{j}{\kappa \max }\left|f^{\prime }\left(U_j^n\right)\right|$，则称Lax-Friedrichs通量为全局Lax-Friedrichs通量；如果 ${\alpha }_{j+\frac{1}{2}}=\kappa \max \left|f^{\prime }\left(U_j^n\right),f^{\prime }\left(U_{j+1}^n\right)\right|$，则称其为局部Lax-Friedrichs通量，其中 $\kappa \ge 1$，是一个常数标量。

使用Lax-Wendroff通量为高阶通量：

$F_{j-\frac{1}{2}}^{n,H}=\frac{f\left(U_{j-1}^n\right)+f\left(U_j^n\right)}{2}-\frac{\Delta t}{2h}{f}^{\prime }\left(\frac{U_{j-1}^n+U_j^n}{2}\right)\left(f\left(U_j^n\right)-f\left(U_{j-1}^n\right)\right)$

$F_{j+\frac{1}{2}}^{n,H}=\frac{f\left(U_j^n\right)+f\left(U_{j+1}^n\right)}{2}-\frac{\Delta t}{2h}{f}^{\prime }\left(\frac{U_j^n+U_{j+1}^n}{2}\right)\left(f\left(U_{j+1}^n\right)-f\left(U_j^n\right)\right)$

下面列举出一些常见的通量限制器 [5]：

1) Chakravarthy-Osher限制器：

$\varphi \left(r\right)=\max \left[0,\min \left(r,\beta \right)\right],1\le \beta \le 2$

当 $\beta =1$ 时，称其为minmod限制器；

2) Koren限制器：

$\varphi \left(r\right)=\max \left[0,\min \left(2r,\left(2+r\right)/3\right),2\right],$

3) Sweby限制器：

$\varphi \left(r\right)=\max \left[0,\min \left(\beta r,1\right),\min \left(r,\beta \right)\right],1\le \beta \le 2.$

当 $\beta =2$ 时，称其为Superbee限制器；

4) OSPRE限制器：

$\varphi \left(r\right)=\frac{3\left(r^2+r\right)}{2\left(r^2+r+1\right)},$

5) van Leer限制器：

$\varphi \left(r\right)=\max \left[0,\frac{2r}{1+r}\right],$

限制器 $\varphi \left(r\right)$ 一般是对称的，即满足以下等式：

$\varphi \left(r^{-1}\right)=\frac{\varphi \left(r\right)}{r}$

4. 多重网格法

本文使用多重网格迭代方法，由于空间系统是强非线性的，我们的多重网格方法将以线性多重网格方法的非线性推广——FAS多重网格方法为基础 [6]。

假设三个正整数 $N,N_0,level$ 满足式：

$N=N_0\times 2^{level-1}$

构建空间网格层次结构： ${\left\{{\Gamma }_k\right\}}_{k=1}^{level}$，假设最细网格 ${\Gamma }_1$ 的网格大小为 $h_1=h$，对于 $k=2,3,\cdots ,level$， ${\Gamma }_k$ 的网格大小为 $h_k=2^{k-1}h$，其中 $N_k=N/2^{k-1}$ 表示 ${\Gamma }_k$ 的网格数量， $k=1,2,\cdots ,level$，下文中与 ${\Gamma }_k$ 有关的量以上标 $h_k$ 的形式书写，例如 $U^{h_k}$。

在多重网格法中，我们通过求解粗网格 ${\Gamma }_k$ 上的离散方程对迭代解进行反复修正：

$A^{h_k}\left(U^{h_k}\right)=r^{h_k}$

其中 $A^{h_k}$ 为非线性算子。

对式(7)采用Richardson迭代方法求解：

${\hat{U}}_i^{h_k}={\tilde{U}}_i^{h_k}+\frac{1}{{\lambda }^{h_k}}\left[b_i-\frac{\Delta t^{\gamma }\Gamma \left(2-\gamma \right)}{h_k}\left(F\left({\tilde{U}}_i^{h_k},{\tilde{U}}_{i+1}^{h_k}\right)-F\left({\tilde{U}}_{i-1}^{h_k},{\tilde{U}}_i^{h_k}\right)\right)-{\tilde{U}}_i^{h_k}\right]$ (8)

其中 $b_i={\displaystyle \underset{k=1}{\overset{n-1}{\sum }}\left(a_{n-k-1}-a_{n-k}\right)U_i^k}+a_{n-1}{U}_i^0$， ${\hat{U}}^{h_k}$ 和 ${\tilde{U}}^{h_k}$ 分别为一次迭代过程中的新近似解和旧近似解， ${\lambda }^{h_k}$ 为一个能够保证迭代方法收敛性的合适的正数。

我们定义延拓和限制算子来实现相邻网格之间的数据传输，其中，延拓算子 $I_{h_{k+1}}^{h_k}$ 实现 $R^{N_{k+1}}\to R^{N_k}$ 的数据传输，即

$U^{h_k}=I_{h_{k+1}}^{h_k}{U}^{h_{k+1}}$

定义为：

$\{\begin{array}{l}{U}_{2i}^{h_k}=U_i^{h_{k+1}}\\ U_{2i+1}^{h_k}=\frac{1}{2}\left(U_i^{h_{k+1}}+U_{i+1}^{h_{k+1}}\right)\end{array}$

限制算子 $I_{h_k}^{h_{k+1}}$ 实现 $R^{N_k}\to R^{N_{k+1}}$ 的数据传输，即

$U^{h_{k+1}}=I_{h_k}^{h_{k+1}}{U}^{h_k}$

定义为：

$U_i^{h_{k+1}}=\frac{1}{4}\left(U_{2i-1}^{h_k}+U_{2i}^{h_k}+U_{2i+1}^{h_k}\right)$

以下提出求解方程的多重网格迭代方法步骤：

定义 ${\hat{U}}_i^{h_k}=FAS\left(k,{\bar{U}}_i^{h_k},r^{h_k}\right)$ ：

步1：在最网格上计算非线性离散方程 $A^{h_k}\left(U^{h_k}\right)=f^{h_k}$，得到方程近似解 ${\bar{U}}^{h_k}$ ；

步2：在粗网格上计算非线性残差方程 $r^{h_k}=f^{h_k}-A^{h_k}\left({\bar{U}}^{h_k}\right)$ ；

步3：得到新的右端项 $f^{h_{k+1}}=A^{h_{k+1}}\left(I_k^{k+1}{\bar{U}}^{h_k}\right)+\left(I_k^{k+1}{r}^{h_k}\right)$ ；

步4：在粗网格上对近似解进行修正 ${\hat{U}}^{h_k}={\bar{U}}^{h_k}+I_{h_{k+1}}^{h_k}\left(U^{h_{k+1}}-I_{h_k}^{h_{k+1}}{\bar{U}}^{h_k}\right)$。

5. 数值实验

例5.1设计算域 $\left[a,b\right]\times \left(0,T\right]=\left[0,1\right]\times \left(0,0.2\right]$，初值条件为 $u\left(x,0\right)=\sin \left(2\pi x\right)$，边界条件为 $u\left(a\right)=u\left(b\right)=0$，取 $N=128$， $M=100$，分别使用OSPRE limiter (见图1)和Koren limiter (见图2)求解在 $\gamma =0.8,0.6,0.4,0.2$ 时方程(2)的解：

Figure 1. OSPRE limiter

图1. OSPRE限制器

Figure 2. Koren limiter

图2. Koren限制器

例5.2设计算域 $\left[a,b\right]\times \left(0,T\right]=\left[-1,3\right]\times \left(0,0.2\right]$，初值条件为 $u\left(x,0\right)=-\sin \left(\pi x\right)$，边界条件为 $u\left(a\right)=u\left(b\right)=0$，取 $N=256$， $M=100$，分别使用Koren limiter (见图3)和van Leer limiter (见图4)求解在 $\gamma =0.8,0.6,0.4,0.2$ 时方程(2)的解：

Figure 3. Koren limiter

图3. Koren限制器

Figure 4. Van Leer limiter

图4. Van Leer限制器

例5.3设计算域 $\left[a,b\right]\times \left(0,T\right]=\left[-1,3\right]\times \left(0,0.2\right]$，初值条件为 $u\left(x,0\right)=\{\begin{array}{l}0,x\ge 0\\ 1,x<0\end{array}$，边界条件为 $u\left(a\right)=1$， $u\left(b\right)=0$，取 $N=256$， $M=128$，分别使用OSPRE limiter (见图5)和van Leer limiter (见图6)求解在 $\gamma =0.8,0.6,0.4,0.2$ 时方程(2)的解：

Figure 5. OSPRE limiter

图5. OSPRE限制器

Figure 6. Van Leer limiter

图6. Van Leer限制器

6. 总结

本文运用多重网格迭代方法求解时间分数阶Burgers方程，对于对流项，引入了通量限制器，在三种不同的初值和边界条件下，选用了两种通量限制器，分别在 $\gamma =0.8,0.6,0.4,0.2$ 时方程的解，得到了较好的结果。