收稿日期:20221118
基金项目：辽宁省自然科学基金资助项目(20180550996
)㊂作者简介：张荣培(1978 )，男，山东泰安人，广东工业大学副教授，博士；通信作者：刘 昊(1995 )
reaction diffusion，男，宁夏吴忠人，沈阳师范大学在读硕士研究生㊂
第41卷 第1期2023年 2月
沈阳师范大学学报(自然科学版)
J o u r n a l o f S h e n y a n g N o r m a lU n i v e r s i t y (
N a t u r a l S c i e n c eE d i t i o n )V o l .41N o .1F e b .202
3
文章编号:16735862(2023)01007305
基于快速傅里叶变换的有限差分方法
求解三维反应扩散方程
张荣培1,
2
，刘 昊2，左韩星3，温学兵1(1.沈阳师范大学数学与系统科学学院，沈阳 110034;2.广东工业大学先进制造学院，广州 510006;
3.沈阳师范大学学前与初等教育学院，沈阳 110034
)摘  要：基于快速傅里叶变换求解齐次N e u m a n n 边界条件下的三维非线性反应扩散方程，应用有限差分方法给出二阶中心差分格式，利用K r o n e c k e r 积的性质将三维拉普拉斯算子的微分矩阵进行对角化处理，得到相对应的特征值与特征向量；在时间离散上采用C r a n k -N i c o l s o n 方法，并采用P
i c a r d 迭代求解离散得到的非线性代数方程组㊂结果发现，利用快速傅里叶变换求解A l l e n -C a h n 方程，随着时间推移，显示出解从初始状态㊁过渡层㊁亚稳态进而到达到稳态的演化过程㊂最后，给出数值算例，验证了所用方法求解三维反应扩散方程可在保持精度的同时，减少存储量，并可大幅度降低计算时间㊂
关 键 词：反应扩散方程；有限差分;C r a n k -N i c o l s o n 方法；快速傅里叶变换中图分类号:O 241.82    文献标志码:A
d o i :10.3969/j .i s s n .16735862.2023.01.012F i n i t ed i f f
e r e n c e m e t h o db a s e do n
f a s tF o u r i e rt r a n s f o r m f o r
s o l v i n g t h r e e d i m e n s i o n a l r e a c t i o nd i f f u s i o n e q
u a t i o n Z HA N GR o n g p e i 1,L I U H a o 2,Z U O H a n x i n g
3
(1.S c h o o l o fM a t h e m a t i c sa n dS y s t e m a t i cS c i e n c e s ,S h e n y a n g N o r m a lU n i v e r s i t y ,S h e n y a n g 1
10034,C h i n a ;2.A d v a n c e d M a n u f a c t u r i n g C o l l e g e ,G u a n g d o n g U n i v e r s i t y o f T e c h n o l o g y ,G u a n g z h o u 510006,C h i n a ;3.S c h o o l o f P r e s c h o o l a n dP r i m a r y E d u c a t i o n ,S h e n y a n g N o r m a lU n i v e r s i t y ,S h e n y a n g 1
10034,C h i n a )A b s t r a c t :T o s o l v e t h e t h r e e -d i m e n s i o n a l n o n l i n e a r r e a c t i o n -
d i f f u s i o n
e q u a t i o n u n d e r h o m o g e n e o u sN e u m a n nb o u n d a r y c o n d i t i o n sb a s e do nF a s tF o u r i e r t r a n s
f o r m (t a k i n
g A
l l e nC a h n e q u a t i o na sa ne x a m p l e ),t h es e c o n do r d e rc e n t r a ld i f f e r e n c es c h e
m e i s g i v e nb y u s i n g t h ef i n i t e d i f f e r e n c em e t h o d .T h e d i f f e r e n t i a lm a t r i x o f t h e t h r e e -d i m e n s i o n a l L a p l a c e o p e r a t o r i s d i a g o n a l i z e d b y u s i n g t h e p r o p e r t i e so fK r o n e c k e r p r o d u c t ,a n dt h ec o r r e s p o n d i n g e i g e n v a l u e sa n de i g
e n v e c t o r s a r e o b t a i n e d ;C r a n kN i c o l s o nm e t h o d i s u s e d
f o r t i m e d i s c r e t i z a t i o n ,a n dP i c a r d i t e r a t i o n i s u s e d t o s o l v e t h e d i s c r e t i z e dn o n l i n e a r a l
g e b r a i ce q u a t i o n s .I t i s f o u n dt
h a t t h e t h r e e -d
i m e n s i o n a l r e a c t i o n -d i f f u s i o ne q u a t i o ns h o w st h ee v o l u t i o n p r o c e s sf r o m t h ei n i t i a ls t a t e ,t h et r a n s i t i o nl a y e r ,t h e m e t a s t a b l e s t a t e t o t h e s t e a d y -s t a t e s o l u t i o nw i t h t h e p a s s a g e o f t i m e .F i n a l l y ,a n e x a m p l e i s g i v e n t h r o u g hn u m e r i c a le x p e r i m e n t s .U s i n g t h ef i n i t ed i f f e r e n c e m e t h o do f f a s tF o u r i e rt r a n s f o r m t o s o l v et h et h r e e -d i m e n s i o n a lr e a c t i o n -d i f f u s i o ne q u a t i o nc
a n n o to n l y m a i n t a i nt h ea c c u r a c y a n d s i m p l i c i t y ,b u t a l s o r e d u c e t h e s t o r a g
e a n d c a l c u l a t i o n t i m e .K e y w o r d s :r e a c t i o nd i
f f u s i o ne q u a t i o n ;f i n i t ed i f f e r e n c e ;C r a n k -
N i c o l s o n m e t h o d ;f a s tF o u r i e r t r a n s f o r m
在本文中，考虑求解如下齐次N e u m a n n 边界的反应扩散方程:
u t =d (u x x +u y y +u z z )+f (u ),(x ,y ,z )ɪΩ,t >0(1)边界和初始条件如下:
u (x ,y ,
z )췍Ω=0u (x ,y ,z )=u 0(x ,y ,z ),(x ,y ,z )ɪR {
3
(2
)其中:f (u )为一般形式;d 为扩散系数;Ω=[a ,b ]ˑ[c ,d ]ˑ[e ,f ]为一个区域;췍Ω代表边界;u 0(x ,y ,
z )为初始值㊂
方程(1
)在物理学㊁化学㊁生物学㊁奇异摄动问题等方面有广泛的应用，尤其是在非线性热传导㊁半导体中的电子和空穴流㊁燃烧理论㊁B e l o u s o v -
Z h a b o t i n s k i i 反应以及神经轴突中电脉冲的传导H o d g k i n -H u x l e y 方程中更为重要㊂由于式(1)右端非线性项的存在，要求出精确解具有很大的难度，所以关于该问题数值解的研究就具有非常重要的意义㊂近年来，众多学者发展了丰富的数值方法求解
该方程㊂田娅等[1]利用上下解方法求解了一类带有变指数非局部项的反应扩散方程；彭红玲和樊明书[2]提出了C a f f a r e l l i -S i l v e s t r e 的延拓方法㊁凸方法；李建军等[3]利用凹函数法求解拟线性分数阶反应扩散方程；邹享明[4]提出了高精度的有限体积元格式；蒋沁纱和陈浩[5]提出了交替分裂预处理迭代方法；王妍[6]发展了B 方法；冯周平[7]提出了变系数分数阶N o n -
F i c k i a n 反应扩散方程的有限差分方法㊂但上述方法仅针对二维情形，同时对三维非线性扩散方程求解比较少，而且由于三维计算复杂，避免不了求解较为庞大的非线性代数方程组，算法冗长且耗费时间较多㊂本文应用K r o n e c k e r 积的性质将三维拉普拉斯算子的微分矩阵进行特征分解，得到其特征解与特征方程，再结合快速傅里叶变换(f a s t
F o u r i e r t r a n s f o r m ,F F T )
方法进行求解，在大大提高精确度的同时，节省了计算时间㊂1 有限差分格式
首先用直线x =x l ,y =y j ,
z =z k ，其中l =0,1, ,N x ,j =0,1, ,N y ,k =0,1, ,N z ，在Ω上打网格，将计算区域Ω=[a ,b ]ˑ[c ,d ]ˑ[e ,f ]
划分为(N x +1)ˑ(N y +1)ˑ(N z +1)个网格点，其中网格点坐标为T h ={(x l ,y j ,z k )=a +l h x ,c +j h y ,
e +k h z },h 2=
f -e N 2
㊂设u (x ,y ,z )在网格点(x l ,y j ,z k )的数值解为u l ,j ,
k ，定义离散解的矩阵U =(u l ,j ,k )1ɤl ɤN x ,1ɤj ɤN y ,1ɤj ɤN z
，由差分格式和边界条件，可得到x ,y ,z 方向的微分矩阵如下:G m =1h 2m 1-100 0-12-10 00-12-1 0︙︙⋱⋱⋱
︙0
00-12-10
000-éëêêêêê
êêêùûú
ú
ú
ú
úú
úú
11N m ˑN m (3
)其中m =x ,y ,
z ㊂将这3个矩阵G m 对角化后分别表示为G x =S x Λx S -1x ,G y =S y
Λy S -1y ,G z =S z Λz S -1
z ，这里的特征向量矩阵S m 为离散傅里叶变换矩阵,
是正交矩阵[8
]，即S -1m =S T m ,Λm 为特征值组成的对角阵，可表示成如下的形式:
Λx =d i a g {λ1,λ2, ,λN x },λl =-2+2c o s l 2πh x Λy =d i a g {λ1,λ2, ,λN y },λy =-2+2c o s j 2πh y
Λz =d i a g {λ1,λ2, ,λN z
},λz =-2+2c o s k 2πh ìîíï
ï
ïïz (4
)(S x )l j k =1N x e -i 2πN x l j k ,(S y )l j k =1N y e -i 2πN y l j k ,(S z )l j
k =1N z
e -i 2π
N z l j k
在M A T L A B 中,
离散傅里叶变换矩阵对向量的乘积相当于对向量进行了快速傅里叶变换实现[10
]㊂设I m 分别为N m 阶单位矩阵,
将解矩阵向量化后得u =V e c (U )=(u 1,1,1, ,u 1,1,N z ,u 1,2,1, ,u 1,2,N z , ,u 1,N y ,N z ,u 2,1,1, ,u N x ,N y ,N z
)T (5
)  利用式(3)以及K r o n e c k e r 积的定义，式(1
)的二阶中心差分格式可以写成矩阵形式47沈阳师范大学学报(自然科学版)              第41卷
d u d t
=(
I z 췍I y 췍G x +I z 췍G y 췍I x +G z 췍I y 췍I
x )u +F (U )(6
)  设计算终止时间为T ，时间步长τ=Δt N
，则第n 层时间步长t n =
n τ,n =0,1,2, ,N ㊂接下来应用C r a n k -N i c o l s o n 格式对式(6)进行求解:U n +1-U n Δt -G U n +1+U n 2=
F n +F n +
12
(7)其中G =I z 췍I y 췍G x +I z 췍G y 췍I x +G z 췍I y 췍I x 是拉普拉斯算子的微分矩阵㊂上式等价于U n +1=G Δt -I æèçöø÷2-1G Δt +I æèçöø÷2U n +G Δt -I æèçöø÷
2-1F n +F n +12
(8)  采用P i c a r d 迭代求解非线性代数方程组得
U k +1n +
1=
I Δt -G æèçöø÷2-1I Δt +G æèçöø÷2U n +I Δt -G æèçöø÷
2-1F (U n )+F (U n +1)k 2
(9
)2 快速求解
在式(8)中需要求矩阵I Δt -G æèçöø÷2T I Δt +G æèçöø÷2以及I Δ
t -G æèçöø÷2T
与向量的乘积,2个矩阵乘积的运算相似，下面以I Δt -췍g æè
çöø÷2-1
I Δt +췍g æèçö
ø÷2v e c (U )
为例进行说明㊂为了快速计算，首先将矩阵G 进行分解㊂下面给出2个定理㊂
定理1 对矩阵G x 进行对角化，定义对角矩阵Λ=d i a g (λ1,λ2, ,λN x ),λl =1h 2x
(-2+2c o s l 2πh x )
,l =1,2, ,N x ，则G x S x =S x Λ，即G x =S x ΛS -1x ，其中(S x )l j
k =1N x
e -i 2π
N x l j k
㊂定理2 K r o n e c k e r 积的性质[9]
:
1)(A 췍B )(C 췍D )=A C 췍B D
2)(A 췍B )-1
=A -1췍B -1
3)(B T 췍A )v e c (X )=v e c (A ˑB )
4)(A 췍B )T
=A T 췍B T
5)d i a g
(v e c (A ))ˑ(v e c (B ))=v e c (A ˑB )接下来利用定理1和2对矩阵G 进行特征分解㊂由I x =S T
x
I x S x ，利用定理2㊁性质1和性质2可得I z 췍I y 췍G x =(S z I z S -1z )췍(S y I y
S -1y )췍(S x Λx S -1x )=(S z 췍S y 췍S x )(I z 췍I y 췍Λx )(S -1z 췍S -1y 췍S -1x )(10
)对I z 췍G y 췍I x ,G z 췍I y 췍
I x 可做类似运算㊂故可将式(6)改写为G =(S z 췍S y 췍S x )췍Λ(S -1z 췍S -1y 췍S -1x )，这里췍Λ=I z 췍I y 췍1Δx 2Λx +I z 췍1Δy
2Λy 췍I x +1Δz
2
Λz 췍I y 췍I x 是一个对角矩阵，则易证明췍Λ=d i a g
(v e c (Λ))(11
)于是三维数组U 可以被视为二维数组上所有这样的一维向量的集合，即U =췍j =1
N y
k =1 N z
U (:,j ,k )，其中U (:,j ,
k )是通过固定最后2个指标j ,k 得到的N x 维向量[11]
㊂由式(11)的分解可得U =G -1+F (U )=(S z 췍S y 췍S x )췍Λ(S -1z 췍S -1y 췍
S -1x )F +F (U )，根据矩阵G 的分解，可得下面的分解
I Δt -g æèçöø÷2-1I Δt +g æèçöø÷2=(S z 췍S y 췍S x )I Δ
t -췍Λæèçöø÷2-1I Δt +췍Λæèçöø÷2(S z 췍S y 췍S x )-1(12)  第1步应用逆离散正弦变换得到向量V 1，即V 1=(S -1z 췍S -1y 췍
S -1
x )F ，若将U 视为二维数组上的一维向量的集合，则有如下形式:
5
7
第1期      张荣培，等：基于快速傅里叶变换的有限差分方法求解三维反应扩散方程
V 1=췍l =1 N x
j =1 N y
i d s t (췍l =1 N x k =1 N z
i d s t (췍j =1
N y
k =1 N z i d s t (F (:,j ,k )))(l ,:,k ))(l ,j ,:)(13
)即V 1=i d s t (i d s t (i d s t (F )))；第2步先定义L ʈ
(l i ,j ,k ),1ɤl ɤN x ,1ɤj ɤN y ,1ɤk ɤN z ，其中l l ,j ,
k =λx l Δx 2+λy j
Δy 2+λz
k Δ
z æèçöø÷2N x ˑN y ˑN z ，且有(Λʈz +Λʈy +Λʈx )(S -1z 췍S -1y 췍S -1x )=V 1./L ʈ，即V 2=V 1./L ʈ㊂其中./为矩阵中每个元素相除，这一步的计算复杂度为N x ˑN y ˑN z 次乘法；第3步可直接应用离散正弦变换
实现，若将U 视为二维数组上的一维向量的集合，则有如下形式
U =췍l =1 N x
j =1 N y
d s t (췍l =1 N x k =1 N z
d s t 췍l =1
N y
k =1 N z
d s t (V 2(:,j ,k ()))(l ,:,k ))(l ,j ,:)即U =d s t (d s t (d s t (V 2)))㊂  在本文提出核心的快速离散正弦变换方法中，第1步和第3步的计算复杂度仅为O (N 2l o g 2N )
，由此可见离散正弦变换方法的计算速度很有优势㊂
3 数值实验
本节分为2个部分，主要进行齐次N e u m a n n 边界条件下的三维A l l e n -
C a h n 方程的稳定性检验实验和精确解的精度测试实验㊂为了便于后文叙述，在2个数值实验中，令h x =h y =h z =h ，应用本文提出的有限差分法求解下面的数值算例㊂3.1 精确解的精度测试
为了对收敛速度进行精确的定量估计，本文在齐次N e u m a n n 边界条件下分别对4组较细的网格㊁5
组不同时间的初始问题进行多次模拟㊂利用本文提出的F F T 方法的有限差分方法,
在均匀的网格和时间步长上进行数值计算㊂选取初值为u 0(x ,y ,z )=c o s (πx )+c o s (πy )+c o s (πz )，时间步长为Δt =10-3s ，网格步长N 分别取8,16,32和64，时间t 分别取0,0.1,0.2和0.6s ㊂表1给出了齐次N e u m a n n 边界条件下精确解的误差精度㊂从表1可以发现,F F T 方法确实在时间上具有二阶精度，同时准确地反映了本文方法能得到收敛性较低的精确解的优点㊂
表1 在Δt =10-3时精确解收敛性测试
T a b l e1 T h ec o n v e r g
e n c e t e s t o
f e x a c t s o l u t i o nw h e nΔt =10-3网格步长t =0.1s C p u
/s t =0.2s
C p u
/s t =0.4s
C p u
/s t =0.8s
C p u
/s 87.44ˑ10-4
5.19
1.51ˑ10-3
5.11
2.91ˑ10-3
10.11
5.60ˑ10-3
20.09
161.89ˑ10-410.233.75ˑ10-436.5397.36ˑ10-439.451.40ˑ10-379.63324.76ˑ10-540.999.43ˑ10-579.981.84ˑ10-4159.823.55ˑ10-4318.7264
1.19ˑ10-5
169.33
2.36ˑ10-5
319.47
4.62ˑ10-5
675.02
8.89ˑ10-5
2854.75
3.2 精确解的精度测试
初值条件选取为u 0(x ,y ,z )=e x p (c o s (1.5πx )c o s (1.5πy )
s i n (πz )s i n (2πz ))，网格步长为64，计算区域Ω=[-1,1]3
，针对具有齐次N e u m a n n 边界条件的方程(1)进行求解㊂时间离散步长为Δt =
10-3,
对不同的时间点t 进行数值模拟㊂在齐次N e u m a n n 边界条件下,
随着时间t 的推移，图1分别显示了三维反应扩散方程的解从初始(a )t =0(b )t =0.1s
67沈阳师范大学学报(自然科学版)              第41卷
(c )t =0.2s (d )t =0.8s
图1 稳定性试验中N e u m a n n 界条件的形态演化与稳态解
F i g .1 M o r p h o l o g i c a l e v o l u t i o na n d s t e a d y -s t a t es o l u t i o no f N e u m a n nb o u n d a r y c o n d i t i o n i n s t a b i l i t y
t e s t 状态到过渡层，再到亚稳态，最后到达稳态的演化过程㊂从图1中可以看到，解初始状态到亚稳态是个快速动态过程，并在亚稳态中形成了2个过渡层㊂图1(b )和图1(c )表明初始状态下的分层较多，过渡层中的2个界面在相分离的情况下，自由边界条件不太严格，进而实现速度更快，为后面达到稳定状态奠定了良好的基础；图1(b )表明在N e u m a n n 边界条件下，在t =0.1s 时，解处于过渡层；图1(c
)表明在t =0.2s 时，解处于亚稳态；图1(d )展示了在求解A l l e n -C a h n 方程的过程中，在t =0.8s 时，解就迅速形成了稳态㊂
4 结  语
本文给出了齐次N e u m a n n 边界条件下三维反应扩散方程的一种新的线性化有限差分格式㊂由于
F F T 方法能够快速地求解经过离散得到的非线性代数方程组,
故在提高计算效率的同时也大大节省了计算时间㊂通过求解方程，可以发现随时间推移数值结果形状的变化规律，即随时间推移，初始状态㊁过渡层㊁亚稳态进而达到稳态的解的演化过程㊂本文的处理方法在保持精度和简洁性的同时，大大降低了算法的复杂度，使得程序运算时间明显缩短，对于进行长时间
尺度动力学模拟具有重要意义㊂参考文献:
[1]田娅，秦瑶，向晶.一类带有变指数非局部项的反应扩散方程解的爆破行为[J ].应用数学和力学,2022,43(10):11771184.
[2]彭红玲，樊明书.一类半线性分数阶反应扩散方程解的性质[J ].四川师范大学学报：自然科学版,2022(6):29. [3]李建军，王看看，徒君.具有奇异势的拟线性分数阶扩散方程解的爆破性[J ].应用数学,2022,35(4):819826.[4]邹享明.基于指数变换的对流占优反应扩散问题的有限体积法[J ].科学技术创新,2022(22):5358.
[5]蒋沁纱，陈浩.空间变系数反应扩散方程的一类交替分裂预处理迭代方法[J ].重庆师范大学学报：自然科学版,2022,39(5):8390.[6]王妍.非线性对流反应扩散方程爆破解的B 方法[J ].纺织高校基础科学学报,2022,35(2):8391.
[7]冯周平.变系数分数阶N o n -F i c k i a n 反应扩散方程的有限差分方法研究[D ].成都：四川师范大学,2022.[8]王云霞.A l l e n -C a h n 方程和C a h n -H i l l i a r d -H e l e -S h a w 方程的有限元数值算法研究[D ].太原：太原理工大学,2022.
[9]G I L B E R T S .C o m p u t a t i o n a lS c i e n c ea n d E n g i n e e r i n g [M ].M a s s a c h u s e t t s :W e l l e s l y -C a m b r i d g e P r e s s ,2007:7882.
[10]张荣培，霍俊蓉，杨程程.基于离散余弦变换的积分因子方法求解非线性A l l e n -C a h n 方程[J ].沈阳师范大学学报：自然科学版,2021,39(2):159163.[11]Z HA I S ,F N G X ,H eE Y.N u m e r i c a l s i m u l a t i o no f t h et h r e ed i m e n s i o n a lA l l e n -C a h ne q u a t i o nb y t h eh i g h -o r d e r c o m p a c tA D Im e t h o d [J ].C P C ,2014,185(10):24492455.7
7
第1期      张荣培，等：基于快速傅里叶变换的有限差分方法求解三维反应扩散方程

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。