半线性抛物问题高效有限体积元法

The Efficient Finite Volume Element Methods for Semilinear Parabolic Equations

摘要: 本文研究了一类半线性抛物方程的有限体积元全离散格式。首先在空间上以插值系数线性有限体积元进行半离散，得到关于时间的一阶非线性常微分方程组的初值问题，然后在时间上采用向后差分方法得到全离散格式。其次讨论了该全离散格式的稳定性和收敛性。最后给出了一个数值例子说明所研究方法的高效性。

Abstract: In this paper, the fully discrete finite volume element schemes for a class of semi-linear parabolic equations are studied. Firstly, the initial value problem of the first-order nonlinear ordinary dif-ferential equations with time is obtained by interpolation coefficient finite volume element in space, and then the fully discrete scheme is obtained by backward difference method in time. Secondly, the stability and convergence of the fully discrete scheme are discussed. Finally, a numerical example is given to illustrate the efficiency of the proposed method.

1. 前言

有限体积元方法(Finite Volume Method)是一种求解偏微分方程的数值方法，国内学者最早由李荣华教授提出 [1]，在国外，该方面文章较早出现于二十世纪八十年代末 [2]。插值系数有限元方法的基本思想最早在1980年由Zlamal等提出 [3] [4]，并用这种思想成功解决拟线性抛物问题，它的基本思路是在用有限元方法求解的过程中对非线性项直接插值处理，这样既简化了计算存储空间，又能保证所求数值解的超收敛性 [5]。插值系数有限元方法对于求解半线线性问题非常有效 [6] [7] [8]。类似地，将插值思想应用于微分方程的有限体积元法中，就得到插值系数有限体积元方法 [9] [10] [11] [12]。

本文将把插值系数有限体积元方法推广到半线性抛物问题中。考虑如下半线性抛物方程的初边值问题的有限体积元全离散格式

$\{\begin{array}{l}{u}_t-{\left(p\left(x\right)u_x\right)}_x+f\left(u\right)=g\left(x,t\right),\left(x,t\right)\in \left(a,b\right)\times \left[0,T\right],\\ u\left(a,t\right)=u\left(b,t\right)=0,t\in \left[0,T\right],\\ u\left(x,0\right)=\phi \left(x\right),x\in \left[a,b\right],\end{array}$ (1.1)

其中 $p\left(x\right),f\left(s\right),g\left(x,t\right)$ 为充分光滑函数，且 $f^{\prime }\left(s\right)>0$。本文将作如下安排：第2节构造出有限体积元的全离散格式，先在空间方向上进行半离散，再在时间方向上分别采用向前差分法、向后差分法得到两种全离散格式。第3节针对本文所给的问题得到的两种全离散格式，讨论其稳定性并给出证明。第4节给出两个数值例子以说明该方法的有效性。

2. 有限体积元格式

2.1. 空间有限体积元半离散

问题(1.1)的变分形式为，对每一个t，寻求 $u\left(x,t\right)\in H_0^1\left(I\right)$，使得

$\{\begin{array}{l}\left(u_t,v\right)+a\left(u,v\right)+\left(f\left(u\right),v\right)-\left(g,v\right)=0,\forall v\in H_0^1\left(I\right),\\ u\left(x,0\right)=\phi \left(x\right),\end{array}$ (2.1)

其中 $H_0^1\left(I\right)=\left\{v\left(x,t\right)|v\left(x,t\right)\in H^1\left(I\right),I=\left[a,b\right],v\left(a,t\right)=v\left(b,t\right)=0\right\},t\in \left[0,T\right]$，

$a\left(u,v\right)={\displaystyle {\int }_a^{b}p{u}_x{v}_x\text{d}x}.$ (2.2)

对区间 $I=\left[a,b\right]$ 作如下剖分：

设单元 $I_j=\left[x_{j-1},x_j\right]$，每个单元的中点记为 $x_{j-1/2}=\left(x_j+x_{j-1}\right)/2$，步长 $$， $h=\underset{1\le j\le n}{\max }{h}_j$，

并假设剖分 $J_h$ 是拟一致的，即存在常数 $C>0$ 使得 $h\le C{h}_j$，另外设 $I_0^{\ast }=\left[x_0,x_{1/2}\right]$， $I_i^{\ast }=\left[x_{i-1/2},x_{i+1/2}\right]$ $\left(i=1,2,\cdots ,n-1\right)$， $I_n^{\ast }=\left[x_{n-1/2},x_n\right]$，那么所有的 $I_i^{\ast }\left(i=0,1,\cdots ,n\right)$ 构成 $J_h$ 的对偶剖分 $J_h^{*}$，每个 $I_i^{\ast }\left(i=0,1,\cdots ,n\right)$ 也称为体积控制元.设基于剖分 $J_h$ 的试探函数空间为 $U_h=\left\{u|u\in C\left(I\right),u在每个单元为线性多项式\right\},$ 对每个t， $u\left(x,t\right)$ 在节点 $x_j$ 上的基函数为：

${\varphi }_j\left(x\right)=\{\begin{array}{l}\frac{x-x_{j-1}}{h_j},x\in I_j,\\ \frac{x_{j+1}-x}{h_{j+1}},x\in I_{j+1},\\ 0,其他,\end{array}j=1,2,\cdots ,n-1,$

设基于对偶剖分 $J_h^{*}$ 的函数空间 $V_h$ 的为分片常数空间，对每个t， $u\left(x,t\right)$ 在节点 $x_j$ 上的基函数为：

${\phi }_j\left(x\right)=\{\begin{array}{l}1,x\in I_j^{\ast },\\ 0,其他,\end{array}j=1,2,\cdots ,n-1,$

对每个t，我们再定义插值算子 ${\text{I}}_h:C\left(I\right)\to U_h$， $v\left(x,t\right)\mapsto {\displaystyle \underset{j=1}{\overset{n-1}{\sum }}v\left(x_j,t\right){\varphi }_j\left(x\right)}$，和插值算子 ${\text{I}}_h^{*}:C\left(I\right)\to V_h$， $v\left(x,t\right)\mapsto {\displaystyle \underset{j=1}{\overset{n-1}{\sum }}v\left(x_j,t\right){\phi }_j\left(x\right)}$。则问题(1.1)的插值系数线性有限体积元解 $u_h\left(x,t\right)\in U_h$ 满足

$\begin{array}{l}\left(u_{h,t}\left(x,t\right),I_h^{*}{v}_h\right)+a\left(u_h\left(x,t\right),I_h^{*}{v}_h\right)+\left(I_{h}f\left(u_h\left(x,t\right)\right),I_h^{*}{v}_h\right)-\left(g\left(x,t\right),I_h^{*}{v}_h\right)=0,\\ \forall v_h\in U_h,\end{array}$ (2.3)

这里 $u_h\left(x,t\right),{\text{I}}_{h}f\left(u_h\left(x,t\right)\right)$ 分别为

$u_h\left(x,t\right)={\displaystyle \underset{j=1}{\overset{n-1}{\sum }}{u}_h\left(x_j,t\right){\varphi }_j\left(x\right)}={\displaystyle \underset{j=1}{\overset{n-1}{\sum }}{u}_j\left(t\right){\varphi }_j\left(x\right)},$

${\text{I}}_{h}f\left(u_h\left(x,t\right)\right)={\displaystyle \underset{j=1}{\overset{n-1}{\sum }}f\left(u_h\left(x_j,t\right)\right){\varphi }_j\left(x\right)}={\displaystyle \underset{j=1}{\overset{n-1}{\sum }}f\left(u_j\left(t\right)\right){\varphi }_j\left(x\right)},$

其中 $u_j\left(t\right)=u_h\left(x_j,t\right),f\left(u_j\left(t\right)\right)=f\left(u_h\left(x_j,t\right)\right)$。则式(2.3)可以等价的表示为

$\begin{array}{l}\left(u_{h,t}\left(x,t\right),{\phi }_j\right)+a\left(u_h\left(x,t\right),{\phi }_j\right)+\left({\text{I}}_{h}f\left(u_h\left(x,t\right)\right),{\phi }_j\right)-\left(g\left(x,t\right),{\phi }_j\right)=0,\\ j=1,2,\cdots ,n-1,\end{array}$ (2.4)

即

$\begin{array}{l}{\displaystyle {\int }_{x_{j-1/2}}^{x_{j+1/2}}{u}_{h,t}\text{d}x}={\displaystyle {\int }_{x_{j-1/2}}^{x_{j+1/2}}{\left(p{u}_{h,x}\right)}_x\text{d}x}+{\displaystyle {\int }_{x_{j-1/2}}^{x_{j+1/2}}{\text{I}}_{h}f\left(u_h\right)\text{d}x}+{\displaystyle {\int }_{x_{j-1/2}}^{x_{j+1/2}}g\left(x,t\right)\text{d}x},\\ j=1,2,\cdots ,n-1,\end{array}$ (2.5)

其中 $u_{h,t}=\frac{\partial u_h\left(x,t\right)}{\partial t},u_{h,x}=\frac{\partial u_h\left(x,t\right)}{\partial x}$。注意到 $u_0\left(t\right)=u_h\left(a,t\right)=0$，且 $u_n\left(t\right)=u_h\left(b,t\right)=0$，(2.5)得到如下计

算格式：

$\{\begin{array}{l}\frac{3}{8}\left(h_1+h_2\right){u^{\prime }}_1\left(t\right)+\frac{1}{8}{h}_2{u^{\prime }}_2\left(t\right)+\left(\frac{p_{1/2}}{h_1}+\frac{p_{3/2}}{h_2}\right)u_1\left(t\right)-\frac{p_{3/2}}{h_2}{u}_2\left(t\right)\\ +\frac{3}{8}\left(h_1+h_2\right)f\left(u_1\left(t\right)\right)+\frac{1}{8}{h}_{2}f\left(u_2\left(t\right)\right)=g_1\left(t\right)+\frac{1}{8}{h}_1{f}_0,\\ \frac{1}{8}{h}_j{u^{\prime }}_{j-1}\left(t\right)+\frac{3}{8}\left(h_j+h_{j+1}\right){u^{\prime }}_j\left(t\right)+\frac{1}{8}{h}_j{u^{\prime }}_{j+1}\left(t\right)-\frac{p_{j-1/2}}{h_j}{u}_{j-1}\left(t\right)+\left(\frac{p_{j-1/2}}{h_j}+\frac{p_{j+1/2}}{h_{j+1}}\right)u_j\left(t\right)\\ -\frac{p_{j+1/2}}{h_{j+1}}{u}_{j+1}\left(t\right)+\frac{1}{8}{h}_{j}f\left(u_{j-1}\left(t\right)\right)+\frac{3}{8}\left(h_j+h_{j+1}\right)f\left(u_j\left(t\right)\right)+\frac{1}{8}{h}_{j}f\left(u_{j+1}\left(t\right)\right)=g_j\left(t\right),j=2,\cdots ,n-2,\\ \frac{1}{8}{h}_{n-1}{u^{\prime }}_{n-2}\left(t\right)+\frac{3}{8}\left(h_{n-1}+h_n\right){u^{\prime }}_{n-1}\left(t\right)-\frac{p_{n-3/2}}{h_{n-1}}{u}_{n-2}\left(t\right)+\left(\frac{p_{n-3/2}}{h_{n-1}}+\frac{p_{n-1/2}}{h_n}\right)u_{n-1}\left(t\right)\\ -\frac{1}{8}{h}_{n-1}f\left(u_{n-2}\left(t\right)\right)+\frac{3}{8}\left(h_{n-1}+h_n\right)f\left(u_{n-1}\left(t\right)\right)=g_{n-1}\left(t\right)+\frac{1}{8}{h}_n{f}_n,\end{array}$ (2.11)

其中

$f_0=f\left(u_h\left(x_0,t\right)\right)=f\left(u\left(a,t\right)\right)=f\left(0\right),$

$f_n=f\left(u_h\left(x_n,t\right)\right)=f\left(u\left(b,t\right)\right)=f\left(0\right).$

$g_j\left(t\right)={\displaystyle {\int }_{x_{j-1/2}}^{x_{j+1/2}}g\left(x,t\right)\text{d}x},j=1,2,\cdots ,n-1,$

如下向量记号：

$u\left(t\right)={\left(\begin{array}{ccccc}{u}_1\left(t\right)& u_2\left(t\right)& u_3\left(t\right)& \cdots & u_{n-1}\left(t\right)\end{array}\right)}^{\text{T}},$

$\stackrel{\dot{}}{u}\left(t\right)={\left(\begin{array}{ccccc}{u^{\prime }}_1\left(t\right)& {u^{\prime }}_2\left(t\right)& {u^{\prime }}_3\left(t\right)& \cdots & {u^{\prime }}_{n-1}\left(t\right)\end{array}\right)}^{\text{T}},$

$F\left(u\right)={\left(\begin{array}{ccccc}f\left(u_1\left(t\right)\right)& f\left(u_2\left(t\right)\right)& f\left(u_3\left(t\right)\right)& \cdots & f\left(u_{n-1}\left(t\right)\right)\end{array}\right)}^{\text{T}},$

$G={\left(\begin{array}{ccccc}{g}_1\left(t\right)+\frac{1}{8}{h}_{1}f\left(0\right)& g_2\left(t\right)& g_3\left(t\right)& \cdots & g_{n-1}\left(t\right)+\frac{1}{8}{h}_{n}f\left(0\right)\end{array}\right)}^{\text{T}}$

相应的 $\stackrel{\dot{}}{u}\left(t\right),u\left(t\right)$ 系数矩阵记为 $A,B$ ,则式(2.11)用矩阵表示为：

$A\stackrel{\dot{}}{u}\left(t\right)+Bu\left(t\right)+AF\left(u\right)-G\left(t\right)=0.$ (2.12)

另记

则方程(1.1)在空间上用线性有限体积元方法离散，得到如下的 $$ 关于t的一阶常微分方程组

$\{\begin{array}{l}A\stackrel{\dot{}}{u}\left(t\right)+Bu\left(t\right)+AF\left(u\left(t\right)\right)-G\left(t\right)=0,\\ u\left(0\right)={\alpha }_0.\end{array}$ (2.13)

2.2. 时间差分全离散

在式(2.13)中，由于 $A$ 是严格对角占优矩阵，故 $A$ 可逆，则式(2.13)可化为

$$ (2.14)

将时间区域 $J=\left[0,T\right]$ 分成N等分，节点为 $t_k=t_0+i\Delta t$，步长为 $\Delta t=T/N$，将式(2.14)在节点处离散化，对充分光滑的解 $u\left(t\right)$，利用向后差分可得

$\{\begin{array}{l}\frac{u\left(t_k\right)-u\left(t_{k-1}\right)}{\Delta t}+A^{-1}Bu\left(t_k\right)+F\left(u\left(t_k\right)\right)-A^{-1}G\left(t_k\right)=0,\\ u\left(0\right)={\alpha }_0.\end{array}$ (2.17)

即得到单步隐式格式：

$\{\begin{array}{l}M{u}^{\left(k\right)}+\Delta tF\left(u^{\left(k\right)}\right)-\Delta t{A}^{-1}G=u^{\left(k-l\right)},\\ u^{\left(0\right)}={\alpha }_0.\end{array}$ (2.18)

其中 $M=I+\Delta t{A}^{-1}B$，， $k=0,1,2,\cdots ,N$。在式(2.18)中，令

$H\left(u^{\left(k\right)}\right)=M{u}^{\left(k\right)}+\Delta tF\left(u^{\left(k\right)}\right)-\Delta t{A}^{-1}{G}^{\left(k\right)}-u^{\left(k-l\right)},$ (2.19)

则式(2.18)可化为

(2.20)

利用Newton迭代法可得

$u^{\left(k\right)，\left(l+1\right)}=u^{\left(k\right),\left(l\right)}-D{H}^{-1}\left(u^{\left(k\right),\left(l\right)}\right)H\left(u^{\left(k\right),\left(l\right)}\right)，k=1,2,3,\cdots ,N,$ (2.21)

这里的 $u^{\left(k\right),\left(l\right)}$ 表示 $$ 的第l次迭代， $DH\left(u^{\left(k\right),\left(l\right)}\right)$ 表示对矩阵求导， $D{H}^{-1}\left(u^{\left(k\right),\left(l\right)}\right)$ 表示的逆矩阵，对式(2.19)求导可得

$DH\left(u\right)=M+\Delta t{F}^{\prime }\left(u\right)，$ (2.22)

其中

$F^{\prime }\left(u\right)=\left[\begin{array}{cccc}{f}^{\prime }\left(u_1\right)& 0& \cdots & 0\\ 0& f^{\prime }\left(u_2\right)& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & f^{\prime }\left(u_{n-1}\right)\end{array}\right].$

利用(2.21)式，就可得到方程(1.1)的数值解。

3. 有限体积元全离散格式的稳定性与收敛性

本节讨论有限体积元全离散格式的稳定性，在这之前，先做一些准备工作。

定义3.1 设 $A\in R^{n\times n}$ 的特征值为 ${\lambda }_i\left(i=1,2,\cdots ,n\right)$，称

$\rho \left(A\right)=\underset{1\le i\le n}{\max }\left|{\lambda }_i\right|$

为 $A$ 的谱半径。

设 $x^{*}$ 是方程组 $x=Bx+f$ 的精确解， $x^{\left(k\right)}$ 是其迭代解，记误差向量 $e^{\left(k\right)}=x^{\left(k\right)}-x^{*}$，则有

即

$e^{\left(k\right)}=B^k{e}^{\left(0\right)},k=0,1,\cdots ,n-1,$

其中与k无关，所以 $x=Bx+f$ 迭代收敛的充要条件为对于任意 $x^{\left(0\right)}\in R^n$，都有

$\underset{k\to \infty }{\lim }{e}^{\left(k\right)}=\underset{k\to \infty }{\lim }{B}^k{e}^{\left(0\right)}=0.$

根据文献 [13]，有如下两个引理：

引理3.1 设矩阵 $B\in R^{n\times n}$，则的充要条件是 $\rho \left(B\right)<1$。

引理3.2 对于任意的初始向量 $x^{\left(0\right)}$ 和右端向量f，解方程组 $x=Bx+f$ 的迭代法 $x^{\left(k+1\right)}=B{x}^{\left(k\right)}+f$ 稳定的充要条件是 $\rho \left(B\right)<1$。

引理3.3 对于本文中的矩阵 $A$，若 $A=O\left(h\right)$，则 $A^{-1}=O\left(h^{-1}\right)$。

证明 由于 $A=O\left(h\right)$，则

$\left|A\right|=O\left(h^{n-1}\right),\left|A_{ij}\right|=O\left(h^{n-2}\right),$

其中 $A_{ij}$ 表示伴随矩阵的元素，从而

$A^{-1}=\frac{A^{*}}{\left|A\right|}=\frac{O\left(h^{n-2}\right)}{O\left(h^{n-1}\right)}=O\left(h^{-1}\right).$

引理3.3得证。

为了方便起见，引入一些记号，令 $u_j^n$ 为方程(1.1)的精确解， $u_j\left(t_n\right)$ 为方程(1.1)的半离散解， $u\left(x_j,t_n\right)$ 为方程(1.1)的全离散解，即有限体积元方法的数值解。

有了上面的准备，下面来讨论有限体积元方法的稳定性及其证明。

本小节中考虑线性有限体积元向后差分全离散格式，由式(1.19)得向后差分全离散迭代格式为

$\{\begin{array}{l}M{u}^{\left(k\right)}=u^{\left(k-1\right)},\\ u^{\left(0\right)}={\alpha }_0,\end{array}$ (3.5)

其中 $M=I+\Delta t{A}^{-1}B$ .

定理3.2对方程(1.1)采用有限体积元向后差分全离散迭代格式(3.5)时，在任何剖分的情况下都是稳定的，即稳定性与时间剖分、空间剖分的密度无关。

证明 由于 $A^{-1}B$ 是正定的，因此特征值 $$，从而 $M=I+\Delta t{A}^{-1}B$ 的特征值 $\lambda \left(M\right)>1$，即 $\rho \left(M\right)>1$，亦即 $\Vert M\Vert >1$，所以 $\Vert {\left(M^{-1}\right)}^k\Vert <1$，即

$\rho \left({\left(M^{-1}\right)}^k\right)<1,$ (3.6)

由(3.5)可得 $u^{\left(k\right)}=M^{-1}{u}^{\left(k-1\right)}={\left(M^{-1}\right)}^k{u}^{\left(0\right)}$，所以

${\epsilon }^{\left(k\right)}={\left(M^{-1}\right)}^k{\epsilon }^{\left(0\right)},$

其中 ${\epsilon }^{\left(k\right)}=u_j^n-u\left(x_j,t_n\right)$，，则由式(3.6)及引理3.2可得迭代格式(3.5)是稳定的，与h无关。定理3.2得证。

4. 数值例子

为验证半线性抛物问题的有限体积元方法的有效性，本节给出下面例子。

例 求解如下半线性抛物问题

$\{\begin{array}{l}{u}_t-u_{xx}+u^3=g\left(x,t\right),\left(x,t\right)\in \left(0,1\right)\times \left[0,T\right],\\ u\left(0,t\right)=u\left(1,t\right)=0,t\in \left[0,T\right],\\ u\left(x,0\right)=x\left(1-x\right),x\in \left[0,1\right].\end{array}$ (4.1)

这里取 $g\left(x,t\right)=\left(x^2-x+2\right){\text{e}}^{-t}+x^3{\left(1-x\right)}^3{\text{e}}^{-3t}$，则其精确解为 $u\left(x,t\right)=x\left(1-x\right){\text{e}}^{-t}$。

分别采用线性有限体积元半离散，时间向后差分全离散格式进行数值计算。依次取均匀空间步长h和均匀时间步长 $\Delta t$。表1给出了当 $t=1$ 时，网格比为 $\Delta t/h^3=1/10$ 的节点在N = 10，20，40，80的绝对误差值。其中u表示精确解， $E_{H1}$ 至表示在N = 10，20，40，80时相应剖分下的绝对误差。表2给出了相应剖分下的绝对误差最大值和节点平均值误差，其中 $e_{H1}$ 表示在不同剖分下的节点绝对值误差最大值， $e_{H2}$ 表示节点平均误差值。

Table 1. Absolute errors of linear finite element with time backward difference for (4.1) under different partitions

表1. 问题(4.1)线性有限元时间向后差分在不同剖分下的绝对误差

Table 2. Problem (4.1) maximum absolute error and average error of linear finite element with time backward difference under different subdivisions

表2. 问题(4.1)线性有限元时间向后差分在不同剖分下的绝对误差最大值和平均误差值

从表1，表2可以看出用线性有限元时间向后差分全离散格式求解半线性抛物问题是可靠的。它的精度不仅与剖分密度有关，还与时间和空间的相对密度有关，也就是说，当剖分密度越大，时间和空间的相对密度越大，则利用该方法的精度越高。

基金项目

本文为国家自然科学基金项目(11571102)资助课题。