二维线性传输方程满足两个守恒律的数值格式

The Numerical Scheme Satisfying Two Conservation Laws for Two Dimensional Linear Advection Equations

摘要: 在本文中，我们针对二维线性传输方程设计了满足两个守恒律的数值格式。该格式不仅能保持数值解守恒，同时能保持数值能量守恒。通过数值算例验证格式的有效性，数值结果表明该格式在远离极值点的区域内具有误差相互抵消的超收敛性质，并且能够很好的保持解的结构。

Abstract: In this paper, we develop an improved numerical scheme satisfying two conservation laws for two-dimensional linear advection equations, which satisfying both the numerical solution and numerical energy conservative. The numerical results show that the scheme has the super-con- vergence property away from the extreme points, and can keep the structure of the solution well.

1. 引言

本文对二维线性传输方程设计能够满足数值解和数值能量同时守恒的数值格式，这种满足多个守恒律的格式与传统的守恒型格式不同，在格式中不仅计算数值解，而且要计算数值能量，并且在格式中有多个数值实体参与格式的计算。

文献 [1] 中对线性传输方程设计了满足两个守恒律的数值格式；文献 [2] [3] 中对KdV方程设计了满足两个守恒律的数值格式；文献 [4] 中针对一维守恒型方程组设计了熵耗散格式；文献 [5] 中对一维线性传输方程设计了满足两个守恒律的数值格式，并将此思想推广到了二维的情况；文献 [6] 中设计了线性传输方程满足三个守恒律的数值格式，数值实验结果表明这类满足多个守恒律的数值格式比传统的格式有更高的精度且长时间的数值模拟效果非常好。

本文将这种满足多个守恒律的思想应用到二维线性传输方程，与文献 [5] 不同的是，本文设计的格式，采用分裂算子法，将二维线性传输方程写成了两个一维方程，对其分别设计了满足两个守恒律的数值格式，在格式中数值解和数值能量都能守恒，并且数值能量不是被动的守恒律，而是作为一个数值实体参与了格式的计算。数值结果表明该格式具有误差互相抵消的超收敛性质，并且在长时间的数值模拟中能很好的保持解的结构。

本文的结构安排如下，第一节是引言；第二节是格式的具体描述；第三节是格式的可行性分析，第四节给出了数值算例，验证了算法的有效性；第五节是结论。

2. 格式的描述

我们考虑二维线性传输方程：

$u_t+u_x+u_y=0$ (2.1)

它的解满足以下形式的守恒律：

$U{\left(u\right)}_t+U{\left(u\right)}_x+U{\left(u\right)}_y=0$ (2.2)

在本文中取 $U\left(u\right)=u^2$，我们称之为能量。

在数值格式中，我们采用正规网格，令 $x_i=i{h}_1$， $t_n=n\tau$， $x_{i\pm \frac{1}{2}}=\left(i\pm \frac{1}{2}\right)h_1$， $y_j=j{h}_2$， $y_{j\pm \frac{1}{2}}=\left(j\pm \frac{1}{2}\right)h_2$。其中 $h_1$ 和 $h_2$ 分别为X轴方向和Y轴方向上的空间步长， $\tau$ 是时间步长， ${\lambda }_1=\frac{\tau }{h_1}$， ${\lambda }_2=\frac{\tau }{h_2}$ 分别为X轴方向和Y轴方向上的网格步长比。

方程(2.1)的精确解为 $u\left(x,y,t\right)$，方程(2.2)的精确解为 $U\left(x,y,t\right)=u^2\left(x,y,t\right)$，在数值格式中我们不仅

计算了数值解，还计算了数值能量。在网格 $\left(x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}\right)\times \left(y_{j-\frac{1}{2}},y_{j+\frac{1}{2}}\right)$ 上， $t_n$ 时刻的数值解记为 $u_{i,j}^n$，数值

能量记为 $U_{i,j}^n$，它们分别为精确解和精确能量的网格平均的近似：

$u_{ij}^n\simeq \frac{1}{h_1{h}_2}{\displaystyle {\int }_{y_{j-\frac{1}{2}}}^{y_{j+\frac{1}{2}}}{\displaystyle {\int }_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}u\left(x,y,t_n\right)\text{d}x\text{d}y}}$ (2.3)

$U_{ij}^n\simeq \frac{1}{h_1{h}_2}{\displaystyle {\int }_{y_{j-\frac{1}{2}}}^{y_{j+\frac{1}{2}}}{\displaystyle {\int }_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}U\left(u\left(x,y,t_n\right)\right)\text{d}x\text{d}y}}$ (2.4)

对于二维方程(2.1)和(2.2)，我们采用Godunov型分裂算子法 [7] [8]，将其写成两部分A和B，

$\text{A:}\{\begin{array}{l}{u}_t+u_x=0;\\ U_t+U_x=0;\\ u\left(x,y,0\right)=u_0\left(x,y\right).\end{array}$ (2.5)

$\text{B:}\{\begin{array}{l}{u}_t+u_y=0;\\ U_t+U_y=0;\\ u\left(x,y,0\right)=u\left(x,y,t\right).\end{array}$ (2.6)

其中 $U=u^2$。

下面，我们将讨论如何根据 $t_n$ 层上的数值解 $u_{i,j}^n$ 和数值能量 $U_{i,j}^n$ 来计算 $t_{n+1}$ 层上的数值解 $u_{i,j}^{n+1}$ 和数值能量 $U_{i,j}^{n+1}$。

2.1. 守恒方程A的数值格式

对于守恒方程A，其数值格式为Godunov型，我们按照重构、发展、网格平均三步设计格式，但是与传统Godunov型格式不同的是，我们的格式不仅计算了数值解，还计算了数值能量。

第一步(重构)：在 $t_n$ 层上，对数值解 $u_{i,j}^n$ 和数值能量 $U_{i,j}^n$ 进行重构：

$R\left(x,t;u_{i,j}^n,U_{i,j}^n\right)=u_{i,j}^n+s_{i,j}^n\left(x-x_i\right),x\in \left(x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}\right)$ (2.7)

显然， $R\left(x,t;u_{i,j}^n,U_{i,j}^n\right)$ 满足：

$\frac{1}{h_1}{\displaystyle {\int }_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}R\left(x,t;u_{i,j}^n,U_{i,j}^n\right)\text{d}x}=u_{i,j}^n$ (2.8)

为了保持两个守恒律同时成立，我们还要求重构函数 $R\left(x,t;u_{i,j}^n,U_{i,j}^n\right)$ 满足：

$\frac{1}{h_1}{\displaystyle {\int }_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}{\left(R\left(x,t;u_{i,j}^n,U_{i,j}^n\right)\right)}^2}\text{d}x=U_{i,j}^n$ (2.9)

即重构函数数值能量的网格平均和数值能量相等。根据(2.9)式，我们可以算出重构函数的斜率 $s_{i,j}^n$ ：

$s_{i,j}^n=\mathrm{sgn}\left(u_{i+1,j}^n-u_{i-1,j}^n\right)\sqrt{\frac{12\left(U_{i,j}^n-{\left(u_{i,j}^n\right)}^2\right)}{{\left(h_1\right)}^2}}$ (2.10)

从这里可以看出，与传统差分格式不同的是，斜率 $s_{i,j}^n$ 作为一个自由度，是通过数值能量的守恒得到的。

第二步(发展)：以 $R\left(x,t;u_{i,j}^n,U_{i,j}^n\right)$ 为 $t_n$ 层上的初值，求解初值问题：

$\{\begin{array}{l}{v}_t+v_x=0;\\ v\left(x,y_j;t_n\right)=R\left(x,t;u_{i,j}^n,U_{i,j}^n\right).\end{array}-\infty <x<\infty ,t_n<t<t_{n+1}$ (2.11)

其精确解为： $v\left(x,y_j;t_n\right)=R\left(x-\left(t-t_n\right),u_{i,j}^n,U_{i,j}^n\right)$。

第三步(网格平均)：在 $t_{n+1}$ 层上，A部分的数值解 $u_{i,j}^{A,n+1}$ 和数值能量 $U_{i,j}^{A,n+1}$ 分别为：

$u_{i,j}^{A,n+1}=\frac{1}{h_1}{\displaystyle {\int }_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}v\left(x,y_j;t_{n+1}\right)\text{d}x}$ (2.12)

$U_{i,j}^{A,n+1}=\frac{1}{h_1}{\displaystyle {\int }_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}{\left(v\left(x,y_j;t_{n+1}\right)\right)}^2}$ (2.13)

在实际的计算中，我们采用积分形式来计算数值解 $u_{i,j}^{A,n+1}$ 和数值能量 $U_{i,j}^{A,n+1}$ (见文献 [3] )。

因此数值解 $u_{i,j}^{A,n+1}$ 的数值格式为：

$u_{i,j}^{A,n+1}=u_{i,j}^n-{\lambda }_1\left({\hat{f}}_{i+\frac{1}{2},j}^n-{\hat{f}}_{i-\frac{1}{2},j}^n\right)$ (2.14)

其中，数值流函数为：

${\hat{f}}_{i+\frac{1}{2},j}^n=\frac{1}{\tau }{\displaystyle {\int }_{t_n}^{t_{n+1}}v\left(x_{i+\frac{1}{2}},y_j,t\right)\text{d}t}=u_{i,j}^n+\frac{s_{i,j}^n}{2}\left(h_1-\tau \right)$ (2.15)

上式中的数值流函数是精确求解。如果方程是非线性传输方程，数值流函数也可以采用数值解法精确求出。

数值能量 $U_{i,j}^{A,n+1}$ 的数值格式为：

$U_{i,j}^{A,n+1}=U_{i,j}^n-{\lambda }_1\left({\hat{F}}_{i+\frac{1}{2},j}^n-{\hat{F}}_{i-\frac{1}{2},j}^n\right)$ (2.16)

其中，数值能量流函数为：

$\begin{array}{c}{\hat{F}}_{i+\frac{1}{2},j}^n=\frac{1}{\tau }{\displaystyle {\int }_{t_n}^{t_{n+1}}{\left(v\left(x_{i+\frac{1}{2}},y_j,t\right)\right)}^2}\text{d}t\\ ={\left(u_{i,j}^n\right)}^2+u_{i,j}^n{s}_{i,j}^n\left(h_1-\tau \right)+{\left(s_{i,j}^n\right)}^2{\left(h_1\right)}^2\frac{4{\lambda }_1{}^2-6{\lambda }_1+3}{12}\end{array}$ (2.17)

2.2. 守恒方程B的数值格式

对于守恒方程B部分的数值格式，由分裂算子法可知，以守恒方程 A 得到的数值解 $u_{i,j}^{A,n+1}$ 和数值能量 $U_{i,j}^{A,n+1}$ 为初值，按照重构、发展、网格平均三步进行设计格式。在这里我们同样不仅计算数值解，也要计算数值能量。

第一步(重构)：对守恒方程A得到的数值解 $u_{i,j}^{A,n\text{+}1}$ 和数值能量 $u_{i,j}^{A,n\text{+}1}$ 进行重构。

$R\left(y;u_{i,j}^{A,n\text{+}1},U_{i,j}^{A,n\text{+}1}\right)=u_{i,j}^{A,n\text{+}1}+{\sigma }_{i,j}^{n\text{+}1}\left(y-y_j\right),y\in \left(y_{j-\frac{1}{2}},y_{j+\frac{1}{2}}\right)$ (2.18)

同样地在我们设计的数值格式中，要求重构函数 $R\left(y;u_{i,j}^{A,n\text{+}1},U_{i,j}^{A,n\text{+}1}\right)$ 满足：

$\{\begin{array}{l}\frac{1}{h_2}{\displaystyle {\int }_{y_{j-\frac{1}{2}}}^{y_{j+\frac{1}{2}}}R\left(y;u_{i,j}^{A,n\text{+}1},U_{i,j}^{A,n\text{+}1}\right)\text{d}y}=u_{i,j}^{A,n\text{+}1},\\ \frac{1}{h_2}{\displaystyle {\int }_{y_{j-\frac{1}{2}}}^{y_{j+\frac{1}{2}}}{\left(R\left(y;u_{i,j}^{A,n\text{+}1},U_{i,j}^{A,n\text{+}1}\right)\right)}^2}\text{d}y=U_{i,j}^{A,n\text{+}1},\end{array}$ (2.19)

其中，(2.19)式的第一个式子保持了重构函数数值解的守恒，第二个式子保持了重构函数数值能量的守恒。

根据(2.19)式的第二个式子，我们可以解出斜率 ${\sigma }_{i,j}^{n\text{+}1}$ ：

${\sigma }_{i,j}^{n\text{+}1}=\mathrm{sgn}\left(u_{i,j+1}^{A,n\text{+}1}-u_{i,j-1}^{A,n\text{+}1}\right)\sqrt{\frac{12\left(U_{i,j}^{A,n\text{+}1}-{\left(u_{i,j}^{A,n\text{+}1}\right)}^2\right)}{{\left(h_2\right)}^2}}$ (2.20)

上式表明，斜率 ${\sigma }_{i,j}^{n\text{+}1}$ 不仅与数值解和数值能量有关，而且数值解和数值能量之间是一种非线性关系。

第二步(发展)：以 $R\left(y;u_{i,j}^{A,n\text{+}1},U_{i,j}^{A,n\text{+}1}\right)$ 为初值，求解初值问题：

$\{\begin{array}{l}{v}_t+v_y=0;\\ v\left(x_i,y;t_n\right)=R\left(y,t;u_{i,j}^{A,n\text{+}1},U_{i,j}^{A,n\text{+}1}\right).\end{array}-\infty <y<\infty ,t_n<t<t_{n+1}$ (2.21)

其精确解为： $v\left(x_i,y;t_n\right)=R\left(y-\left(t-t_n\right),u_{i,j}^{A,n\text{+}1},U_{i,j}^{A,n\text{+}1}\right)$。

第三步(网格平均)：在 $t_{n+1}$ 层上，我们采用积分形式计算数值解 $u_{i,j}^{n+1}$ 和数值能量 $U_{i,j}^{n+1}$。因此，数值解 $u_{i,j}^{n+1}$ 和数值能量 $U_{i,j}^{n+1}$ 的数值格式分别为：

$\{\begin{array}{l}{u}_{i,j}^{n+1}=u_{i,j}^{A,n+1}-{\lambda }_2\left({\hat{g}}_{i,j+\frac{1}{2}}^n-{\hat{g}}_{i,j-\frac{1}{2}}^n\right),\\ U_{i,j}^{n+1}=U_{i,j}^{A,n+1}-{\lambda }_2\left({\hat{G}}_{i,j+\frac{1}{2}}^n-{\hat{G}}_{i,j-\frac{1}{2}}^n\right),\end{array}$ (2.22)

其中，数值流函数和数值能量流函数分别为：

${\hat{g}}_{i,j+\frac{1}{2}}^{n\text{+}1}=\frac{1}{\tau }{\displaystyle {\int }_{t_n}^{t_{n+1}}v\left(x_i,y_{j+\frac{1}{2}},t\right)\text{d}t}=u_{i,j}^{A,n\text{+}1}+\frac{{\sigma }_{i,j}^{n\text{+}1}}{2}\left(h_2-\tau \right),$ (2.23)

$\begin{array}{c}{\hat{G}}_{i,j+\frac{1}{2}}^{n\text{+}1}=\frac{1}{\tau }{\displaystyle {\int }_{t_n}^{t_{n+1}}{\left(v\left(x_i,y_{j+\frac{1}{2}},t\right)\right)}^2}\text{d}t\\ ={\left(u_{i,j}^{A,n\text{+}1}\right)}^2+u_{i,j}^{A,n\text{+}1}{\sigma }_{i,j}^{n\text{+}1}\left(h_2-\tau \right)+{\left({\sigma }_{i,j}^{n\text{+}1}\right)}^2{\left(h_2\right)}^2\frac{4{\lambda }_2{}^2-6{\lambda }_2+3}{12}\end{array}$ (2.24)

综上所述，对于方程(2.1)和(2.2)，由A部分的数值格式(2.14)、(2.16)和B部分的数值格式(2.22)可知，数值解和数值能量的格式可以写成如下形式：

$u_{i,j}^{n+1}=u_{i,j}^n-{\lambda }_1\left({\hat{f}}_{i+\frac{1}{2},j}^n-{\hat{f}}_{i-\frac{1}{2},j}^n\right)-{\lambda }_2\left({\hat{g}}_{i,j+\frac{1}{2}}^n-{\hat{g}}_{i,j-\frac{1}{2}}^n\right),$ (2.25)

$U_{i,j}^{n+1}=U_{i,j}^n-{\lambda }_1\left({\hat{F}}_{i+\frac{1}{2},j}^n-{\hat{F}}_{i-\frac{1}{2},j}^n\right)-{\lambda }_2\left({\hat{G}}_{i,j+\frac{1}{2}}^n-{\hat{G}}_{i,j-\frac{1}{2}}^n\right).$ (2.26)

3. 格式的可行性

在数值格式中，重构函数的斜率 $s_{i,j}^n$， ${\sigma }_{i,j}^n$ 是通过计算(2.10)式和(2.20)式得到的，从中我们可以看到计算过程中涉及到开根号运算，定理3.1保证了开根号运算的可行性。

定理 3.1：如果 $U_{i,j}^0\ge {\left(u_{i,j}^0\right)}^2$， $U_{i,j}^{A,0}\ge {\left(u_{i,j}^{A,0}\right)}^2$，那么通过格式(2.25)和(2.26)式计算所得的数值解和数值

能量满足关系： $U_{i,j}^n\ge {\left(u_{i,j}^n\right)}^2,U_{i,j}^{A,n}\ge {\left(u_{i,j}^{A,n}\right)}^2$。

证明：我们采用数学归纳法来证明 $U_{i,j}^{A,n}\ge {\left(u_{i,j}^{A,n}\right)}^2$。

当n = 0时，根据已知条件有：

$U_{i,j}^0\ge {\left(u_{i,j}^0\right)}^2$ (3.1)

当n > 0时，假设t_n层上满足关系 $U_{i,j}^{A,n}\ge {\left(u_{i,j}^{A,n}\right)}^2$，下面证明 $t_{n+1}$ 层上满足关系 $U_{i,j}^{A,n+1}\ge {\left(u_{i,j}^{A,n+1}\right)}^2$。

根据上述假设，可知(2.10)式具有实根，因此重构函数 $R\left(x,t;u_{i,j}^n,U_{i,j}^n\right)$ 是存在的，根据(2.14)到(2.17) 式可以计算出 $u_{i,j}^{A,n+1}$ 和 $U_{i,j}^{A,n+1}$，即：

$\{\begin{array}{l}{u}_{i,j}^{A,n+1}=\frac{1}{h_1}{\displaystyle {\int }_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}v\left(x,y_j;t_{n+1}\right)\text{d}x},\\ U_{i,j}^{A,n+1}=\frac{1}{h_1}{\displaystyle {\int }_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}{\left(v\left(x,y_j;t_{n+1}\right)\right)}^2}\text{d}x,\end{array}$ (3.2)

根据Jessen不等式可得：

$\frac{1}{h_1}{\displaystyle {\int }_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}U\left(v\left(x,y_j,t_{n+1}\right)\right)\text{d}x}\ge U\left(\frac{1}{h_1}{\displaystyle {\int }_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}v\left(x,y_j,t_{n+1}\right)}\right)\text{d}x,$ (3.3)

因此，

$U_{i,j}^{A,n+1}\ge {\left(u_{i,j}^{A,n+1}\right)}^2$ (3.4)

同理，可以得到：

$U_{i,j}^{n+1}\ge {\left(u_{i,j}^{n+1}\right)}^2$ (3.5)

证毕。

4. 数值算例

在这一部分，我们将通过数值算例来观察格式的精度和数值模拟效果。

算例 4.1: 对于二维线性传输方程 $u_t+u_x+u_y=0$，令初值为：

$u\left(x,y,0\right)=\sin \left(2\pi \left(x+y\right)\right),0\le x\le 1$ (4.1)

周期边界条件。该方程的精确解为 $u\left(x,y,t\right)=\sin \left(2\pi \left(x+y-2t\right)\right)$。令网格步长比为 ${\lambda }_1={\lambda }_2=0.2$，分别用40，80，160，320，640，1280个网格计算到t = 1时刻的数值解，并对数值结果做 $L^1$ 误差和 $L^{\infty }$ 误差及收敛阶的分析。

通过表1可以看出， $L^1$ 误差是二阶的，而 $L^{\infty }$ 误差是介于一阶和二阶之间的，这是因为在极值点附近，为了保持解的稳定性，精度会降低一阶。

为了研究在远离极值点的区域内数值误差的精度，我们考虑远离极值点的集合 $\Omega =\left[0,\text{1}/6\right]\cup \left[\text{1}/3,\text{2}/3\right]\cup \left[\text{5}/6,\text{1}\right]$，同样对数值结果做 $L^1$ 误差和 $L^{\infty }$ 误差及收敛阶的分析。

通过表2可以看出，在远离极值点的区域内，格式可以达到三阶精度。这是因为在计算过程中，虽然采用线性重构，但是由于本文设计的格式同时计算了数值能量，因此数值解的误差不是简单的累加，误差在每一步的计算中会相互抵消。因此在远离极值点的区域内格式的精度能达到三阶，表现出超收敛性质。

Table 1. The error and convergence order of interval [0, 1] ×[0, 1]

表1. [0, 1] ×[0, 1]区间的误差和收敛阶

Table 2. The error and convergence order of interval Ω

表2. Ω区间的误差和收敛阶

算例4.2：对于方程(2.1)式，取初值为：

$\{\begin{array}{l}u\left(x,y,0\right)=\exp \left(-\frac{1}{1-16{\left(x-\frac{1}{2}\right)}^2}\right)\times \exp \left(-\frac{1}{1-16{\left(y-\frac{1}{2}\right)}^2}\right)\\ \left(x,y\right)\in \left(\frac{1}{4},\frac{3}{4}\right)\times \left(\frac{1}{4},\frac{3}{4}\right)\\ \\ u\left(x,y,0\right)=0,\left(x,y\right)\in \left(\frac{1}{4},\frac{3}{4}\right)\times \left(\frac{1}{4},\frac{3}{4}\right)\end{array}$ (4.2)

无穷边界，网格步长为1/80，网格步长比为0.1。图1中的(a)图为初始时刻的解，(b)图为计算了t = 0.2 时刻的数值解。从图1(b)可以看出，我们设计的格式能很好的保持数值解的结构。

5. 结论

本文采用分裂算子法对二维线性传输方程设计了满足两个守恒律的数值格式，与传统的Godunov型格式不同的是，我们的格式不仅能够保持数值解守恒，而且能够保持数值能量守恒。在这里，数值解不是被动的守恒，而是作为一个数值实体参与了格式的计算。

(a) t = 0时刻的解 (b) t = 0.2时刻的数值解

Figure 1. Results of example 4.2 at different time

图1. 格式计算算例4.2在不同时刻的结果

通过数值算例表明在远离极值点处，我们的格式具有误差相互抵消的超收敛性质，在数值模拟中能够很好的保持解的结构。