利用Laplace变换求解时间分数阶Swift-Hohenberg方程

Solving Time Fractional Swift-Hohenberg Equation Using Laplace Transform

摘要: 针对时间分数阶Swift-Hohenberg方程，本文提出了基于Laplace变换的高效数值算法。首先利用Laplace变换将原Caputo型分数阶方程转化为整数阶方程，然后利用算子分裂法进一步将其分解成线性方程和非线性方程，其中，非线性方程通过积分法近似求解，线性方程通过Crank-Nicolson差分格式求解，最后通过数值实验验证了所给格式的有效性。

Abstract: For the time fractional Swift-Hohenberg equation, this paper proposes an efficient numerical algorithm based on Laplace transform. First, the Laplace transform is used to transfrom the original Caputo fractional equation into an integer-order equation, and then the operator splitting method is used to further decompose the equation into linear equation and nonlinear equation. The nonlinear equation is approximately solved by the integral method, and the linear equation is solved by the Crank-Nicolson scheme and central difference. Finally, the validity of the given scheme is verified through numerical experiments.

1. 引言

Swift-Hohenberg方程是1977年Swift和Hohenberg通过Rayleigh-Bénard不稳定球面波的模型 [1] 给出的，近年来相关研究涉及激光器 [2]、平面几何 [3]、时滞动力学 [4]、图灵分叉 [5] 等多方面，在物理学领域发挥着重要的作用，特别在流体力学和激光的对流问题中有着广泛应用。

关于Swift-Hohenberg方程的求解，中外研究者们已经提出了很多可行的解析方法和数值方法，例如指数法 [6]、同伦摄动法 [7] (Homotopy perturbation Method)、微分变换法 [7] (Differential Transform Method)、分数变分迭代法 [8] (Fractional Variational Iteration Method)、同伦分析法 [9] (Homotopy Analysis Method)等。

分数阶微积分诞生后的很长时间内，研究都局限在理论数学领域内，直到二十世纪七十年代，研究人员发现自然科学与工程领域的分形几何、记忆过程等相关现象与分数阶微积分存在密切的联系，才掀起了研究热潮。整数阶导数往往局限于某个物理过程在特定时刻的变化和性质，而分数阶导数模型既克服了经典整数阶微分理论模型实验结果吻合度低的缺点，又能较好地体现系统函数发展的历史依赖过程，为我们研究复杂物理力学问题提供了很大的帮助。

Riemann-Liouville型和Caputo型是目前最常用的两种分数阶导数的定义方式 [10]。观察发现，时间分数阶导数的计算常使用Caputo定义，而空间分数阶导数的计算更多使用Riemann-Liouville定义和级数定义，当然分数阶导数的不同定义也可以理解为不同的方法分别求解分数阶导数。

本文研究如下的时间分数阶Swift-Hohenberg方程：

$\{\begin{array}{l}{}_0{}^{c}D{}_t^{\alpha }u\left(x,t\right)=-\left[{\left(1+\frac{{\partial }^2}{\partial x^2}\right)}^{2}u+f\left(u\right)\right],a\le x\le b,0\le t\le T,0<\alpha <1\\ u\left(a,t\right)=0,u\left(b,t\right)=0,0\le t\le T\\ u\left(x,0\right)=u_0\left(x\right),a\le x\le b\end{array}$ (1)

其中， ${}_0{}^{c}D{}_t^{\alpha }u\left(x,t\right)$ 是Caputo分数阶导数， $f\left(u\right)$ 是函数 $F\left(u\right)$ 关于u的导，参数 $\eta >0$， $\epsilon \ge 0$ 都是具有物理意义的常量。

${}_0{}^{c}D{}_t^{\alpha }u\left(x,t\right)=\frac{1}{\Gamma \left(1-\alpha \right)}{\displaystyle {\int }_0^t\frac{u_{\tau }\left(x,t\right)}{{\left(t-\tau \right)}^{\alpha }}\text{d}\tau },0<\alpha <1$ (2)

$F\left(u\right)=\frac{1}{4}{u}^4-\frac{\eta }{3}{u}^3-\frac{\epsilon }{2}{u}^2$ (3)

$f\left(u\right)=F^{\prime }\left(u\right)=u^3-\eta u^2-\epsilon u$ (4)

首先，利用Laplace变换逼近时间分数阶Swift-Hohenberg方程的分数阶项，将其转化为整数阶问题，然后，利用算子分裂法将得到的整数阶方程分解为线性部分和非线性部分，将解算子分别记为 $S^1$ 和 $S^2$，最后通过以下格式求解，得到方程的数值解。其中，非线性部分利用左矩形法近似求解，线性部分利用Crank-Nicolson差分格式求解。

$u\left(t+\Delta t\right)=S^1\left(\frac{\Delta t}{2}\right)S^2\left(\Delta t\right)S^1\left(\frac{\Delta t}{2}\right)u\left(t\right)+O\left({\left(\Delta t\right)}^2\right)$ (5)

2. 分数阶Swift-Hohenberg方程的算子分裂法

2.1. Laplace变换

Laplace变换 [11] 是工程数学中常用的一种积分变化，在科学研究领域有广泛的应用。

$F\left(s\right)=L\left[f\left(t\right)\right]={\displaystyle {\int }_0^{\infty }f\left(t\right){\text{e}}^{-st}\text{d}t}$ (6)

运用分部积分法和Laplace变换公式可以得到以下微分性质：

$L\left[f^1\left(t\right)\right]=sL\left[f\left(t\right)\right]-f\left(0\right)=sF\left(s\right)-f\left(0\right)$ (7)

将该微分性质推广到一般情况，可以得到：

$L\left[f^n\left(t\right)\right]=s^{n}F\left(s\right)-s^{n-1}f\left(0\right)-s^{n-2}{f}^1\left(0\right)-\cdots -s{f}^{n-2}\left(0\right)-f^{n-1}\left(0\right)$ (8)

Laplace逆变换就是实现以上过程的逆向变换。

$f\left(t\right)=L^{-1}\left[F\left(s\right)\right]$ (9)

研究发现 [12] [13]，利用Laplace变换可以很好地将分数阶方程转化为整数阶方程，进而降低求解难度。现利用Laplace变换逼近Caputo型分数阶导数：

$L\left[{}_0{}^{c}D{}_t^{\alpha }u\left(x,t\right)\right]=s^{\alpha }L\left[u\left(x,t\right)\right]-s^{\alpha -1}u\left(x,0\right)$ (10)

利用Lagrange插值实现 $s^{\alpha }$ 线性化，得到：

$s^{\alpha }\approx \alpha s^1+\left(1-\alpha \right)s^0=\alpha s+\left(1-\alpha \right)$ (11)

将式(11)代入式(10)，得到：

$L\left[{}_0{}^{c}D{}_t^{\alpha }u\left(x,t\right)\right]\approx \left[\alpha s+\left(1-\alpha \right)\right]\left[L\left[u\left(x,t\right)\right]-s^{-1}u\left(x,0\right)\right]$ (12)

对上式进行Laplace逆变换，整理后得到：

${}_0{}^{c}D{}_t^{\alpha }u\left(x,t\right)\approx \left(1-\alpha \right)\left[u\left(x,t\right)-u\left(x,0\right)\right]+\alpha \frac{\partial u\left(x,t\right)}{\partial t}$ (13)

于是，时间分数阶Swift-Hohenberg方程可以转化为如下的整数阶方程：

$\frac{\partial u\left(x,t\right)}{\partial t}=-\frac{{\left(1+\text{Δ}\right)}^{2}u}{\alpha }+\frac{\alpha -1}{\alpha }\left[u\left(x,t\right)-u\left(x,0\right)\right]-\frac{f\left(u\right)}{\alpha }$ (14)

2.2. 算子分裂法

算子分裂法 [14] 是一种偏微分方程数值解法，可以把复杂的算子分裂成几个可分别求解的简单子算子之积，实现复杂问题的简单化处理。

运用算子分裂法分解成两个子问题时，算子的顺序不影响计算结果，具体求解过程如下：

(1) 计算前半步长的算子 $S^1$，其中 $t\in \left[t_k,t_{k+1/2}\right]$；

(2) 令 $u^{t_k}=u^{t_{k+1/2}}$，计算一个步长的算子 $S^2$，其中 $t\in \left[t_k,t_{k+1}\right]$；

(3) 令 $u^{t_{k+1/2}}=u^{t_{k+1}}$，计算后半步长的算子 $S^1$，其中 $t\in \left[t_{k+1/2},t_{k+1}\right]$。

取空间节点数为M，则空间步长 $h=\left(b-a\right)/M$；取时间节点数为N，则时间步长 $\tau =T/N$。

空间方向的节点： $x_i=a+ih,i=0,1,2,\cdots ,M$

时间方向的节点： $t_k=k\tau ,k=0,1,2,\cdots ,N$

利用算子分裂法将转化后的整数阶方程进行如下分解：

$\{\begin{array}{l}\frac{\partial u}{\partial t}=-\frac{1}{\alpha }f\left(u\right)=-\frac{1}{\alpha }\left(u^3-\eta u^2-\epsilon u\right)(15)\\ \frac{\partial u}{\partial t}=-\frac{{\left(1+\frac{{\partial }^2}{\partial x^2}\right)}^{2}u}{\alpha }+\frac{\alpha -1}{\alpha }\left[u\left(x,t\right)-u\left(x,0\right)\right](\; 16\; )\end{array}$

对于非线性方程式(15)，方程两边同除以u，然后两边同时积分，可得：

$u^{k+1}=u^k\exp \left[-\frac{\tau }{\alpha }\left({\left(u^k\right)}^2-\eta u^k-\epsilon \right)\right]$ (17)

对于线性方程式(16)，利用Crank-Nicolson格式在时间方向的(k+1/2)层处逼近，可得：

$u^{k+\frac{1}{2}}=\frac{u^k+u^{k+1}}{2}$ (18)

$\frac{u^{k+1}-u^k}{\tau }=-\frac{1}{\alpha }\left(\frac{u_{x^4}^k+u_{x^4}^{k+1}}{2}+u_{x^2}^k+u_{x^2}^{k+1}+\frac{u^k+u^{k+1}}{2}\right)+\frac{\alpha -1}{\alpha }\left(\frac{u^k+u^{k+1}}{2}-u^0\right)$ (19)

注意， $u_{-1}$ 和 $u_{N+1}$ 利用线性插值给出，即：

$\{\begin{array}{l}{u}_{-1}=2u_0-u_1\\ u_{N+1}=2u_N-u_{N-1}\end{array}$ (20)

两端边值条件均为零。

将已知的k 层和待求的k + 1层分别整理到等式的两边，得到矩阵：

$\left[\begin{array}{ccccc}c-a& b& a& & \cdots \\ b& c& b& a& \\ & \ddots & \ddots & \ddots & \\ & a& b& c& b\\ \cdots & & a& b& c-a\end{array}\right]\left[\begin{array}{c}{u}_1^{k+1}\\ u_2^{k+1}\\ ⋮\\ u_{M-2}^{k+1}\\ u_{M-1}^{k+1}\end{array}\right]=\left[\begin{array}{ccccc}d+a& -b& -a& & \cdots \\ -b& -d& b& -a& \\ & \ddots & \ddots & \ddots & \\ & -a& -b& d& -b\\ \cdots & & -a& -b& d+a\end{array}\right]\left[\begin{array}{c}{u}_1^k\\ u_2^k\\ ⋮\\ ⋮\\ u_{M-1}^k\end{array}\right]+\frac{\left(1-\alpha \right)\tau }{\alpha }\left[\begin{array}{c}{u}_1^0\\ u_2^0\\ ⋮\\ ⋮\\ u_{M-1}^0\end{array}\right]$ (21)

其中 $a=\frac{\tau }{2\alpha h^4}$， $b=a\left(2h^2-4\right)$， $c=1+\frac{\tau }{2\alpha }+\frac{\left(1-\alpha \right)\tau }{2\alpha }+a\left(6-4h^2\right)$， $d=1-\frac{\tau }{2\alpha }-\frac{\left(1-\alpha \right)\tau }{2\alpha }-a\left(6-4h^2\right)$。

记左右两边矩阵分别为A和B，则待求的 $k+1$ 层表示如下：

$u_i^{k+1}=A^{-1}\left[B{u}_i^k+\frac{\left(1-\alpha \right)\tau }{\alpha }{u}_i^0\right]$ (22)

将式(17)和式(22)代入式(5)，得到原方程的数值解：

$\{\begin{array}{l}\left[t_k,t_{k+\frac{1}{2}}\right]{\left(u_i^k\right)}^{\left(1\right)}=u_i^{k+\frac{1}{2}}=u_i^k\exp \left[-\frac{\tau }{2\alpha }\left({\left(u_i^k\right)}^2-\eta u_i^k-\epsilon \right)\right]\\ \left[t_k,t_{k+1}\right]{\left(u_i^{k+1}\right)}^{\left(1\right)}=A^{-1}\left[B{\left(u_i^k\right)}^{\left(1\right)}+\frac{\left(1-\alpha \right)\tau }{\alpha }{u}_i^0\right]\\ \left[t_{k+\frac{1}{2}},t_{k+1}\right]u_i^{k+1}={\left(u_i^{k+1}\right)}^{\left(1\right)}\exp \left[-\frac{\tau }{2\alpha }\left({\left({\left(u_i^{k+1}\right)}^{\left(1\right)}\right)}^2-\eta {\left(u_i^{k+1}\right)}^{\left(1\right)}-\epsilon \right)\right]\end{array}$ (23)

3. 数值算例

3.1. 精度测试

验证该数值格式的精确性时，在原方程的基础上引入右端项 $g\left(x,t\right)$。

$\{\begin{array}{l}{}_0{}^{c}D{}_t^{\alpha }u\left(x,t\right)=-\left[{\left(1+\frac{{\partial }^2}{\partial x^2}\right)}^{2}u+f\left(u\right)\right]+g\left(x,t\right)(24)\\ u\left(x,t\right)=t^2\sin \left(\pi x\right)(\; 25\; )\end{array}$

特殊的初值函数求分数阶导，有如下求解公式：

${}_0{}^{c}D{}_t^{\alpha }{t}^n=\frac{n!}{\text{Γ}\left(n+1-\alpha \right)}{t}^{n-\alpha }$ (26)

代入公式求解式(26)的分数阶项：

${}_0{}^{c}D{}_t^{\alpha }{t}^2\sin \left(\pi x\right)=\frac{2}{\text{Γ}\left(3-\alpha \right)}{t}^{2-\alpha }\sin \left(\pi x\right)$ (27)

移项整理得到 $g\left(x,t\right)$ 的表达式：

$g\left(x,t\right)=\frac{2}{\text{Γ}\left(3-\alpha \right)}{t}^{2-\alpha }\sin \left(\pi x\right)+{\left(t^2\sin \left(\pi x\right)\right)}^3-\eta {\left(t^2\sin \left(\pi x\right)\right)}^2+\left[{\left({\pi }^2-1\right)}^2-\epsilon \right]t^2\sin \left(\pi x\right)$ (28)

选取区域范围 $\left[0,1\right]\times \left[0,1\right]$，网格剖分 $M=20$， $N=1000$， $\alpha =0.5$， $\eta =0.5$， $\epsilon =0.2$，给出数值解图像和误差图像，分别如图1，图2所示。

Figure 1. Numerical solution image

图1. 数值解图像

Figure 2. Error image

图2. 误差图像

选定时间方向的最后一层 $t=T$，给出最大相对误差 $E\left(M,N\right)$ 的计算公式：

$E\left(M,N\right)=\underset{0\le i\le M-1}{\max }\frac{\left|u\left(x_i,T\right)-u_i^N\right|}{\left|u\left(x_i,T\right)\right|}$ (29)

其中， $u\left(x_i,t_k\right)$ 和 $u_i^k$ 分别对应节点 $\left(x_i,t_k\right)$ 的精确解和数值解。

取定不同的时间和空间网格剖分，并改变参数值，列表观察最大相对误差变化，如表1所示。

Table 1. The maximum relative error of different time-space meshing and different parameter values

表1. 不同时空网格剖分和不同参数值的最大相对误差

观察发现，数值解逼近精确解具有较高的精度，并且空间方向网格剖分越密， $\alpha$ 越小，精度越高。验证了本文所采用格式的有效性和精确性。

3.2. 长时间动力学行为

考虑如下初值问题：

$u\left(x,0\right)=\sin \left(\pi x\right),x\in \left[0,6\right]$ (30)

两端的边值条件均为零，取定参数 $\eta =0.5$， $\epsilon =0.2$， $t\in \left[0,100\right]$，时空剖分 $\left[M,N\right]=\left[40,1000\right]$，给出 $\alpha =0.1$ 的数值解图像，时间方向最后一层 $t=T$ 的数值解曲线，分别如图3，图4所示。

Figure 3. Numerical solution image ( $\alpha =0.1$ )

图3. 数值解图像( $\alpha =0.1$ )

Figure 4. Numerical solution curve of $t=T$

图4. $t=T$ 的数值解曲线

4. 结论

基于Laplace变换，将时间分数阶Swift-Hohenberg方程转化为整数阶Swift-Hohenberg方程，然后，利用算子分裂法将其进一步分解成两个子问题分别求解，最后通过数值算例验证格式的有效性。本文对Swift-Hohenberg方程的求解方法，为今后该方程的深入研究提供了有力的支持。

基金项目

2019年安徽省高校学科(专业)拔尖人才学术资助项目(gxbjZD38)。