改进Douglas分裂方法求解反应扩散方程的全离散误差分析①

李平，陈浩

重庆师范大学 数学科学学院， 重庆 401331

反应扩散方程的应用非常广泛， 比如斑图动力学[1]、图灵结构[2]、非线性波[3]、孤子或螺旋波[4]、时空混沌[5]. 然而， 由于这类系统存在刚性项和非线性项的耦合， 所以想要有效和精确地模拟这些系统是十分困难的. 到目前为止， 关于反应扩散方程数值解法的研究已经有了很多[6-18]， 例如： 隐显式方法[6]、隐显式Runge-Kutta方法[7]、线性隐式Runge-Kutta方法[8]、指数时间差分法[9]、混合Runge-Kutta方法[10]、指数积分法[11]等.

本文主要利用改进Douglas分裂方法[12]求解反应扩散方程， 该方法具有良好的稳定性且计算速度很快. 但在该文献中， 作者只对线性反应扩散方程进行了时间半离散误差分析得到了二阶收敛的结论. 对于非线性反应扩散方程的全离散误差分析在文献中尚未提到. 本文的主要目的是对非线性反应扩散方程进行全离散误差分析.

本文的结构安排如下： 第一节主要介绍二维非线性反应扩散方程空间、时间离散方法； 第二节给出了所构造的数值格式的全离散误差分析； 第三节给出了几个空间二维及三维的数值算例； 最后在第四节给出简短的总结.

1 离散方法

1.1 空间离散

我们考虑配备齐次Dirichlet边界条件或齐次Neumann边界条件下的二维非线性反应扩散方程

(1)

从而得到半离散系统

(2)

其中：uj，k(t)为u(xj，yk，t)的近似解，fj，k(t，uj，k(t))： =f(xj，yk，t，uj，k(t)). 结合齐次Dirichlet边界条件， 半离散化系统(2)可改写成矩阵形式

U′=AU+F(t，U)

(3)

U=[u1， 1(t)， …，uN， 1(t)，u1， 2(t)， …，uN， 2(t)， …，u1， N(t)， …，uN， N(t)]T

F(t，U)=[f1， 1(t，u1， 1(t))， …，fN， 1(t，uN， 1(t))， …，f1， N(t，u1， N(t))， …，fN， N(t，uN， N(t))]T

若为齐次Neumann边界条件时，

1.2 时间离散

在进行时间离散时， 利用改进Douglas分裂方法[12]求解半离散系统(3)， 将其改写为

U′=A1U+A2U+F(t，U)

(4)

2 收敛性分析

由半离散系统(3)进行分析得方程(1)的精确解格式

(5)

(6)

其中， eτA是矩阵指数. 令s=tn+v， 利用中点公式， 则等式(6)可化为

(7)

由指数函数有理逼近得

所以可得到精确解表达为

(8)

消去式(4)中的中间变量得到数值解表达式为

(9)

假设1假设F(t，U)满足Lipschitz条件， 即存在L>0， 对∀U1，U2∈RN2， ∀t∈[t0，T]， 有

‖F(t，U1)-F(t，U2)‖≤L‖U1-U2‖

引理1对∀τ>0有

证因为A1对称负定， 则A1的特征值μk<0，k=1，…，N2. 于是，

同理可证A2的情况.

定理1若假设1成立， 则Modified Douglas Splitting方法是二阶收敛的， 即

其中c为正常数.

证将精确解表达式(8)与数值解表达式(9)作差后取范数， 由引理1和Lipschitz条件得

(10)

其中c2为正常数.由

得到

因为

所以

则有

(11)

其中：ω=‖A1+A2‖，c3为正常数.将式(11)代入式(10)可得

en+1≤(1+βτ)en+c1τh2+c3τ3+c1τ2h2+c2τ3

依此类推得

3 数值实验

下面给出反应扩散方程的数值算例， 主要对收敛性分析进行数值验证， 实验是通过Matlab实现的.

算例1考虑二维反应扩散方程 (1)， 其中，Ω=[0， 1]2，t∈[0， 1]. 边界条件为齐次Dirichlet边界条件， 以及初始条件

u(x，y， 0)=sin(2πx)sin(2πy)

其中， 扩散项为

f(x，y，t，u)=u-u3+e-t(8π2-2)sin(2πx)sin(2πy)+e-3tsin3(2πx)sin3(2πy)

表1 二维反应扩散方程计算结果

通过表1的结果， 我们可以得出误差关于h，τ是二阶精度的， 与收敛性分析是一致的.

算例2考虑三维反应扩散方程

其中：Ω=[0， 1]3，t∈[0， 1]. 边界条件为齐次Dirichlet边界条件， 以及初始条件

u(x，y，z， 0)=sin(2πx)sin(2πy)sin(2πz)

其中， 扩散项为

f(x，y，z，t，u)=u-u3+e-t(12π2-2)sin(2πx)sin(2πy)sin(2πz)+e-3tsin3(2πx)sin3(2πy)sin3(2πz)

表2 三维反应扩散方程计算结果

通过上表的结果， 我们可以得到三维反应扩散方程可以利用同样的方法进行求解， 得到的误差关于h，τ也是二阶精度的.

算例3考虑二维Schnackenberg方程组

其中：Ω=[0，L]2，t∈[0，T]. 边界条件为齐次Neumann边界条件， 以及初始条件

该方程组没有精确解， 利用该方法进行求解.取N=101，M=8 000时， 其中，a=0.130 5，b=0.769 5，γ=100，Ku=0.05，Kv=1， 当L=1，T=0.5， 1， 2， 3时， 得到数值解u(x，y，t)的图像如图1所示.

图1 数值解

4 结论

参考文献[12]中针对反应扩散方程提出了一类二阶改进Douglas分裂方法， 该方法具备良好的稳定性且计算速度快的特点， 但作者仅仅给出了该方法针对线性问题的半离散误差分析. 而本文主要采用空间二阶中心差分方法对空间方向进行离散， 利用改进Douglas方法对时间方向进行离散得到相应的Modified Douglas Splitting全离散格式. 对该格式进行收敛性分析， 证明该全离散格式关于空间和时间步长是二阶收敛的结论. 最后借助相关数值实验算例进行收敛性验证.