A hybrid method for stiff reaction–diffusion equations

Yuchi Qiu; Weitao Chen; Qing Nie

doi:10.3934/dcdsb.2019144

Article Contents

2019, Volume 24, Issue 12: 6387-6417. Doi: 10.3934/dcdsb.2019144

This issue Previous Article Dynamics of a stochastic hepatitis C virus system with host immunity Next Article Global bounded and unbounded solutions to a chemotaxis system with indirect signal production

A hybrid method for stiff reaction–diffusion equations

1.
Department of Mathematics, University of California, Irvine, Irvine, CA 92697, USA
2.
Department of Mathematics, University of California, Riverside, Riverside, CA 92507, USA
3.
Department of Mathematics, Department of Developmental and Cell Biology, University of California, Irvine, Irvine, CA 92697, USA

^* Corresponding author: Qing Nie
^* Corresponding author: Qing Nie

Received: October 2018

Revised: January 2019

Early access: July 2019

Published: September 2019

Abstract / Introduction Full Text(HTML) Figure(3) / Table(9) Related Papers Cited by

Abstract

The second-order implicit integration factor method (IIF2) is effective at solving stiff reaction–diffusion equations owing to its nice stability condition. IIF has previously been applied primarily to systems in which the reaction contained no explicitly time-dependent terms and the boundary conditions were homogeneous. If applied to a system with explicitly time-dependent reaction terms, we find that IIF2 requires prohibitively small time-steps, that are relative to the square of spatial grid sizes, to attain its theoretical second-order temporal accuracy. Although the second-order implicit exponential time differencing (iETD2) method can accurately handle explicitly time-dependent reactions, it is more computationally expensive than IIF2. In this paper, we develop a hybrid approach that combines the advantages of both methods, applying IIF2 to reaction terms that are not explicitly time-dependent and applying iETD2 to those which are. The second-order $\underline {\text{h}} {\text{ybrid}}$ ${\text{I}}\underline {{\text{IF}}} - \underline {\text{E}} {\text{TD}}$ method (hIFE2) inherits the lower complexity of IIF2 and the ability to remain second-order accurate in time for large time-steps from iETD2. Also, it inherits the unconditional stability from IIF2 and iETD2 methods for dealing with the stiffness in reaction–diffusion systems. Through a transformation, hIFE2 can handle nonhomogeneous boundary conditions accurately and efficiently. In addition, this approach can be naturally combined with the compact and array representations of IIF and ETD for systems in higher spatial dimensions. Various numerical simulations containing linear and nonlinear reactions are presented to demonstrate the superior stability, accuracy, and efficiency of the new hIFE method.

Keywords:

Mathematics Subject Classification: Primary: 65M06, 35K57; Secondary: 65M12.

Citation:

Full Text(HTML)

Figure 1. Plots of the numerical error at $ T = 1 $ after applying IIF2, iETD2, and hIFE2 to the scalar equation in (8) with $ u(0) = 1 $ for various $ \Delta t $. Plots are shown for (A) $ f(u,t) = t^2 $ with $ \alpha = -10^1,\ -10^2,\ -10^3,\ -10^4,\ -10^5 $, and $ -10^6 $; (B) $ f(u,t) = -u $ with $ \alpha = -8,\ -16,\ -32,\ -64 $, and $ -128 $; and (C) $ f(u,t) = -u+t^2 $ with $ \alpha = -10^2,\ -10^3,\ -10^4,\ -10^5 $, and $ -10^6 $. The curves for iETD2 and hIFE2 are identical in (A), and those for IIF2 and hIFE2 are identical in (B). We see that for the time-dependent reactions (A, C), the error in IIF2 increases as $ -\alpha $ increases while the error in iETD2 and hIFE2 decreases

Download: Full-size image PowerPoint slide

Figure 2. The temporal errors at $ T = 1 $ in the maximum norm when solving the semi-discrete form (16) of (27) for different reactions with the IIF, iETD2, and hIFE2 methods. In all simulations, the reaction coefficient $ d = 1 $. (A) IIF2 for $ F(U,t) = t^2 $; (B) iETD2 for $ F(U,t) = t^2 $; (C) hIFE2 for $ F(U,t) = t^2 $; (D) IIF2 for $ F(U,t) = -U $; (E) iETD2 for $ F(U,t) = -U $; (F) hIFE2 for $ F(U,t) = -U $; (G) IIF2 for $ F(U,t) = -U+t^2 $; (H) iETD2 for $ F(U,t) = -U+t^2 $; (I) hIFE2 for $ F(U,t) = -U+t^2 $. Different colors represent the number of points, $ N $, in the spatial discretization, where $ N = 32,\ 64,\ 128,\ 256,\ 512 $, and $ 1024 $. Subfigures in same row share the same $ y $-axis while subfigures in same column share the same $ x $-axis. Panels (B) and (C) are identical because hIFE2 treats time-dependent terms with iETD2, and panels (D) and (F) are identical since hIFE2 treats autonomous terms with IIF2

Download: Full-size image PowerPoint slide

Figure 3. Plots of the numerical error at $ T = 1 $ in maximum norm after applying hIFE2 to (27) with Neumann, Dirichlet, and mixed boundary conditions for various $ \Delta t $ and fixed $ N $. The hIFE2 is applied to both original and transformed (Section 3.2) equations. Plots are shown for hIFE2 on: (A) the original equation with Neumann boundary; (B) the original equation with Dirichlet boundary; (C) the original equation with mixed boundary; (D) the transformed equation with Neumann boundary; (E) the transformed equation with Dirichlet boundary; (F) the transformed equation with mixed boundary. Different colors represent different spatial mesh sizes $ N $, where $ N = 32,\ 64,\ 128,\ 256,\ 512 $, and $ 1024 $

Download: Full-size image PowerPoint slide

Table 1. The truncation errors of IIF2, iETD2, and hIFE2 when applied to (8) with different reactions

Method	Reaction $ f $	Truncation error
IIF2	$ g(t) $	$ -\frac{1}{12}\Delta t^2(\alpha ^2g_n-2\alpha g'_n+g''_n) $
	$ g(t) $	$ -\frac{1}{24}\Delta t^3(\alpha ^3g_n-\alpha ^2 g'_n-\alpha g''_n+g'''_n)+\mathcal{O}(\Delta t^4) $
	$ ru $	$ -\frac{1}{12}\Delta t^2r^3u_n $
	$ ru $	$ -\frac{1}{24}\Delta t^3(2\alpha r^3+r^4)u_n+\mathcal{O}(\Delta t^4) $
	$ ru+g(t) $	$ -\frac{1}{12}\Delta t^2 \big[ \alpha^2g_n+\alpha(-rg_n-2g_n')+(r^3u_n+r^2g_n+rg_n'+g_n'') \big] $
		$ -\frac{1}{24}\Delta t^3\big[ \alpha^3g_n+\alpha^2(-rg_n-g_n')+\alpha(2r^3u_n+r^2g_n-g_n'') $
		$ +(r^4u_n+r^3g_n+r^2g_n'+rg_n''+g_n''')\big]+\mathcal{O}(\Delta t^4) $
iETD2	$ g(t) $	$ -\frac{1}{12}\Delta t^2g''_n $
	$ g(t) $	$ -\frac{1}{24}\Delta t^3(\alpha g''_n+g'''_n)+\mathcal{O}(\Delta t^4) $
	$ ru $	$ -\frac{1}{12}\Delta t^2(\alpha^2r+2\alpha r^2+r^3)u_n $
	$ ru $	$ -\frac{1}{24}\Delta t^3(2\alpha^3 r+5\alpha^2 r^2+4\alpha r^3+r^4)u_n+\mathcal{O}(\Delta t^4) $
	$ ru+g(t) $	$ -\frac{1}{12}\Delta t^2 \big[\alpha^2ru_n+\alpha(2r^2u_n+rg_n)+(r^3u_n+r^2g_n+rg_n'+g_n'')\big] $
		$ -\frac{1}{24}\Delta t^3\big[2\alpha^3ru_n-\alpha^2(5r^2u_n+2rg_n)+\alpha(4r^3u_n+3r^2g_n+2rg_n'+g_n'') $
		$ +(r^4u_n+r^3g_n+r^2g_n'+rg_n''+g_n''')\big]+\mathcal{O}(\Delta t^4) $
hIFE2	$ g(t) $	equivalent to iETD2
	$ ru $	equivalent to IIF2
	$ ru+g(t) $	$ -\frac{1}{12}\Delta t^2 \big[-\alpha rg_n+(r^3u_n+r^2g_n+rg_n'+g_n'')\big] $
		$ -\frac{1}{24}\Delta t^3\big[-\alpha^2rg_n+\alpha(2r^3u_n+r^2g_n+g_n'') $
		$ +(r^4u_n+r^3g_n+r^2g_n'+rg_n''+g_n''')\big]+\mathcal{O}(\Delta t^4) $

| Show Table

DownLoad: CSV

Table 2. Eigenvalues of $ A $, $ \lambda_j $, under different spatial resolutions, where $ d = 1 $, $ a = 0 $, $ b = \pi/2 $, $ j = 1,\ 5,\ N/2,\ N $

	32	64	128	256	512	1024
1	$ - $1.00	$ - $1.00	$ - $1.00	$ - $1.00	$ - $1.00	$ - $1.00
5	$ - $7.97e+1	$ - $8.07e+1	$ - $8.09e+1	$ - $8.09e+1	$ - $8.10e+1	$ - $8.10e+1
$ N/2 $	$ - $7.89e+2	$ - $3.23e+03	$ - $1.31e+04	$ - $5.28e+04	$ - $2.12e+05	$ - $8.49e+05
$ N $	$ - $1.66e+03	$ - $6.64e+03	$ - $2.66e+04	$ - $1.06e+05	$ - $4.25e+05	$ - $1.70e+06

| Show Table

DownLoad: CSV

Table 3. Different boundary conditions in (27), and their corresponding $ A $ and $ B(t) $ in the semi-discrete form (26)

	Neumann	Dirichlet	Mixed
BCs	$ u_x\vert_{x=0}=e^{-2t}\cos \frac{\pi}{6} $	$ u\vert_{x=0}=e^{-2t}\sin \frac{\pi}{6}, $	$ u_x\vert_{x=0}=e^{-2t}\cos \frac{\pi}{6} $
BCs	$ u_x\vert_{x=\frac{\pi}{2}}=e^{-2t}\cos \frac{2\pi}{3} $	$ u\vert_{x=\frac{\pi}{2}}=e^{-2t}\sin \frac{2\pi}{3} $	$ u\vert_{x=\frac{\pi}{2}}=e^{-2t}\sin \frac{2\pi}{3} $
$ B(t) $	$e^{-2t}{ \begin{bmatrix}-\frac{2\cos \frac{\pi}{6}}{\Delta x}\\0\\ \vdots \\0 \\\frac{2\cos \frac{2\pi}{3}}{\Delta x}\end{bmatrix}_{N+1}}$	$e^{-2t}{\begin{bmatrix}\frac{\sin \frac{\pi}{6}}{\Delta x^2}\\ 0\\ \vdots \\0 \\ \frac{\sin \frac{2\pi}{3}}{\Delta x^2}\end{bmatrix}_{N-1}}$	$e^{-2t}{\begin{bmatrix}-\frac{2\cos \frac{\pi}{6}}{\Delta x} \\0\\ \vdots \\0\\ \frac{\sin \frac{2\pi}{3}}{\Delta x^2}\end{bmatrix}_N}$
$ A $	$\frac{1}{\Delta x^2} \begin{bmatrix} -2&2&&\\ 1&-2&1\\ &\ddots&\ddots&\ddots\\ &&1&-2&1\\ &&&2&-2\\ \end{bmatrix}_{(N+1)^2}$	$\frac{1}{\Delta x^2} \begin{bmatrix} -2&1&&\\ 1&-2&1\\ &\ddots&\ddots&\ddots\\ &&1&-2&1\\ &&&1&-2\\ \end{bmatrix}_{(N-1)^2}$	$\frac{1}{\Delta x^2} \begin{bmatrix} -2&2&&\\ 1&-2&1\\ &\ddots&\ddots&\ddots\\ &&1&-2&1\\ &&&1&-2\\ \end{bmatrix}_{N^2}$

| Show Table

DownLoad: CSV

Table 4. Numerical errors in terms of the maximum norm and CPU time for the various methods on the example in Section 5.1 at $ T = 1 $ with diffusion coefficient $ d = 2 $. Here $ N $ is the number of grid points in the spatial discretization ($ \Delta x = \pi/2N $), and the time step $ \Delta t = 0.1\Delta x $. "CPU time 1" is the CPU time for initializing the matrices (Appendix C), "CPU time 2" is the CPU time for the iterations, and "CPU time" is the sum of the two

	$ N $	$ L^{\infty} $ error	Order	CPU time (s)	CPU time 1 (s)	CPU time 2 (s)
IIF2	8	0.00228	-	0.09	0.05	0.04
	16	0.000591	1.95	0.04	0.02	0.02
	32	0.000198	1.58	0.07	0.03	0.04
	64	7.81e-05	1.34	0.13	0.04	0.09
	128	0.000108	$ - $0.46	0.54	0.07	0.47
	256	5.18e-05	1.06	1.26	0.23	1.03
	512	1.83e-05	1.50	4.00	1.39	2.61
	1024	2.07e-05	$ - $0.18	28.30	7.75	20.55
	2048	1.07e-05	0.96	168.12	42.10	126.02
	4096	5.35e-06	1.00	1148.42	265.35	883.07
	$ N $	$ L^{\infty} $ error	Order	CPU time (s)	CPU time 1 (s)	CPU time 2 (s)
iETD2	8	0.00216	-	0.07	0.04	0.03
	16	0.000539	2.00	0.07	0.04	0.03
	32	0.000135	2.00	0.12	0.06	0.06
	64	3.37e-05	2.00	0.80	0.07	0.73
	128	8.41e-06	2.00	3.78	0.16	3.62
	256	2.1e-06	2.00	22.99	0.54	22.45
	512	5.26e-07	2.00	289.66	2.70	286.96
	1024	1.32e-07	2.00	2841.66	14.65	2827.01
	2048	3.31e-08	1.99	35348.32	91.84	35256.48
	4096	-	-	$ \text{too long} $	-	-
	$ N $	$ L^{\infty} $ error	Order	CPU time (s)	CPU time 1 (s)	CPU time 2 (s)
hIFE2	8	0.00217	-	0.12	0.09	0.03
	16	0.000544	1.99	0.06	0.04	0.02
	32	0.000137	1.99	0.08	0.05	0.03
	64	3.42e-05	2.00	0.16	0.08	0.08
	128	8.75e-06	1.97	0.76	0.17	0.59
	256	2.21e-06	1.99	1.85	0.54	1.31
	512	5.53e-07	2.00	9.17	2.61	6.56
	1024	1.49e-07	1.89	61.82	14.20	47.62
	2048	3.93e-08	1.93	419.24	89.49	329.75
	4096	1.12e-08	1.81	3096.23	603.04	2493.19
	$ N $	$ L^{\infty} $ error	Order	CPU time (s)	CPU time 1 (s)	CPU time 2 (s)
fEIF2	8	0.00216	-	0.37	0.37	0.00
	16	0.00054	2.00	0.04	0.04	0.00
	32	0.000135	2.00	0.07	0.07	0.00
	64	3.38e-05	2.00	0.09	0.08	0.01
	128	8.44e-06	2.00	0.54	0.18	0.36
	256	2.11e-06	2.00	1.41	0.69	0.72
	512	5.28e-07	2.00	11.62	3.01	8.61
	1024	1.32e-07	2.00	84.11	16.11	68.00
	2048	3.31e-08	1.99	613.91	101.12	512.79
	4096	8.89e-09	1.90	4700.11	707.64	3992.47

| Show Table

DownLoad: CSV

Table 5. Numerical errors and CPU time for the test in Section 5.2 at time $ T = 1 $. We set the diffusion coefficient $ d = 0.1 $ and the coefficients of the reactions $ a = 500 $ and $ b = -2 $. For each simulation, we fix the number of grid points $ N = 1024 $ ($ \Delta x = \pi/2N $), and run the simulation for $ K $ time steps ($ \Delta t = T/K $). The error $ e $ is measured in the maximum norm, and the relative error is defined by $ e / \max\{ \| U_K\|_{\infty},\| V_K\|_{\infty} \} $, where $ U_K $ and $ V_K $ are the numerical solutions after $ K $ time steps. "CPU time 1" is the CPU time for initialization (Appendix C), "CPU time 2" is the CPU time for the iterations, and "CPU time" is the sum of the two

	$ K $	$ L^{\infty} $ error	Relative error	Order	CPU time (s)	CPU time 1 (s)	CPU time 2 (s)
IIF2	20	10	0.00381	-	5.32	5.25	0.07
	40	4.81	0.00182	1.06	5.07	4.91	0.16
	80	2.32	0.000881	1.05	5.09	4.78	0.31
	160	1.12	0.000425	1.05	5.07	4.44	0.63
	320	0.534	0.000203	1.07	5.24	3.90	1.34
	640	0.251	9.51e-05	1.09	5.92	3.40	2.52
	1280	0.115	4.34e-05	1.13	7.90	2.92	4.98
	2560	0.0503	1.91e-05	1.19	12.84	2.55	10.29
	$ K $	$ L^{\infty} $ error	Relative error	Order	CPU time (s)	CPU time 1 (s)	CPU time 2 (s)
iETD2	20	3.99	0.00151	-	19.63	10.88	8.75
	40	0.994	0.000377	2.00	28.60	10.80	17.80
	80	0.248	9.41e-05	2.00	46.92	10.76	36.16
	160	0.0617	2.34e-05	2.01	80.10	10.41	69.69
	320	0.0152	5.76e-06	2.02	148.60	9.80	138.80
	640	0.00366	1.39e-06	2.05	285.20	9.27	275.93
	1280	0.000872	3.31e-07	2.07	567.11	8.94	558.17
	2560	0.000227	8.61e-08	1.94	1140.59	8.49	1132.10
	$ K $	$ L^{\infty} $ error	Relative error	Order	CPU time (s)	CPU time 1 (s)	CPU time 2 (s)
hIFE2	20	4.19	0.000397	-	11.11	0.00	0.28
	40	1.05	0.000397	2.00	11.96	11.39	0.57
	80	0.261	9.91e-05	2.00	11.61	10.70	0.91
	160	0.0652	2.47e-05	2.00	12.36	10.37	1.99
	320	0.0162	6.14e-06	2.01	13.90	9.84	4.06
	640	0.00397	1.51e-06	2.03	17.65	9.43	8.22
	1280	0.000971	3.68e-07	2.03	25.08	8.88	16.20
	2560	0.000256	9.72e-08	1.92	40.83	8.45	32.38
	$ K $	$ L^{\infty} $ error	Relative error	Order	CPU time (s)	CPU time 1 (s)	CPU time 2 (s)
fEIF2	20	1.49e+29	4.43e+25	-	12.42	11.96	0.46
	40	2.9e+48	8.65e+44	$ - $64.08	12.61	11.78	0.83
	80	6.04e+73	1.8e+70	$ - $84.11	13.07	11.46	1.61
	160	2.27e+96	6.77e+92	$ - $74.99	14.43	11.20	3.23
	320	1.92e+79	5.71e+75	56.72	17.20	10.59	6.61
	640	0.251	7.48e-05	265.37	23.47	9.93	13.54
	1280	0.119	3.54e-05	1.08	35.96	9.57	26.39
	2560	0.0603	1.8e-05	0.98	62.05	9.08	52.97

| Show Table

DownLoad: CSV

Table 6. Numerical a priori error in applying hIFE2 to a one-dimensional reaction–diffusion system with (A) $ f(u,x,t) = \cos u+t $ for the decomposition (48) and (49) and (B) $ f(u,x,t) = (t+1)\cos (xu)+xe^t $ for the decomposition (51) and (52). The a priori error is defined by $ \|u^N-u^{N/2}\|_{\infty} $, where $ N $ is the number of grid points in the spatial discretization. The simulations are run through time $ T = 1 $ with $ \Delta x = \frac{\pi}{2N} $ and $ \Delta t = 0.1\Delta x $

(A)		Decomposition (48)		Decomposition (49)
	$ N $	A priori error	Order	A priori error	Order
	16	0.00103	-	0.00102	-
	32	0.000532	0.95	0.000255	2.00
	64	0.000328	0.70	6.37e-05	2.00
	128	8.21e-05	2.00	1.59e-05	2.00
	256	0.000196	$ - $1.25	3.98e-06	2.00
	512	0.000106	0.88	9.95e-07	2.00
	1024	8.69e-06	3.61	2.49e-07	2.00
	2048	3.37e-05	$ - $1.96	6.19e-08	2.01
	4096	1.75e-05	0.95	1.42e-08	2.12
(B)		Decomposition (51)		Decomposition (52)
	N	A priori error	Order	A priori error	Order
	16	0.00979	-	0.00977	-
	32	0.00245	2.00	0.00245	1.99
	64	6.67e-04	1.88	0.000614	2.00
	128	1.67e-04	2.00	0.000153	2.00
	256	3.94e-04	-1.24	3.83e-05	2.00
	512	2.13e-04	0.89	9.59e-06	2.00
	1024	1.76e-05	3.60	2.4e-06	2.00
	2048	6.75e-05	-1.94	5.99e-07	2.00
	4096	3.50e-05	0.95	1.48e-07	2.02

| Show Table

DownLoad: CSV

Table 7. Numerical errors in the maximum norm for hIFE2 applied to the example in Section 5.4. The spatial resolution is $ \Delta x = \frac{\pi}{2N} $ in all three dimensions, the time step is $ \Delta t = 0.1\Delta x $, the ending time is $ T = 1 $, and the coefficients are $ d_1 = d_2 = d_3 = 1 $ and $ r = -1 $

$ N\times N\times N $	$ L^{\infty} $ error	Order
$ 4\times 4 \times 4 $	1.33e-03	-
$ 8\times 8 \times 8 $	3.28e-04	2.02
$ 16\times 16 \times 16 $	8.17e-05	2.01
$ 32\times 32\times 32 $	2.04e-05	2.00
$ 64\times 64 \times 64 $	5.10e-05	2.00
$ 128\times 128 \times 128 $	1.27e-06	2.00

| Show Table

DownLoad: CSV

Table 8. A summary of the four methods: for their A-stability, and the restriction on $ \Delta t $ to exhibit second order, with explicitly time-dependent reactions or nonhomogeneous boundary conditions

Method	A-stability	$ \Delta t $ to exhibit second-order accuracy
Method	A-stability	Time-dependent reactions	Nonhomogeneous BCs
IIF2	Yes	$ \mathcal{O}(\Delta x^2) $	$ \leq \mathcal{O}(\Delta x^2) $
iETD2	Yes	$ \mathcal{O}(1) $	-
fEIF2	No	$ \mathcal{O}(1) $	$ \mathcal{O}(1) $
hIFE2	Yes	$ \mathcal{O}(1) $	$ \mathcal{O}(1) $
hIFE2 (transformed)	Yes	$ \mathcal{O}(1) $	$ \mathcal{O}(1) $

| Show Table

DownLoad: CSV

Table 9. A comparison of the computational complexity between the IIF2, iETD2, hIFE2, and fEIF2 methods

	Operations per iteration	Total complexity (ratio)
IIF2	$\mathcal{O}(N^2)$	1
iETD2	$\mathcal{O}(kN^3)$	$\mathcal{O}(kN)$
hIFE2	$\mathcal{O}(3N^2)$	3
fEIF2	$\mathcal{O}(5N^2)$	5

| Show Table

DownLoad: CSV

Related Papers

Cited by

References

[1]	D. Alonso, F. Bartumeus and J. Catalan, Mutual interference between predators can give rise to turing spatial patterns, Ecology, 83 (2002), 28-34.
[2]	G. Beylkin, J. M. Keiser and L. Vozovoi, A new class of time discretization schemes for the solution of nonlinear PDEs, Journal of Computational Physics, 147 (1998), 362-387. doi: 10.1006/jcph.1998.6093.
[3]	S. Chen and Y.-T. Zhang, Krylov implicit integration factor methods for spatial discretization on high dimensional unstructured meshes: application to discontinuous Galerkin methods, Journal of Computational Physics, 230 (2011), 4336-4352. doi: 10.1016/j.jcp.2011.01.010.
[4]	S. M. Cox and P. C. Matthews, Exponential time differencing for stiff systems, Journal of Computational Physics, 176 (2002), 430-455. doi: 10.1006/jcph.2002.6995.
[5]	P. D. Dale, J. A. Sherratt and P. K. Maini, Role of fibroblast migration in collagen fiber formation during fetal and adult dermal wound healing, Bulletin of mathematical biology, 59 (1997), 1077-1100. doi: 10.1007/BF02460102.
[6]	Q. Du and W. Zhu, Stability analysis and application of the exponential time differencing schemes, Journal of Computational Mathematics, 22 (2014), 200-209.
[7]	Q. Du and W. Zhu, Analysis and applications of the exponential time differencing schemes and their contour integration modifications, BIT Numerical Mathematics, 45 (2005), 307-328. doi: 10.1007/s10543-005-7141-8.
[8]	A. Eldar, R. Dorfman, D. Weiss, H. Ashe, B.-Z. Shilo and N. Barkai, Robustness of the bmp morphogen gradient in drosophila embryonic patterning, Nature, 419 (2002), 304-308. doi: 10.1038/nature01061.
[9]	E. Gallopoulos and Y. Saad, Efficient solution of parabolic equations by Krylov approximation methods, SIAM Journal on Scientific and Statistical Computing, 13 (1992), 1236-1264. doi: 10.1137/0913071.
[10]	M. Garvie and C. Trenchea, Analysis of two generic spatially extended predator-prey models, Nonlinear Anal. Real World Appl.
[11]	M. R. Garvie, Finite-difference schemes for reaction–diffusion equations modeling predator–prey interactions in matlab, Bulletin of mathematical biology, 69 (2007), 931-956. doi: 10.1007/s11538-006-9062-3.
[12]	A. Gierer and H. Meinhardt, A theory of biological pattern formation, Biological Cybernetics, 12 (1972), 30-39. doi: 10.1007/BF00289234.
[13]	M. Hochbruck and C. Lubich, On Krylov subspace approximations to the matrix exponential operator, SIAM Journal on Numerical Analysis, 34 (1997), 1911-1925. doi: 10.1137/S0036142995280572.
[14]	M. Hochbruck and A. Ostermann, Explicit exponential Runge–Kutta methods for semilinear parabolic problems, SIAM Journal on Numerical Analysis, 43 (2005), 1069-1090. doi: 10.1137/040611434.
[15]	M. Hochbruck and A. Ostermann, Exponential Runge–Kutta methods for parabolic problems, Applied Numerical Mathematics, 53 (2005), 323-339. doi: 10.1016/j.apnum.2004.08.005.
[16]	M. Hochbruck and A. Ostermann, Exponential integrators, Acta Numerica, 19 (2010), 209-286. doi: 10.1017/S0962492910000048.
[17]	E. E. Holmes, M. A. Lewis, J. Banks and R. Veit, Partial differential equations in ecology: spatial interactions and population dynamics, Ecology, 75 (1994), 17-29. doi: 10.2307/1939378.
[18]	E. Isaacson and H. B. Keller, Analysis of Numerical Methods, John Wiley & Sons, Inc., New York-London-Sydney, 1966.
[19]	T. Jiang and Y.-T. Zhang, Krylov implicit integration factor WENO methods for semilinear and fully nonlinear advection–diffusion–reaction equations, Journal of Computational Physics, 253 (2013), 368-388. doi: 10.1016/j.jcp.2013.07.015.
[20]	T. Jiang and Y.-T. Zhang, Krylov single-step implicit integration factor WENO methods for advection–diffusion–reaction equations, Journal of Computational Physics, 311 (2016), 22-44. doi: 10.1016/j.jcp.2016.01.021.
[21]	L. Ju, X. Liu and W. Leng, Compact implicit integration factor methods for a family of semilinear fourth-order parabolic equations, Discrete & Continuous Dynamical Systems-Series B, 19 (2014), 1667-1687. doi: 10.3934/dcdsb.2014.19.1667.
[22]	L. Ju, J. Zhang, L. Zhu and Q. Du, Fast explicit integration factor methods for semilinear parabolic equations, Journal of Scientific Computing, 62 (2015), 431-455. doi: 10.1007/s10915-014-9862-9.
[23]	A.-K. Kassam and L. N. Trefethen, Fourth-order time-stepping for stiff PDEs, SIAM Journal on Scientific Computing, 26 (2005), 1214-1233. doi: 10.1137/S1064827502410633.
[24]	A. Kicheva, P. Pantazis, T. Bollenbach, Y. Kalaidzidis, T. Bittig, F. Jülicher and M. Gonzalez-Gaitan, Kinetics of morphogen gradient formation, Science, 315 (2007), 521-525. doi: 10.1126/science.1135774.
[25]	J. D. Lambert, Numerical Methods for Ordinary Differential Systems: The Initial Value Problem, John Wiley & Sons, Inc., 1991.
[26]	A. D. Lander, Q. Nie and F. Y. Wan, Do morphogen gradients arise by diffusion?, Developmental cell, 2 (2002), 785-796.
[27]	S. Larsson and V. Thomée, Partial Differential Equations with Numerical Methods, vol. 45, Springer-Verlag, Berlin, 2003.
[28]	X. Liu and Q. Nie, Compact integration factor methods for complex domains and adaptive mesh refinement, Journal of Computational Physics, 229 (2010), 5692-5706. doi: 10.1016/j.jcp.2010.04.003.
[29]	D. Lu and Y.-T. Zhang, Krylov integration factor method on sparse grids for high spatial dimension convection–diffusion equations, Journal of Scientific Computing, 69 (2016), 736-763. doi: 10.1007/s10915-016-0216-7.
[30]	D. Lu and Y.-T. Zhang, Computational complexity study on Krylov integration factor WENO method for high spatial dimension convection–diffusion problems, Journal of Scientific Computing, 73 (2017), 980-1027. doi: 10.1007/s10915-017-0398-7.
[31]	M. Machen and Y.-T. Zhang, Krylov implicit integration factor methods for semilinear fourth-order equations, Mathematics, 5 (2017), 63. doi: 10.3390/math5040063.
[32]	R. E. Mickens, Nonstandard finite difference schemes for reaction-diffusion equations, Numerical Methods for Partial Differential Equations, 15 (1999), 201-214. doi: 10.1002/(SICI)1098-2426(199903)15:2<201::AID-NUM5>3.0.CO;2-H.
[33]	Q. Nie, F. Y. Wan, Y.-T. Zhang and X. Liu, Compact integration factor methods in high spatial dimensions, Journal of Computational Physics, 227 (2008), 5238-5255. doi: 10.1016/j.jcp.2008.01.050.
[34]	Q. Nie, Y.-T. Zhang and R. Zhao, Efficient semi-implicit schemes for stiff systems, Journal of Computational Physics, 214 (2006), 521-537. doi: 10.1016/j.jcp.2005.09.030.
[35]	S. V. Petrovskii and H. Malchow, A minimal model of pattern formation in a prey-predator system, Mathematical and Computer Modelling, 29 (1999), 49-63. doi: 10.1016/S0895-7177(99)00070-9.
[36]	S. V. Petrovskii and H. Malchow, Wave of chaos: new mechanism of pattern formation in spatio-temporal population dynamics, Theoretical population biology, 59 (2001), 157-174. doi: 10.1006/tpbi.2000.1509.
[37]	H. Risken, The Fokker-Planck Equation. Methods of Solution and Applications, Springer Series in Synergetics, 18. Springer-Verlag, Berlin, 1984. doi: 10.1007/978-3-642-96807-5.
[38]	Y. Saad, Analysis of some Krylov subspace approximations to the matrix exponential operator, SIAM Journal on Numerical Analysis, 29 (1992), 209-228. doi: 10.1137/0729014.
[39]	J. C. Schulze, P. J. Schmid and J. L. Sesterhenn, Exponential time integration using Krylov subspaces, International Journal for Numerical Methods in Fluids, 60 (2009), 591-609. doi: 10.1002/fld.1902.
[40]	C. Ta, D. Wang and Q. Nie, An integration factor method for stochastic and stiff reaction–diffusion systems, Journal of Computational Physics, 295 (2015), 505-522. doi: 10.1016/j.jcp.2015.04.028.
[41]	J. W. Thomas, Numerical Partial Differential Equations: Finite Difference Methods, Texts in Applied Mathematics, 33. Springer-Verlag, New York, 1999. doi: 10.1007/978-1-4612-0569-2.
[42]	A. M. Turing, The chemical basis of morphogenesis, Bulletin of mathematical biology, 52 (1990), 153-197.
[43]	C. Van Loan, Computational Frameworks for the Fast FOurier Transform, vol. 10, SIAM, 1992. doi: 10.1137/1.9781611970999.
[44]	D. Wang, W. Chen and Q. Nie, Semi-implicit integration factor methods on sparse grids for high-dimensional systems, Journal of Computational Physics, 292 (2015), 43-55. doi: 10.1016/j.jcp.2015.03.033.
[45]	D. Wang, L. Zhang and Q. Nie, Array-representation integration factor method for high-dimensional systems, Journal of Computational Physics, 258 (2014), 585-600. doi: 10.1016/j.jcp.2013.11.002.
[46]	O. Wartlick, A. Kicheva and M. González-Gaitán, Morphogen gradient formation, Cold Spring Harbor perspectives in biology, 1 (2009), a001255. doi: 10.1101/cshperspect.a001255.
[47]	A.-M. Wazwaz, Partial Differential Equations, CRC Press, 2002.
[48]	A. Wiegmann, Fast Poisson, fast Helmholtz and fast linear elastostatic solvers on rectangular parallelepipeds, Lawrence Berkeley National Laboratory. doi: 10.2172/982430.
[49]	S. R. Yu, M. Burkhardt, M. Nowak, J. Ries, Z. Petrášek, S. Scholpp, P. Schwille and M. Brand, Fgf8 morphogen gradient forms by a source-sink mechanism with freely diffusing molecules, Nature, 461 (2009), 533-536. doi: 10.1038/nature08391.
[50]	L. Zhang, A. D. Lander and Q. Nie, A reaction–diffusion mechanism influences cell lineage progression as a basis for formation, regeneration, and stability of intestinal crypts, BMC Systems Biology, 6 (2012), 93. doi: 10.1186/1752-0509-6-93.
[51]	S. Zhao, J. Ovadia, X. Liu, Y.-T. Zhang and Q. Nie, Operator splitting implicit integration factor methods for stiff reaction–diffusion–advection systems, Journal of Computational Physics, 230 (2011), 5996-6009. doi: 10.1016/j.jcp.2011.04.009.
[52]	L. Zhu, L. Ju and W. Zhao, Fast high-order compact exponential time differencing runge–kutta methods for second-order semilinear parabolic equations, Journal of Scientific Computing, 67 (2016), 1043-1065. doi: 10.1007/s10915-015-0117-1.
[53]	Y.-L. Zhu, X. Wu, I.-L. Chern and Z.-Z. Sun, Derivative Securities and Difference Methods, Springer, 2004. doi: 10.1007/978-1-4757-3938-1.