基于快速正弦离散变换的有限差分方法求解半线性抛物型方程

Finite Difference Method for Solving Semi Linear Parabolic Equations Based on Fast Discrete Sine Transform

摘要: 针对带有Dirichlet边界条件的二维半线性抛物方程给出二阶中心差分格式，利用Kronecker积写出二维拉普拉斯算子的微分矩阵。进而应用Crank-Nicolson方法进行时间离散，采用Picard迭代求解离散得到的非线性代数方程组。具体实现过程中结合快速离散正弦变换，本文方法的优点是减少了存储量并大幅度降低了计算时间。数值算例验证本文的方法可以更好地捕捉解的爆破现象。

Abstract: In this paper, the second-order finite difference scheme is applied for the numerical solution of the two-dimensional semi-linear parabolic equations with Dirichlet boundary conditions. We construct the differentiation matrix of two-dimensional Laplacian operator by Kronecker product. The time discretization method is chosen as Crank-Nicolson method. In every time level we solve the nonlinear algebraic equations by Picard iteration method. The fast discrete Sine transform is applied in the process of implementation. The major feature of the proposed method is that the memory requirement and CPU time are reduced obviously. Numerical examples show that the proposed method can better capture the blow up phenomenon of the solution.

1. 引言

半线性抛物与线性抛物方程的一个重要区别在于，随着时间发展其解具有奇异性。也就是说，尽管初始值是光滑的，但是在某个时间点T，方程的解也可能会趋向于无穷大，这种解的奇异性称之为爆破，T称之为爆破时间点。在本文中，我们考虑求解如下半线性抛物型方程：

$\{\begin{array}{l}{u}_t=u_{xx}+u_{yy}+f\left(u\right),\left(x,y\right)\in \Omega ,t>0\\ f\left(u\right)=u^p\end{array}$ (1)

边界和初始条件如下：

${u\left(x,y\right)|}_{\partial \Omega }=0,u\left(x,y,0\right)=u_0\left(x,y\right),\left(x,y\right)\in {\text{R}}^2$ (2)

其中， $f\left(u\right)=u^p$， $p>1$ 的常数， $\Omega =\left[a,b\right]\times \left[c,d\right]$ 为一个矩形区域， $\partial \Omega$ 代表边界， $u_0\left(x,y\right)$ 为初始值。关于方程(1)的爆破解有许多学者进行了研究 [1]，文献 [2] 给出了下面的定理：

定理1：令 $p_c\left(n\right)=1+\frac{2}{n}$，如果 $1<p<p_c$， $p_c$ 是问题的临界指数。对于(1)的任何非平凡解，存在一个爆破时间T，

$\underset{t\to T}{\lim }\left(\underset{x\in {ℝ}^n}{\sup }u\left(x,t\right)\right)=+\infty$

方程(1)在核反应动力学、奇异摄动问题等方面有广泛的应用 [3]。由于右端非线性项的存在，很难得到精确解，所以关于该问题数值解的研究具有非常重要的意义。近年来众多学者发展了丰富的数值方法求解该方程。F. De La Hoz等人提出了谱方法和指数时间差分法 [4]；王俊俊等人发展了混合有限元方法 [5]；张荣培等人将Chebyshev谱方法应用于求解此类方程 [6]；Ju等人发展了紧致指数积分因子法 [7]；熊之光等人提出了高效有限体积元法 [8]。但上述方法在求解过程中需要求解大型的非线性代数方程组，计算复杂，耗费时间较多。本文应用Kronecker积的性质将二维拉普拉斯算子的微分矩阵进行特征分解，结合快速正弦离散变换方法，节省了计算时间。

本文框架如下：第一部分给出半线性抛物方程的有限差分格式，应用Kronecker积写出二维拉普拉斯算子的微分矩阵K，时间离散上采用Crank-Nicolson方法。第二部分给出实现过程，应用Kronecker积的性质将微分矩阵K进行特征分解，结合快速正弦离散变换方法求解矩阵与向量的乘积，使得计算时间和内存大幅度减少。最后给出数值算例，验证本文给出的方法可以更加精确地捕捉爆破解。

2. 有限差分格式

首先用直线 $x=x_i$， $y=y_i$ 在 $\Omega$ 上打网格，其中网格节点为 $x_i=a+i\Delta x$， $y=c+j\Delta y$， $i=1,2,\cdots ,N_x$， $j=1,2,\cdots ,N_y$，网格步长为 $\Delta x=\frac{b-a}{N_x+1}$， $\Delta y=\frac{d-c}{N_y+1}$。设u在网格节点 $\left(x_i,y_j\right)$ 的数值解为 $u_{i,j}$，由Dirichlet边界条件，在网格内部点 $\left(x_i,y_j\right)$ 采用二阶中心差分格式：

$\frac{\text{d}{u}_{i,j}}{\text{d}t}+\frac{1}{\Delta x^2}\left(-u_{i-1,j}+2u_{i,j}-u_{i+1,j}\right)+\frac{1}{\Delta y^2}\left(-u_{i,j-1}+2u_{i,j}-u_{i,j+1}\right)=f\left(u_{ij}\right)$ (3)

定义离散解的矩阵U为：

$U={\left[\begin{array}{ccccc}{u}_{1,1}& u_{1,2}& \cdots & u_{1,N_{y-1}}& u_{1,N_y}\\ u_{2,1}& u_{2,2}& \cdots & u_{2,N_{y-1}}& u_{2,N}\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ u_{N_{x-1},1}& u_{N_{x-1},2}& \cdots & u_{N_{x-1},N_{y-1}}& u_{N_{x-1},N_y}\\ u_{N_x,1}& u_{N_x,2}& \cdots & u_{N_x,N_{y-1}}& u_{N_x,N_y}\end{array}\right]}_{N_x\times N_y}$

定义如下差分方程的微分矩阵为 $K_m$ ：

$K_m={\left[\begin{array}{cccccc}2& -1& 0& 0& \cdots & 0\\ -1& 2& -1& 0& \cdots & 0\\ 0& -1& 2& -1& \cdots & 0\\ ⋮& ⋮& \ddots & \ddots & \ddots & ⋮\\ 0& 0& \cdots & -1& 2& -1\\ 0& 0& \cdots & 0& -1& 2\end{array}\right]}_{N_m\times N_m}$

其中， $m=x,y$。

设 $I_m$ 分别为 $N_m$ 阶单位矩阵，将解矩阵向量化后得

$\mathcal{U}=vec\left(U\right)={\left(u_{1,1},\cdots ,u_{N_x,1},u_{1,2},\cdots ,u_{1,N_y},u_{1,N_y}\cdots ,u_{N_x,N_y}\right)}^{\text{T}}$ (4)

利用(3)式以及Kronecker积的定义，(1)式的二阶中心差分格式可以写成矩阵形式

$\frac{\text{d}\mathcal{U}}{\text{d}t}\text{+}\left(I_y\otimes \frac{1}{\Delta x^2}{K}_x+\frac{1}{\Delta y^2}{K}_y\otimes I_x\right)\mathcal{U}=\mathcal{F}\left(\mathcal{U}\right)$ (5)

其中 $\mathcal{F}\left(\mathcal{U}\right)=\mathcal{U}{.}^p$， ${.}^p$ 表示对 $\mathcal{U}$ 每个元素做运算。

计算终止时间为T，时间步长 $\tau =\frac{\Delta t}{N}$，则第n层时间步长 $t_n=n\tau ,n=0,1,2,\cdots ,N$。接下来应用Crank-Nicolson格式对(5)式进行求解，

$\frac{{\mathcal{U}}_{n+1}-{\mathcal{U}}_n}{\Delta t}-K\frac{{\mathcal{U}}_{n+1}+{\mathcal{U}}_n}{2}=\frac{{\mathcal{F}}_n+{\mathcal{F}}_{n+1}}{2}$ (6)

其中 $K=-\left(I_y\otimes \frac{1}{\Delta x^2}{K}_x+\frac{1}{\Delta y^2}{K}_y\otimes I_x\right)$，是拉普拉斯算子的微分矩阵。上式等价于

${\mathcal{U}}_{n+1}={\left(\frac{I}{\Delta t}-\frac{K}{2}\right)}^{-1}\left(\frac{I}{\Delta t}+\frac{K}{2}\right){\mathcal{U}}_n+{\left(\frac{I}{\Delta t}-\frac{K}{2}\right)}^{-1}\frac{\mathcal{F}\left({\mathcal{U}}_n\right)+{\mathcal{F}}_{\left(u_{n+1}\right)}}{2}$ (7)

下面采用Picard迭代求解非线性代数方程组

${\mathcal{U}}_{n+1}^{k+1}={\left(\frac{I}{\Delta t}-\frac{K}{2}\right)}^{-1}\left(\frac{I}{\Delta t}+\frac{K}{2}\right){\mathcal{U}}_n+{\left(\frac{I}{\Delta t}-\frac{K}{2}\right)}^{-1}\frac{{\mathcal{F}}_{\left(u_n\right)}+{\mathcal{F}}_{\left(u_{n+1}\right)}^k}{2}$ (8)

具体的Picard迭代算法如下：

$\begin{array}{l}{\mathcal{U}}_0={\mathcal{U}}_0\left(x_i,y_j\right)\\ \text{for}n=0,\text{1,}\cdots ,N\\ {\mathcal{U}}_{n+1}^0={\mathcal{U}}_n\\ \text{for}k=1,2,\cdots \\ \begin{array}{cc}& \end{array}{\mathcal{U}}_{n+1}^{k+1}={\left(\frac{I}{\Delta t}-\frac{K}{2}\right)}^{-1}\left(\frac{I}{\Delta t}+\frac{K}{2}\right){\mathcal{U}}_n+{\left(\frac{I}{\Delta t}-\frac{K}{2}\right)}^{-1}\cdot \frac{{\mathcal{F}}_{\left(u_n\right)}+{\mathcal{F}}_{\left(u_{n+1}\right)}^k}{2}\\ \begin{array}{cc}& \end{array}\text{while}{\Vert {\mathcal{U}}_{n+1}^{k+1}-{\mathcal{U}}_{n+1}^k\Vert }_{\infty }<\epsilon ={10}^{-13}\\ \begin{array}{c}\end{array}\text{continue}\\ \text{end}\\ {\mathcal{U}}_n={\mathcal{U}}_{n+1}^{k+1}\\ \text{end}\end{array}$

上述算法可分为内循环和外循环，外循环 $n=0,1,\cdots ,N$ 是关于时间的循环， $k=1,2,\cdots$ 为内循环，当 ${\Vert {\mathcal{U}}_{n+1}^{k+1}-{\mathcal{U}}_{n+1}^k\Vert }_{\infty }<\epsilon ={10}^{-13}$ 时，跳出内循环，外循环一直计算到时间终止。

3. 快速求解

在式(7)中需要求矩阵 ${\left(\frac{I}{\Delta t}-\frac{K}{2}\right)}^{-1}\left(\frac{I}{\Delta t}+\frac{K}{2}\right)$ 以及 ${\left(\frac{I}{\Delta t}-\frac{K}{2}\right)}^{-1}$ 与向量的乘积，两个矩阵乘积的运算相似，下面以 ${\left(\frac{I}{\Delta t}-\frac{\tilde{\mathcal{K}}}{2}\right)}^{-1}\left(\frac{I}{\Delta t}+\frac{\tilde{\mathcal{K}}}{2}\right)vec\left(U_n\right)$ 为例进行说明。为了快速计算，首先将矩阵K进行分解。下面给出两个定理。

定理2：矩阵 $K_x$ 进行对角化，定义对角矩阵 $\Lambda =diag\left({\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_{N_x}\right)$， ${\lambda }_j=2-2\cos \left(\frac{j\pi }{N_x+1}\right)$， $j=1,2,\cdots ,N_x$，则 $K_x{S}_x=S_x\Lambda$，即

$K_x=S_x\Lambda S_x^{-1}$ (9)

其中 ${\left(S_x\right)}_{jk}=\sin \left(\frac{jk\pi }{N_x+1}\right)$。

矩阵 $S_x$ 为离散正弦变换矩阵，是一个正交矩阵，即 $S_x=S_x^{-1}$，在Matlab中可以由命令dst(N)生成。该定理证明参考文献 [9]。对 $K_y$ 也可得类似的结果。

定理3：Kronecker积的性质 [9]：

1) $\left(A\otimes B\right)\left(C\otimes D\right)=AC\otimes BD$

2) ${\left(A\otimes B\right)}^{-1}=A^{-1}\otimes B^{-1}$

3) $\left(B^{\text{T}}\otimes A\right)Vec\left(X\right)=Vec\left(AXB\right)$

4) ${\left(A\otimes B\right)}^{\text{T}}=A^{\text{T}}\otimes B^{\text{T}}$

5) $diag\left(vec\left(A\right)\right)\cdot \left(vec\left(B\right)\right)=vec\left(A\cdot B\right)$

接下来利用定理2和3对矩阵K进行特征分解，由 $I_x=S_x^{\text{T}}{I}_x{S}_x$，利用定理3性质1和性质2可得：

$I_y\otimes K_x=\left(S_y\otimes S_x\right)\left(I_y\otimes {\Lambda }_x\right){\left(S_y\otimes S_x\right)}^{-1}$ (10)

对 $K_y\otimes I_x$ 可做类似运算。有下式成立

$K=\left(S_y\otimes S_x\right)\left(I_y\otimes \frac{1}{\Delta x^2}{\Lambda }_x+\frac{1}{\Delta y^2}{\Lambda }_y\otimes I_x\right){\left(S_y\otimes S_x\right)}^{-1}=\left(S_y\otimes S_x\right)\tilde{\Lambda }{\left(S_y\otimes S_x\right)}^{-1}$ (11)

这里 $\tilde{\Lambda }=I_y\otimes \frac{1}{\Delta x^2}{\Lambda }_x+\frac{1}{\Delta y^2}{\Lambda }_y\otimes I_x$ 是一个对角矩阵。

定义矩阵

$\Lambda ={\left[\begin{array}{cccc}\frac{{\lambda }_1^x}{\Delta x^2}+\frac{{\lambda }_1^y}{\Delta y^2}& \frac{{\lambda }_1^x}{\Delta x^2}+\frac{{\lambda }_2^y}{\Delta y^2}& \cdots & \frac{{\lambda }_1^x}{\Delta x^2}+\frac{{\lambda }_{N_y}^y}{\Delta y^2}\\ \frac{{\lambda }_2^x}{\Delta x^2}+\frac{{\lambda }_1^y}{\Delta y^2}& \frac{{\lambda }_2^x}{\Delta x^2}+\frac{{\lambda }_2^y}{\Delta y^2}& \cdots & \frac{{\lambda }_2^x}{\Delta x^2}+\frac{{\lambda }_{N_y}^y}{\Delta y^2}\\ ⋮& ⋮& \ddots & ⋮\\ \frac{{\lambda }_{N_x}^x}{\Delta x^2}+\frac{{\lambda }_1^y}{\Delta y^2}& \frac{{\lambda }_{N_x}^x}{\Delta x^2}+\frac{{\lambda }_2^y}{\Delta y^2}& \cdots & \frac{{\lambda }_{N_x}^x}{\Delta x^2}+\frac{{\lambda }_{N_y}^y}{\Delta y^2}\end{array}\right]}_{N_x\times N_y}$

则易证明

$\tilde{\Lambda }=diag\left(vec\left(\Lambda \right)\right)$ (12)

根据矩阵K的分解，可得下面的分解

${\left(\frac{I}{\Delta t}-\frac{\mathcal{K}}{2}\right)}^{-1}\left(\frac{I}{\Delta t}+\frac{\mathcal{K}}{2}\right)=\left(S_y\otimes S_x\right){\left(\frac{I}{\Delta t}-\frac{\tilde{\Lambda }}{2}\right)}^{-1}\left(\frac{I}{\Delta t}+\frac{\tilde{\Lambda }}{2}\right){\left(S_y\otimes S_x\right)}^{-1}$ (13)

基于上述矩阵可分解性质，下面分三步利用利用快速正弦变换完成矩阵与向量乘积的运算实现 ${\left(\frac{I}{\Delta t}-\frac{\tilde{\mathcal{K}}}{2}\right)}^{-1}\left(\frac{I}{\Delta t}+\frac{\tilde{\mathcal{K}}}{2}\right)vec\left(U_n\right)$。

Step1

$\left(S_y\otimes S_x\right)vec\left(U_n\right)=vec\left(S_x{U}_n{S}_y\right)$ (14)

Step2

定义 $L=\left(l_{ij}\right)$，其中 $l_{ij}=\left(\frac{1}{\Delta t}+\frac{l_i^x+l_j^y}{2}\right)/\left(\frac{1}{\Delta t}-\frac{l_i^x+l_j^y}{2}\right)$，则由(12)可得 ${\left(\frac{I}{\Delta t}-\frac{\tilde{\Lambda }}{2}\right)}^{-1}\left(\frac{I}{\Delta t}+\frac{\tilde{\Lambda }}{2}\right)=diag\left(vec\left(L\right)\right)$，从而

${\left(\frac{I}{\Delta t}-\frac{\tilde{\Lambda }}{2}\right)}^{-1}\left(\frac{I}{\Delta t}+\frac{\tilde{\Lambda }}{2}\right)vec\left(S_x{U}_n{S}_y\right)=vec\left(S_x{U}_n{S}_y.\ast L\right)$ (15)

其中 $.\ast$ 为矩阵中每个元素相乘。

Step3

$\left(S_y^{-1}\otimes S_x^{-1}\right)vec\left(S_x{U}_n{S}_y.\ast L\right)=vec\left(S_x^{-1}\left(S_x{U}_n{S}_y.\ast L\right)S_y^{-1}\right)$ (16)

注：Step1中，在Matlab中利用 $U_1=idst\left(idst{\left(U_n^{\text{T}}\right)}^{\text{T}}\right)$ 实现；Step3中，在Matlab中利用 $U=dst\left(dst{\left(U_2^{\text{T}}\right)}^{\text{T}}\right)$ 实现。

4. 数值实验

应用本文提出的有限差分法求解下面的数值算例。考虑二维计算区域 $\Omega =\left[-10,10\right]$，针对具有齐次狄利克雷边界条件的方程(1)进行求解。网格剖分128 × 128，时间离散步长时间步长 $\Delta t={10}^{-3}$，因为临界指数 $p_c=2$，所以在数值实验中选择 $p=1.5$，对不同的时间点t进行数值模拟。

情形1. 初值选取为 $u_0\left(x,y\right)=10\exp \left(-10\left(x^2+y^2\right)\right)$。三个时间点 $t=2$， $t=6.5$， $t=7.01$ 的计算结果在图1中列出。初始条件在原点峰值为10，t取如图1中的各值时，可以看到随着时间的进行，单峰图像逐渐变陡，向原点不断聚拢，最终在 $t=7.01$ 达到爆破的峰值 $u=15000$。

Figure 1. The numerical simulation results of case 1 when $t=2$, $t=6.5$, $t=7.01$

图1. 情形1在计算时间 $t=2$， $t=6.5$， $t=7.01$ 下的数值模拟结果

情形2. 初值选取 $u_0\left(x,y\right)=10\exp \left(-10\left({\left(x+1\right)}^2+{\left(y+1\right)}^2\right)\right)+10\exp \left(-10\left({\left(x-1\right)}^2+{\left(y-1\right)}^2\right)\right)$，其他参数与情形1相同，三个时间点 $t=0.3$， $t=1$， $t=4.89$ 的计算结果如图2所示。初始条件在原点的峰值为10，从图像可以看出它们是关于原点对称的，随着时间的进行，双峰逐渐聚合，在 $t=1$ 时效果较为明显，在 $t=4.89$ 时爆破值达到4 × 10⁴。

Figure 2. The numerical simulation results of case 2 when $t=0.3$, $t=1$, $t=4.89$

图2. 情形2在计算时间 $t=0.3$， $t=1$， $t=4.89$ 下的数值模拟结果

情形3. 初值选取

$\begin{array}{l}{u}_0\left(x,y\right)=10\exp \left(-10\left({\left(x+5\right)}^2+{\left(y+5\right)}^2\right)\right)+10\exp \left(-10\left({\left(x-2\right)}^2+{\left(y+2\right)}^2\right)\right)\\ \begin{array}{cc}\begin{array}{cc}\begin{array}{cc}& \end{array}& \end{array}& +10\exp \left(-10\left({\left(x+2\right)}^2+{\left(y-2\right)}^2\right)\right)+10\exp \left(-10\left({\left(x-5\right)}^2+{\left(y-5\right)}^2\right)\right)\end{array}\end{array}$

其他参数与情形2相同，三个时间点 $t=0.\text{2}$， $t=\text{4}\text{.4}$， $t=\text{5}\text{.96}$ 的计算结果如图3所示。初始条件的峰值为10。可以发现随着时间的进行，四个峰的聚拢方向依然指向原点，最后实现爆破。这三个实验的结果非常相似，无论初值有多少极值点，在一定的时间过程后，最终都会向原点聚拢，进而达到爆破的结果。数值模拟结果验证了定理1的理论。

Figure 3. The numerical simulation results of case 3 when $t=2$, $t=4.4$, $t=5.96$

图3. 情形3在计算时间 $t=2$， $t=4.4$， $t=5.96$ 下的数值模拟结果

为了显示本文方法计算速度快的优点，表1给出使用快速离散正弦变换和不使用快速离散正弦变换计算式(6)的时间比较。从表1可以发现，不使用快速离散正弦变换的计算时间较长，而应用本文的数值方法进行模拟的运行时间大大减少。

Table 1. Time comparison of three cases of DST and No DST

表1. 三种情形下使用速离散正弦变换和不使用速离散正弦变换的时间比较

5. 总结

本文应用快速离散正弦变换求解半线性抛物线型方程，该方法能够快速地求解离散后得到的非线性代数方程组，提高了计算效率。通过求解齐次狄利克雷边界条件下的半线性抛物线型方程，可以发现随时间推移数值结果形状的变化规律，即在一定时间后，数值解发生了爆破。在本文问题求解的过程中，当 $N=N_x=N_y$ 时，当利用高斯消去法解代数方程组时，计算复杂度高达 ${\rm O}\left(N^4\right)$，但通过本文所应用的快速离散正弦变换方法，复杂度仅有 ${\rm O}\left(N^2{\log }_{2}N\right)$，大大降低了算法的复杂度，使得程序运算时间明显缩短。下一步的工作是将该方法推广至求解四阶的抛物线型方程。

NOTES

^*通讯作者。