Research Article  Open Access
On the Long Time Simulation of ReactionDiffusion Equations with Delay
Abstract
For a consistent numerical method to be practically useful, it is widely accepted that it must preserve the asymptotic stability of the original continuous problem. However, in this study, we show that it may lead to unreliable numerical solutions in long time simulation even if a classical numerical method gives a larger stability region than that of the original continuous problem. Some numerical experiments on the reactiondiffusion equations with delay are presented to confirm our findings. Finally, some open problems on the subject are proposed.
1. Introduction
Reactiondiffusion equations with delay have attracted a significant interest in the last several decades due to their frequent appearance in actual applications [1–11]. They are usually proposed as models for neural networks, population ecology, and cell biology. As far as we know, a delay term could not only make the exact solutions difficult to be obtained in general but also change the long time dynamic properties of a system. That is the reason why such equations became the focus of researchers in numerical analysis and simulation. Up to now, most of the numerical methods available to approximate the diffusion term are based on some classical numerical methods, for example, central finite difference method [3–7], compact finite difference method [8, 9], discontinuous Galerkin methods [10], spectral collocation method [11], and so forth. Many important properties of these numerical methods such as asymptotic stability and convergence have been extensively investigated. However, little research has been conducted in studying the long time dynamical properties of these numerical methods.
The long time dynamical properties of different numerical methods for PDEs are attractive (e.g., [12, 13]). As a typical example in the delay field, we refer the readers to the works [14, 15], where the authors considered the following reactiondiffusion equation with delay:
As it is pointed out in the papers [14, 15], the spatial discretization based on standard central difference method (CDM) unavoidably leads to a reduction of the asymptotic stability region of the corresponding PDE problem. In other words, the long time numerical solutions derived by CDM may be unstable even if the original problem (1) is stable. As a result, it is widely accepted that a numerical method should preserve stability of the original continuous problems in actual application (e.g., [14–20]). This naturally leads to the following question: are the long time numerical solutions reliable if the discrete problem has a larger stability region than that of the continuous one?
In this study, it is shown that a negative answer to the above question can be obtained. The conclusion is obtained when we try to find some highorder numerical methods, which are sufficient to preserve stability of the continuous problem. It is well known that the finite element method (FEM) can be easily designed for high order of accuracy in space. We have proved that onedegree FEM gives a larger stability region than that of the continuous problem. It means that the numerical method is sufficient to preserve the asymptotic stability of the original problems. Numerical experiments indicate that the long time numerical solutions decay in time. However, after some more numerical tests, we find that the long time numerical solutions derived by FEM may not be reliable for the following two reasons. First, the numerical solutions derived by the FEM decay more quickly than the ones by the alternative central difference method (ACDM). The latter method is proved to be a good choice for solving the problem (1) in [15]. Second, if we reduce the stepsize, the long time FEM solutions at every fixed time have changed, while the solutions derived by ACDM nearly remain the same.
The rest of the paper is organized as follows. Section 2 describes in detail the FEM and stability analysis of the problem. Section 3 shows experimental studies for verifying the proposed results. Finally, conclusions and discussions for this paper are summarized in Section 4.
2. Main Results
In this section, after a statement of the classical FEM, stability of the discrete reaction diffusion equations with delay is investigated.
For the spatial discretization of system (1), we divide the interval with a mesh: with the space stepsize . Let be a finite dimensional space; then the semidiscrete FEM is to find , such that for all test functions and , satisfying
For the classical onedegree FEM, the piecewise linear function is chosen as follows: , , and
Let ; one can derive the following delay differential equations in the matrix form:
Here with denoting an approximation to , , and , where with its eigenvalues , .
In what follows, we will turn our attention to the stability analysis of the semidiscrete PDE (4). First, we introduce some lemmas, which will assist in the proof of our main result.
Lemma 1. Suppose that is an eigenvector of the matrix with the corresponding eigenvalue . Then, is also an eigenvector of the matrix with the corresponding eigenvalue , where matrixes and are defined in (4). Moreover, the maximum eigenvalue of the matrix is given by
Proof. It follows from that
Hence,
Therefore,
Meanwhile, noting that , we have
When , we have the maximum denominator and the minimum numerator . Thus, the maximum eigenvalue is
Moreover, substituting the Maclaurin series expansions
into (11) yields
which completes the proof.
Lemma 2 (cf. [14]). Suppose that all roots of the following equation: have negative real part for any given , . Then for any , all the roots of the equation have negative real part too.
The following lemma can be derived from [14, Corollary 2.12] or [16, Proposition 1].
Lemma 3. The following equation: has negative real part if and only if and , where is the root of such that .
Theorem 4. Assume . Then the zero solution of the semidiscrete delay PDE (2) is asymptotically stable if and only if , where is defined in (6) and is the root of such that .
Proof. The zero solution to semidiscrete delay PDE (4) is asymptotically stable if all the roots of the corresponding characteristic equation
have negative real parts.
According to Lemma 1, the eigenvalues of the matrix are . Therefore, the asymptotic stability of system is equivalent to the condition that all the roots of each of the equations
have the negative real part.
Together with Lemma 2 and the maximum eigenvalue (it is proved in Lemma 1), we claim that the asymptotic stability of (2) is equivalent to the condition that the root of the following equation:
has the negative real part.
Finally, application of Lemma 3 gives the result.
The asymptotical stability of the original PDE (1) is shown as follows.
Lemma 5 (cf. [15]). The zero solution of (1) is asymptotically stable if and only if and , where is the root of such that .
The preceding theorem and Lemma 5 immediately yield the following result.
Theorem 6. The spatial discretization based on standard FEM gives a larger asymptotic stability region than that of the corresponding PDE.
Proof. Together with (it is proved in Lemma 1), Lemma 5, and Theorem 4, we have the conclusion.
Theorem 6 implies that the FEM, combined with an appropriate time discretization, can be sufficient to preserve the asymptotic stability of the continuous delay PDEs (1). However, in the next section, we will show that onedegree FEM may lead to unreliable long time numerical solutions.
3. Numerical Experiments
In this section, we will present several numerical experiments to illustrate our findings, where the diffusion term is discreted by the FEM, the CDM and the ACDM (cf. [14, 15]) respectively, where is a numerical approximation to . Moreover, the time discretization is done by the trapezoidal rule, which is proved as a good choice for long term numerical simulation of in [15].
We firstly show that the stability region of the discrete delay PDE derived by FEM is larger than that of the original continuous problem. To do this, we set , , , the time interval [0,1200], and the following initial condition:
According to Lemma 1, the original problem is not asymptotically stable. Now, combined with the previously mentioned time discretization, we set space stepsize and time stepsize and apply the CDM, the ACDM, and the FEM to discrete the equation, respectively. The numerical solutions at different time are listed in Table 1. As it is shown, the numerical solutions derived by finite element method decay in time. It indicates that the stability region of discrete equation derived by FEM is larger than that of the original problem.

Next, we will show that the FEM solutions may not be reliable. Let , , and . It is easy to check that the original delay PDE is asymptotically stable. Now, we still set space stepsize and time stepsize and apply the different numerical method to solve problem (1) with the initial condition (22). Some statistics of the numerical results are presented in Table 2. Clearly, the numerical solutions derived by the CDM increase in time. And there is a decay of the numerical solutions derived by the other two methods. However, we find that the numerical solutions derived by the FEM decay more quickly than that by the ACDM. In order to distinguish which numerical solutions are reliable, we reduce the stepsize and let space stepsize and time stepsize . The numerical results with different methods can be found in Table 3. It is shown that the numerical solutions derived by ACDM remain the same, while the others have changed, which confirm our findings.


4. Conclusions and Discussions
It is pointed out in [14, 15] that the spatial discretization based on standard CDM unavoidably leads to a reduction of the asymptotic stability region of the corresponding continuous PDE problem. Therefore, for a numerical method to be practically useful for the long time simulation, it must preserve the asymptotic stability of the original continuous problem. It is reasonable but not enough. In this study, it is shown that the FEM is sufficient to preserve the asymptotically stability of the corresponding PDE problem, but the method may give long time unreliable solutions. The unreliable numerical solutions usually appear when the set of parameters is located near the boundary of the stability region of the continuous problem. As a result, the optimal numerical method for long time simulation may be the one whose stability region agrees with that of the original continuous PDEs.
The following problems on the subject remain to be considered in the future:(i)long time stability analysis of some classical numerical methods for solving the typical test model (1);(ii)construction of highorder numerical schemes, whose stability region agrees with that of the original continuous problems;(iii)extension of the long time stability analysis for the other partial differential equations with delay.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This work is supported by NSFC (Grant nos. 11201161, 11171125, and 91130003) and HSF (Grant no. 2011CDB289).
References
 J. Wu, Theory and Applications of Partial Functional Differential Equations, Springer, New York, NY, USA, 1996.
 D. Li, C. Zhang, and W. Wang, “Long time behavior of nonFickian delay reactiondiffusion equations,” Nonlinear Analysis—Real World Applications, vol. 13, no. 3, pp. 1401–1415, 2012. View at: Publisher Site  Google Scholar
 B. ZubikKowal, “Stability in the numerical solution of linear parabolic equations with a delay term,” BIT Numerical Mathematics, vol. 41, no. 1, pp. 191–206, 2001. View at: Google Scholar
 B. ZubikKowal and S. Vandewalle, “Waveform relaxation for functionaldifferential equations,” SIAM Journal on Scientific Computing, vol. 21, no. 1, pp. 207–226, 1999. View at: Google Scholar
 H. Tian, D. Zhang, and Y. Sun, “Delayindependent stability of Euler method for nonlinear onedimensional diffusion equation with constant delay,” Frontiers of Mathematics in China, vol. 4, no. 1, pp. 169–179, 2009. View at: Publisher Site  Google Scholar
 S. Wu and S. Gan, “Analytical and numerical stability of neutral delay integrodifferential equations and neutral delay partial differential equations,” Computers and Mathematics with Applications, vol. 55, no. 11, pp. 2426–2443, 2008. View at: Publisher Site  Google Scholar
 S.L. Wu, C.M. Huang, and T.Z. Huang, “Convergence analysis of the overlapping Schwarz waveform relaxation algorithm for reactiondiffusion equations with time delay,” IMA Journal of Numerical Analysis, vol. 32, no. 2, pp. 632–671, 2012. View at: Publisher Site  Google Scholar
 Z.Z. Sun and Z.B. Zhang, “A linearized compact difference scheme for a class of nonlinear delay partial differential equations,” Applied Mathematical Modelling, vol. 37, no. 3, pp. 742–752, 2013. View at: Publisher Site  Google Scholar
 Q. Zhang and C. Zhang, “A compact difference scheme combined with extrapolation techniques for solving a class of neutral delay parabolic differential equations,” Applied Mathematics Letters, vol. 26, no. 2, pp. 306–312, 2013. View at: Publisher Site  Google Scholar
 D. Li, C. Zhang, and H. Qin, “LDG method for reactiondiffusion dynamical systems with time delay,” Applied Mathematics and Computation, vol. 217, no. 22, pp. 9173–9181, 2011. View at: Publisher Site  Google Scholar
 Z. Jackiewicz and B. ZubikKowal, “Spectral collocation and waveform relaxation methods for nonlinear delay partial differential equations,” Applied Numerical Mathematics, vol. 56, no. 34, pp. 433–443, 2006. View at: Publisher Site  Google Scholar
 S. Gottlieb, F. Tone, C. Wang, X. Wang, and D. Wirosoetisno, “Long time stability of a classical effcient scheme for twodimensional Navierstokes equations,” SIAM Journal on Numerical Analysis, vol. 50, no. 1, pp. 126–150, 2012. View at: Google Scholar
 Z. Chen and X. Wu, “Longtime stability and convergence of the uniaxial perfectly matched layer method for timedomain acoustic scattering problems,” SIAM Journal on Numerical Analysis, vol. 50, no. 5, pp. 2632–2655, 2012. View at: Publisher Site  Google Scholar
 C. Huang and S. Vandewalle, “An analysis of delaydependent stability for ordinary and partial differential equations with fixed and distributed delays,” SIAM Journal on Scientific Computing, vol. 25, no. 5, pp. 1608–1632, 2004. View at: Publisher Site  Google Scholar
 C. Huang and S. Vandewalle, “Unconditionally stable difference methods for delay partial differential equations,” Numerische Mathematik, vol. 122, no. 3, pp. 579–601, 2012. View at: Publisher Site  Google Scholar
 C. Huang, “Delaydependent stability of high order RungeKutta methods,” Numerische Mathematik, vol. 111, no. 3, pp. 377–387, 2009. View at: Publisher Site  Google Scholar
 D. Li and C. Zhang, “Nonlinear stability of discontinuous Galerkin methods for delay differential equations,” Applied Mathematics Letters, vol. 23, no. 4, pp. 457–461, 2010. View at: Publisher Site  Google Scholar
 W. Wang and C. Zhang, “Preserving stability implicit Euler method for nonlinear Volterra and neutral functional differential equations in Banach space,” Numerische Mathematik, vol. 115, no. 3, pp. 451–474, 2010. View at: Publisher Site  Google Scholar
 W. Wang, C. Zhang, and D. Li, “Asymptotic stability of exact and discrete solutions for neutral multidelayintegrodifferential equations,” Applied Mathematical Modelling, vol. 35, no. 9, pp. 4490–4506, 2011. View at: Publisher Site  Google Scholar
 S. Diop, I. Kolmanovsky, P. E. Moraal, and M. van Nieuwstadt, “Preserving stability/performance when facing an unknown timedelay,” Control Engineering Practice, vol. 9, no. 12, pp. 1319–1325, 2001. View at: Publisher Site  Google Scholar
Copyright
Copyright © 2014 Dongfang Li and Chengjian Zhang. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.