基于快速傅里叶变换的有限差分方法求解三维反应扩散方程

摘要：基于快速傅里叶变换求解齐次Neumann边界条件下的三维非线性反应扩散方程，应用有限差分方法给出二阶中心差分格式，利用Kronecker积的性质将三维拉普拉斯算子的微分矩阵进行对角化处理，得到相对应的特征值与特征向量;在时间离散上采用Crank-Nicolson方法，并采用Picard迭代求解离散得到的非线性代数方程组。结果发现，利用快速傅里叶变换求解Allen-Cahn方程，随着时间推移，显示出解从初始状态、过渡层、亚稳态进而到达到稳态的演化过程。最后，给出数值算例，验证了所用方法求解三维反应扩散方程可在保持精度的同时，减少存储量，并可大幅度降低计算时间。

关 键 词：反应扩散方程; 有限差分; Crank-Nicolson方法; 快速傅里叶变换

中图分类号：O241.82 文献标志码：A

doi：10.3969/j.issn.1673－5862.2023.01.012

Finite difference method based on fast Fourier transform for solving three dimensional reaction diffusion equation

ZHANG Rongpei1， LIU Hao 2， ZUO Hanxing3

（1. School of Mathematics and Systematic Sciences， Shenyang Normal University， Shenyang 110034， China;

2. Advanced Manufacturing College， Guangdong University of Technology， Guangzhou 510006， China ;

3. School of Preschool and Primary Education， Shenyang Normal University， Shenyang 110034， China）

Abstract：To solve the three-dimensional nonlinear reaction-diffusion equation under homogeneous Neumann boundary conditions based on Fast Fourier transform（taking Allen Cahn equation as an example）， the second order central difference scheme is given by using the finite difference method. The differential matrix of the three-dimensional Laplace operator is diagonalized by using the properties of Kronecker product， and the corresponding eigenvalues and eigenvectors are obtained; Crank Nicolson method is used for time discretization， and Picard iteration is used to solve the discretized nonlinear algebraic equations. It is found that the three-dimensional reaction-diffusion equation shows the evolution process from the initial state， the transition layer， the metastable state to the steady-state solution with the passage of time. Finally， an example is given through numerical experiments. Using the finite difference method of fast Fourier transform to solve the three-dimensional reaction-diffusion equation can not only maintain the accuracy and simplicity， but also reduce the storage and calculation time.

Key words：reaction diffusion equation; finite difference; Crank-Nicolson method; fast Fourier transform

3 数值实验

本节分为2个部分，主要进行齐次Neumann边界条件下的三维Allen-Cahn方程的稳定性检验实验和精确解的精度测试实验。为了便于后文叙述，在2个数值实验中，令hx=hy=hz=h，应用本文提出的有限差分法求解下面的数值算例。

3.1 精确解的精度测试

为了对收敛速度进行精确的定量估计，本文在齐次Neumann边界条件下分别对4组较细的网格、5组不同时间的初始问题进行多次模拟。利用本文提出的FFT方法的有限差分方法，在均匀的网格和时间步长上进行数值计算。选取初值为u0（x，y，z）=cos（πx）+cos（πy）+cos（πz），时间步长为Δt=10-3s，网格步长N分别取8，16，32和64，时间t分别取0， 0.1， 0.2和0.6 s。表1给出了齐次Neumann边界条件下精确解的误差精度。从表1可以发现，FFT方法确实在时间上具有二阶精度，同时准确地反映了本文方法能得到收敛性较低的精确解的优点。

3.2 精确解的精度测试

初值条件选取为u0（x，y，z）=exp（cos（1.5πx）cos（1.5πy）sin（πz）sin（2πz）），网格步长为64，计算区域Ω=［-1，1］3，针对具有齐次Neumann边界条件的方程（1）进行求解。时间离散步长为Δt=10-3，对不同的时间点t进行数值模拟。

在齐次Neumann边界条件下，随着时间t的推移，[JY]图1分别显示了三维反应扩散方程的解从初始状态到过渡层，再到亚稳态，最后到达稳态的演化过程。从图1中可以看到，解初始状态到亚稳态是个快速动态过程，并在亚稳态中形成了2个过渡层。图1（b）和图1（c）表明初始状态下的分层较多，过渡层中的2个界面在相分离的情况下，自由边界条件不太严格，进而实现速度更快，为后面达到稳定状态奠定了良好的基础;图1（b）表明在Neumann边界条件下，在t=0.1 s时，解处于过渡层;图1（c）表明在t=0.2 s时，解处于亚稳态;图1（d）展示了在求解Allen-Cahn方程的过程中，在t=0.8 s时，解就迅速形成了稳态。

4 结 语

本文给出了齐次Neumann边界条件下三维反应扩散方程的一种新的线性化有限差分格式。由于FFT方法能够快速地求解经过离散得到的非线性代数方程组，故在提高计算效率的同时也大大节省了计算时间。通过求解方程，可以发现随时间推移数值结果形状的变化规律，即随时间推移，初始状态、过渡层、亚稳态进而达到稳态的解的演化过程。本文的处理方法在保持精度和简洁性的同时，大大降低了算法的复杂度，使得程序运算时间明显缩短，对于进行长时间尺度动力学模拟具有重要意义。

参考文献：

［ 1 ］田娅，秦瑶，向晶. 一类带有变指数非局部项的反应扩散方程解的爆破行为［J］. 应用数学和力学， 2022，43（10）：1177－1184.

［ 2 ］彭红玲，樊明书. 一类半线性分数阶反应扩散方程解的性质［J］. 四川师范大学学报：自然科学版， 2022（6）：2－9.

［ 3 ］李建军，王看看，徒君. 具有奇异势的拟线性分数阶扩散方程解的爆破性［J］. 应用数学， 2022，35（4）：819－826.

［ 4 ］邹享明. 基于指数变换的对流占优反应扩散问题的有限体积法［J］. 科学技术创新， 2022（22）：53－58.

［ 5 ］蒋沁纱，陈浩. 空间变系数反应扩散方程的一类交替分裂预处理迭代方法［J］. 重庆师范大学学报：自然科学版， 2022，39（5）：83－90.

［ 6 ］王妍. 非线性对流反应扩散方程爆破解的B方法［J］. 纺织高校基础科学学报， 2022，35（2）：83－91.

［ 7 ］冯周平. 变系数分数阶Non-Fickian反应扩散方程的有限差分方法研究［D］. 成都： 四川师范大学， 2022.

［ 8 ］王云霞. Allen-Cahn方程和Cahn-Hilliard-Hele-Shaw方程的有限元数值算法研究［D］. 太原： 太原理工大学， 2022.

［ 9 ］GILBERT S. Computational Science and Engineering［M］. Massachusetts： Wellesly-Cambridge Press， 2007：78－82.

［10］张荣培，霍俊蓉，杨程程. 基于离散余弦变换的积分因子方法求解非线性Allen-Cahn方程［J］. 沈阳师范大学学报：自然科学版， 2021，39（2）：159－163.

［11］ZHAI S，FNG X，He E Y. Numerical simulation of the three dimensional Allen-Cahn equation by the high-order compact ADI method［J］. CPC， 2014，185（10）：2449－2455.

收稿日期：2022－11－18

基金项目：辽宁省自然科学基金资助项目（20180550996）。

作者简介：张荣培（1978—），男，山东泰安人，广东工业大学副教授，博士;

通信作者：刘 昊（1995—），男，宁夏吴忠人，沈阳师范大学在读硕士研究生。