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 reaction-diffusion equations with delay are presented to confirm our findings. Finally, some open problems on the subject are proposed.
1. Introduction
Reaction-diffusion 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 reaction-diffusion 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 high-order 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 one-degree 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 one-degree 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 one-degree 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 high-order 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).