Our main objective is to construct NSFD schemes preserving the positivity and
stability of (2.5.1). It is worth noting that the continuous model under consideration
possesses the equilibrium points which are not only LAS but also GAS. As mentioned
above, the construction of NSFD schemes preserving the GAS of continuous models
is a very important task and is also a great challenge in general. To the best of our
knowledge, up to date, there have been only a few works concerning NSFD schemes
preserving the GAS of continuous models. Many results on NSFD schemes only
consider the preservation of the LAS of equilibrium points (see, e.g., [21, 25, 27–29,
56, 57])
190 trang |
Chia sẻ: tueminh09 | Ngày: 22/01/2022 | Lượt xem: 584 | Lượt tải: 0
Bạn đang xem trước 20 trang tài liệu Development of nonstandard finite difference methods for some classes of differential equations, để xem tài liệu hoàn chỉnh bạn click vào nút DOWNLOAD ở trên
cos(jθλ) + 2 cos(θλ).
(3.2.6)
Obviously, P2p−1(ϕ, λ) is a polynomial of ϕ having the highest degree 2p− 1 with the
coefficients depending on rλ and θλ, and the free term is 2 cos(θλ).
Consider two cases of the stability of y∗:
Case 1. y∗ is a LAS equilibrium point of (3.2.1). Then, Re(λi) < 0 for any λi ∈ σ(J)
(i = 1, n). The necessary and sufficient condition for y∗ to be a LAS equilibrium point
of (3.2.3) is |µi| < 1 for any µi ∈ σ(Jˆ) (i = 1, n). Since Re(λ) < 0, there holds
cos(θλ) < 0. Since limϕ→0P2p−1(ϕ, λ) = 2 cos(θλ) < 0, the definition of limit of a
function follows that there exists a number ϕ1 > 0 such that P2p−1(ϕ, λ) < 0 for all
ϕ ∈ (0, ϕ1), where P2p−1(ϕ, λ) is defined by (3.2.6).
Case 2. y∗ is a linearly unstable equilibrium of (3.2.1). Then, there exists λi ∈ σ(J)
for some 1 ≤ i ≤ n such that Re(λi) > 0. The necessary and sufficient condition for
y∗ to be a linearly unstable equilibrium point of (3.2.3) is the existence of some j such
143
that |µj| > 1. Suppose for i = l there holds Re(λl) > 0. The necessary and sufficient
condition for |µl| > 1 is P2p−1(ϕ, λl) > 0, where P2p−1(ϕ, λ) is defined by (3.2.6).
Since Re(λl) > 0, we have cos(θλl) > 0. Therefore, there exists a number ϕ2 > 0
such that P2p−1(ϕ, λl) > 0 for all ϕ ∈ (0, ϕ2).
Denote Ω+ =
⋃
y∗∈Γ+ σ(J(y
∗)), Ω− =
{
ξ ∈ ⋃y∗∈Γ− σ(J(y∗)) : Re(ξ) > 0},
where Γ+ and Γ− are the set of LAS equilibria and the set of linearly unstable equilibria
of (3.2.1), respectively. Set ϕ∗ = min
{
ϕ+, ϕ−
}
, where
ϕ+ = min
{
ϕ+(λ) : λ ∈ Ω+
}
, ϕ− = min
{
ϕ−(λ) : λ ∈ Ω−
}
,
ϕ+(λ) = sup
ϕ+>0
{
ϕ+ : P2p−1(ϕ, λ) < 0, ∀ϕ ∈ (0, ϕ+), λ ∈ Ω+
}
,
ϕ−(λ) = sup
ϕ−>0
{
ϕ− : P2p−1(ϕ, λ) > 0,∀ϕ ∈ (0, ϕ−), λ ∈ Ω−
}
.
Then, from Case 1 and Case 2, we have that if 0 < ϕ < ϕ∗, then (3.2.3) is elementary
stable. The proof is complete.
The following theorem for the case p < s is stated and proved in a similar way
as Theorem 3.6.
Theorem 3.7. Suppose that the original ESRK method (A, bT , h) has the order p <
s. Then, there exists a number ϕ∗ := ϕ∗(p,Ω, A, bT , s) > 0 such that the ENRK
method(3.2.3) is elementary stable if the function ϕ(h) satisfies the condition 0 <
ϕ(h) 0.
Proof. Since p < s, the eigenvalues µi ∈ Jˆ (i = 1, n) corresponding to the eigenvalues
λi ∈ J are defined by
µi =
p∑
j=0
zj
j!
+
s∑
j>p
zjbTAj−11 = 1+z+. . .+
zp
p!
+
s∑
j>p
ajz
j, z = ϕλi, aj := aj(p,A, b
T , s).
(3.2.7)
Notice that the right hand side of (3.2.7) is the stability function of an ESRK method
having the order p < s [3, Section 4.4], [19, Section IV.2]. Repeating the proof of
Theorem 3.6 with the attention that now P2p−1(ϕ, λ) defined by (3.2.6) is replaced by
the polynomial P2s−1(ϕ, λ) of the degree 2s− 1, we obtain to the conclusion of the
theorem.
144
Remark 3.1. It is well-known that ESRK methods have to use step sizes small enough to
guarantee the LAS, i.e., ESRK methods are conditionally stable. Meanwhile, according
to Theorems 3.6 and 3.7, ENRK methods with the appropriately chosen denominator
function are unconditionally stable. Notice that if p = s, then ϕ∗ depends only on p
but if p < s, then it depends also on A, bT and s. From the proofs of the theorems, it
is possible to design an algorithm to determine the number ϕ∗ based on the finding
of the minimal positive root of the polynomials P2p−1(ϕ, λ) and P2s−1(ϕ, λ). In the
cases p = s = 1 and p = s = 2, the explicit formulas for ϕ∗ are given in [28, 29].
Therefore, Theorems 3.6 and 3.7 can be considered as a generalization of the results
constructed in [28, 29].
3.2.2. Positive ENRK methods
In this subsection, we propose a method to construct positive ENRK methods
based on results of the positivity of Runge-Kutta methods presented in Subsection
1.3.3.
Suppose that the right-hand side f(y) of (3.2.1) satisfies conditions such that
its solutions are nonnegative for all y0 ≥ 0, i.e, f ∈ P . Now suppose that the right-
hand side f belongs to one of the sets F∗, F∗(ρ), F∗∞(ρ), F∗∞,w(ρ) or Pα. Then, for
explicit Runge-Kutta methods with R(A, b) > 0, we can determine a positivity step
size thresholds depending on R(A, b). Suppose this number is H > 0. Combining this
fact with Theorem 3.6 and Theorem 3.7, we deduce that ESRK methods (3.2.3) are
PES for (3.2.1) if
ϕ(h) 0. (3.2.8)
Notice that the above conditions for existence of a positivity step size thresholdsH > 0
is only sufficient conditions. For many Runge-Kutta methods, although the number
H > 0 cannot be determined by this way, there exists a number H∗ > 0 such that the
method is positive for any h 0 determined by the
mentioned method may be not strict positivity step size thresholds (see [94]).
Next, based on results of the number R(A, b) formulated in [95, Section 9], we
choose some ESRK having R(A, b) > 0 for determining positivity step size thresholds
as follows.
145
(i) For 1-stage methods, only the Euler method has R(A, b) = 1.
(ii) For 2-stage methods, we choose the RK2 methods having R(A, b) = 1.
(iii) Among 4-stage methods, we choose the 4-stage method of order 3 (RK43) having
maximalR(A, b). This is the method defined in [95, Section 9] withR(A, b) = 2:
b1 = b2 = b3 =
1
6
, b4 =
1
2
, a21 =
1
2
, a31 = a32 =
1
2
, a41 = a42 = a43 =
1
6
.
(iv) According to [95, Theorem 9.6], there does not exist explicit 4-stage coefficient
scheme (A, b) with classical order p = 4 and R(A, b) > 0. Therefore, we
choose 5-stage method with p = 4 (RK54). This method has maximal R(A, b) ≈
1.50818 (see [95, p. 522, Section 9]).
Remark 3.2. The condition R(A, b) > 0 narrows down possible ESRK methods. For
example, two well-known 4-stage methods, namely, the classical Runge-Kutta method
and the 3/8-rule do not belong to the methods having R(A, b) > 0. However, in
numerical simulations, it will be seen that these sufficient conditions may be freed.
3.2.3. The choice of the denominator function
In this section, we analyze the influence of denominator functions to select
appropriate functions for ENRK methods.
First, notice that the standard denominator function ϕ(h) = h does not satisfy
(3.2.8). It is easy to choose a function ϕ(h) satisfying (3.2.8), for example (see [21,
25, 28, 29])
ϕ1(h) =
1− e−τ1h
τ1
, τ1 > 0. (3.2.9)
However, the change of the denominator function may cause the decrease of the order
of accuracy of ENRK methods, i.e., it does not preserve the order of accuracy of the
original ESRK. Therefore, it is important to choose the function ϕ(h) so that the order
of accuracy of the original ESRK methods is preserved.
Theorem 3.8. If the original ESRK method (A, bT , h) has the order of accuracy p,
then the local truncation error of ENRK method (A, bT , ϕ) is O(hp), whenever
ϕ(h) = h+O(hp+1), h→ 0. (3.2.10)
146
Proof. First, note that if ϕ(h) satisfies (3.2.10), then ϕ′(0) = 1 and ϕ(0) = ϕ′′(0) =
ϕ(3)(0) = . . . = ϕ(p)(0) = 0. Applying the method for constructing order conditions
of ESRK methods based on the Taylor expansion [19, Chapter II], it is easy to deduce
that the order conditions for (A, bT , h) and (A, bT , ϕ(h)) are the same if ϕ(h) satisfies
(3.2.10). Thus, the theorem is proved.
Remark 3.3. As will be seen later in Tables 3.4-3.8, the condition (3.2.10) is needed
to guarantee that the ENRK methods (A, bT , ϕ) also have the order p.
Now, to choose the function ϕ(h) satisfying simultaneously (3.2.8) and (3.2.10),
we consider the class of functions
ϕ2(h) = he
−τ2hm, m ∈ Z+, m ≥ p, τ2 > 0. (3.2.11)
Clearly, ϕ2(h) satisfies (3.2.10) and reaches the maximal value at h∗ = m
√
1/mτ2,
i.e., for any h > 0, we have ϕ2(h) ≤ ϕ(h∗) = e−1/m m
√
1/mτ2 → 0 as τ2 →∞. This
means that there always exists τ2 > 0 such that ϕ2(h) satisfies (3.2.8). Here, it suffices
to choose τ2 > τ 2opt := [me(τ
∗)m]−1. Equation (3.2.9) follows that the condition for
ϕ1(h) to satisfy (3.2.8) is τ1 > τ 1opt := (τ
∗)−1.
It is easy to see that ϕi(τi)→ h as τi → 0 (i = 1, 2). This means that when h
and τi are small, then ENRK methods (with the denominator functions ϕi(h)) have the
order of accuracy equals to that of the original ESRK method. In other words, in order
to ensure the order of accuracy, τi must be chosen as small as possible. However, the
choice of τi depends on the value of τ iopt.
It is seen that the constraint τ1 > τ 1opt := (τ
∗)−1 does not allow to choose τ1
arbitrarily small, especially when τ ∗ is very small. This is the reason leading to the
decrease of the order of accuracy of the original ESRK when choosing the function
ϕ1(h) defined by (3.2.9). The numerical simulations in the next section will clearly
demonstrate this fact.
Concerning the functionϕ2(h), we see that if τ ∗ ≥ 1, then τ 2opt = [me(τ ∗)m]−1 →
0 asm→∞. Therefore, if m is chosen large enough, then τ 2opt is very small. When
τ ∗ < 1 then τ 2opt →∞ asm→∞, and hence the choicem = p is best. If τ ∗ is very
small (e.g. for stiff problems [3, 19]) then τ 2opt is also very large. Then, ESRK methods
should use the step-size h < τ ∗ to guarantee the stability and the order of accuracy.
147
Since ϕ2(h) satisfies (3.2.10), ENRK methods also have the same order of accuracy
as the original ESRK when h < τ ∗. When h ≥ τ ∗, ESRK methods become unstable,
while, ENRK methods is still stable. Thus, if h is small, then ϕ2(h) gives the accuracy
better than ϕ1(h).
However, since ϕ2(h)→ 0 as h→∞, whenever h is large, ϕ2(h) ≈ 0. Then,
numerical solutions obtained at all steps are slightly different from the initial value.
Therefore, it is impossible to estimate the asymptotic behavior of the solutions. In
other words, it is not recommended to use the function ϕ2(h) for large h. In this case,
the function ϕ1(h) is advantageous because it is a monotonically increasing and upper
bounded function. As h→∞, there holds ϕ1(h)→ 1/τ1 < τ ∗. Therefore, although h
is large enough, the numerical solutions remain to have the LAS as the exact solutions
(see [28, 29]).
In order to overcome the shortcoming of ϕ2(h), we propose a class of new
functions ϕ(h) to guarantee the order when h is small as well as the asymptotic
behavior of numerical solutions when h is large. It is the class of functions of the form
ϕ3(h) = θ(h)ϕ2(h) +
(
1− θ(h))ϕ1(h), (3.2.12)
where the function θ(h) has the property θ(h) = 1 +O(hp) as h→ 0, 0 < θ(h) < 1
for all h > 0, limh→0 θ(h) = 1 and limh→∞ θ(h) = 0.
Obviously, ϕ3(h) satisfies (3.2.10), and if ϕ1(h) and ϕ2(h) satisfy (3.2.8), then
so does ϕ3(h). Especially, when h is small, then ϕ3(h) is equivalent to ϕ2(h), so this
guarantees the best error, conversely, when h is large, then ϕ3(h) is equivalent to
ϕ1(h) and this guarantees the asymptotic behavior of numerical solutions. It is best
to choose θ(h) = e−h
k
, where k ∈ Z+. Figure 3.3 depicts the graphs of the function
ϕi(h)(i = 1, 2, 3 for some particular values of the parameters.
148
h
0 0.5 1 1.5 2 2.5 3
0
0.2
0.4
0.6
0.8
1
ϕ
∗
ϕ1
ϕ2
ϕ3
h
0 0.5 1 1.5 2 2.5 3
0
0.2
0.4
0.6
0.8
1
ϕ
∗
ϕ1
ϕ2
ϕ3
Figure 3.3. Graphs of the functions ϕi(h) in two cases of the paramerters. In the
upper figure: ϕ∗ = 1, ϕ1 = 1 − e−h, ϕ2 = he−0.12h4 , ϕ3 = (1 − e−h3)ϕ1 + e−h3ϕ2.
In the lower figure: ϕ∗ = 1/1.2, ϕ1 = (1 − e−1.2h)/1.2, ϕ2 = he−0.2h5 , ϕ3 =
(1− e−h4)ϕ1 + e−h4ϕ2.
3.3. Some applications of the ENRK methods
In this section, we apply the constructed ENRK methods to some important
mathematical models. Additionally, numerical experiments are performed to confirm
the validity of the theoretical results.
3.3.1. ENRK methods for a predator-prey system
We consider the following predator-prey system with Beddington-DeAngelis
functional response [107], which was considered in [29]
dx
dt
= x− Axy
1 + x+ y
,
dy
dt
=
Exy
1 + x+ y
−Dy, (x(0), y(0)) = (1, 1.6), (3.3.1)
where x and y represent the prey and predator population sizes, respectively, and the
values of the constants are A = 2, D = 1 and E = 10. A mathematical analysis of
149
System (3.3.1) shows that there exist two equilibria
E0 := (0, 0), E∗ :=
( AD
AE − E − AD,
E
AE − E − AD
)
= (0.25, 1.25),
where the equilibrium (0.25, 1.25) is GAS in the interior of the first quadrant, and the
equilibrium (0, 0) is unstable. The eigenvalues of J(0, 0) are λ1 = 1 and λ2 = −1,
and the eigenvalues of J(0.25, 1.25) are
λ3,4 = −1
5
± 3
5
i = r
(
cos(θ)± i sin(θ)
)
, r =
√
2
5
, θ = 0.6024pi.
Since it is impossible to find the exact solution of the system, we use the numerical
solution obtained by a 11-stage Runge-Kutta method of the order 8 (RK8) [108]
with h = 10−5 as a benchmark solution. Here, error = maxk
{|xk − Xk| + |yk −
Yk|
}
is used as a measure for the accuracy of ENRK methods, where {(xk, yk)} and
{(Xk, Yk)} are the solutions obtained by ENRK methods and the benchmark solution,
respectively. Besides, rate := logh1/h2(error(h1)/error(h2)) (see [3, Example 4.1])
is an approximation for the order of accuracy of the methods.
In this example, we will consider ENRK1 (based on the Euler method), ENRK2
(based on the second Heun method), ENRK43 (based on the RK43 method), ENRK54
(based on the RK54 method). The numbers R(A, b) for these methods are 1, 1, 2 and
1.50818, respectively. Besides, ENRK4 (based on the classical 4-stage RK method) is
also considered although this method does not possess R(A, b) > 0.
It is easy to verify that the right-hand side f of the model belongs to the set
Pα := {f |f(t, v) + αv ≥ 0 for all t, v ≥ 0} with α = max{A − 1, D} = 1, and
therefore, it is easy to determine the positivity step size thresholds for ENRK methods.
Moreover, it is not difficult to determine the polynomialsP2p−1 andP2s−1(ϕ, λi)
(i = 1, 2, 3, 4) corresponding to ENRK methods. Therefore, it is easy to determine
the number ϕ∗. The numbers ϕ∗, τ iopt (i = 1, 2) and the functions ϕj(h) (j = 1, 2, 3)
for the methods ENRK1, ENRK2, ENRK43, ENRK54, ENRK4 are given in Table
3.3. The errors and the rates of the methods for small h are reported in Tables 3.3-3.8,
where for short, the columns with the headings ϕ(h) and ϕi stand for the errors of the
methods with the denominator functions ϕ(h) and ϕi, respectively; the columns with
the headings ratei stand for the rates of the methods corresponding to ϕi.
150
Notice that it is impossible to determine the positivity threshold for ENRK4
method because R(A, b) = 0. In this case, it is only possible to determine its el-
ementary stability threshold. However, many numerical simulations show that the
elementary stability threshold is also the positivity threshold. This fact is also seen
from the numerical simulations in [29].
Table 3.3. The values τ iopt and the denominator functions ϕi(h) (i = 1, 2, 3) of the
ENRK methods
Method (s, p) ϕ∗ H τ∗ τ1opt ϕ1(h) τ2opt ϕ2(h) θ(h) in (3.2.12)
ENRK1 (1, 1) 0.9998 1 0.9998 1.0002
1− e−1.0005h
1.0005
0.0919 he−0.095h
4
e−0.01h
2
ENRK2 (2, 2) 2.6604 1 1 1 1− e−h 0.0920 he−0.095h4 e−0.01h4
ENRK43 (4, 3) 4.7332 2 2 0.5
1− e−0.55
0.55
9.5802e-004 he−0.001h
6
e−h
6
ENRK54 (5, 4) 5.0631 1.50818 1.50818 0.6631
1− e−0.68
0.68
1.7179e-003 he−0.002h
8
e−h
8
ENRK4 (4, 4) 4.4476 * 4.4476 0.2248
1− e−0.25
0.25
7.9214e-006 he−0.0001h
6
e−0.01h6
Table 3.4. The errors and rates of ENRK1 methods
h ϕ(h) = h ϕ1 rate1 ϕ2 rate2 ϕ3 rate3
0.2 0.4303 0.6056 0.4304 0.4304
0.1 0.2032 0.2937 1.0443 0.2032 1.0827 0.2032 1.0827
0.05 0.0986 0.1444 1.0239 0.0986 1.0439 0.0986 1.0439
0.01 0.0192 0.0285 1.0091 0.0192 1.0156 0.0192 1.0156
0.005 0.0096 0.0142 1.0027 0.0096 1.0045 0.0096 1.0045
0.001 0.0019 0.0028 1.0009 0.0019 1.0016 0.0019 1.0016
151
Table 3.5. The errors and rates of ENRK2 methods
h ϕ(h) = h ϕ1 rate1 ϕ2 rate2 ϕ3 rate3
0.2 7.3223e-003 4.1755e-001 7.1013e-003 7.0992e-003
0.1 1.7189e-003 2.1136e-001 0.9823 1.7052e-003 2.0581 1.7051e-003 2.0578
0.05 4.1773e-004 1.0622e-001 0.9926 4.1687e-004 2.0323 4.1686e-004 2.0322
0.01 1.6354e-005 2.1321e-002 0.9978 1.6352e-005 2.0121 1.6352e-005 2.0121
0.005 4.0770e-006 1.0665e-002 0.9994 4.0769e-006 2.0040 4.0769e-006 2.0040
0.001 1.6271e-007 2.1337e-003 0.9998 1.6271e-007 2.0014 1.6271e-007 2.0014
Table 3.6. The errors and rates of ENRK43 methods
h ϕ(h) = h ϕ1 rate1 ϕ2 rate2 ϕ3 rate3
0.2 5.8286e-004 1.9063e-001 5.8275e-004 5.7796e-004
0.1 7.1911e-005 9.5672e-002 0.9946 7.1910e-005 3.0186 7.1872e-005 3.0075
0.05 8.9428e-006 4.7924e-002 0.9973 8.9428e-006 3.0074 8.9425e-006 3.0067
0.01 7.1300e-008 9.5989e-003 0.9991 7.1300e-008 3.0021 7.1300e-008 3.0021
0.005 8.9081e-009 4.8003e-003 0.9997 8.9081e-009 3.0007 8.9081e-009 3.0007
0.001 7.1181e-011 9.6021e-004 0.9999 7.1181e-011 3.0007 7.1181e-011 3.0007
Table 3.7. The errors and rates of ENRK54 methods
h ϕ(h) = h ϕ1 rate1 ϕ2 rate2 ϕ3 rate3
0.2 3.1359e-005 2.8632e-001 3.1368e-005 3.1665e-005
0.1 2.0695e-006 1.4419e-001 0.9897 2.0695e-006 3.9219 2.0700e-006 3.9352
0.05 1.3274e-007 7.2338e-002 0.9952 1.3274e-007 3.9626 1.3274e-007 3.9629
0.01 2.1686e-010 1.4502e-002 0.9985 2.1686e-010 3.9871 2.1686e-010 3.9871
0.005 1.3706e-011 7.2531e-003 0.9996 1.3706e-011 3.9839 1.3706e-011 3.9839
0.001 1.9159e-012 1.4510e-003 0.9999 1.9159e-012 1.2225 1.9159e-012 1.2225
Table 3.8. The errors and rates of ENRK4 methods
h ϕ(h) = h ϕ1 rate1 ϕ2 rate2 ϕ3 rate3
0.2 1.9481e-005 1.0622e-001 1.9488e-005 2.1385e-005
0.1 1.1945e-006 5.3233e-002 0.9967 1.1946e-006 4.0280 1.2044e-006 4.1502
0.05 7.3021e-008 2.6646e-002 0.9984 7.3022e-008 4.0321 7.3099e-008 4.0423
0.01 1.1429e-010 5.3336e-003 0.9995 1.1429e-010 4.0137 1.1430e-010 4.0143
0.005 7.2312e-012 2.6671e-003 0.9998 7.2312e-012 3.9824 7.2312e-012 3.9824
0.001 1.9159e-012 5.3346e-004 0.9999 1.9159e-012 0.8253 1.9159e-012 0.8253
152
From Tables 3.4-3.8, we see that the replacement of the denominator function
ϕ(h) = h by the function ϕ1(h) decreases the order of accuracy of the original
ESRK methods. Namely, ENRK methods have only the order 1. Meanwhile, ENRK
methods with the appropriate denominator functions ϕi(h) (i = 2, 3) have the orders
of accuracy equal to those of the original ESRK methods. In other words, the orders of
accuracy of the original ESRK methods are preserved.
In the columns rate2 and rate3 of Tables 3.7 and 3.8, we see an unexpected
phenomenon, namely, the rates decrease when h is small. A similar phenomenon
also was indicated in [3, Example 4.1] when studying explicit standard Runge-Kutta
methods. The reason of this is that the rounding errors generally increase as h decreases.
Next, we consider ENRK54 methods for large h, specifically h = 4. Then
ϕ2(4) ≈ 4.7667e−057, therefore it is possible to accept ϕ2(4) = 0. The computational
experiments show that the numerical solutions obtained when using the function ϕ2(h)
at all grid nodes are constant and equal to the initial values. Meanwhile, the numerical
solutions obtained by the methods with the use of the functions ϕ1(h) and ϕ3(h) have
the asymptotic behavior similar to that of the exact solution. Figure 3.4 depicts the
numerical solution obtained by ENRK54 method with
ϕ3(h) = e
−h8he−0.002h
6
+ (1− e−h8)1− e
−0.68h
0.68
, h = 4.
The experiments for the other methods give similar results. This completely agrees
with the analysis made in Subsection 3.2.3. Notice that from Figure 3.4, we see that
the GAS of the model is also preserved. In the figure, each curve (blue, red, yellow,
. . . ) represents a trajectory (or orbit) corresponding to a specific initial data.
153
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8
0
1
2
3
4
5
6
7
8
x
y
Figure 3.4. Phase planes for the model (3.3.1) with some different inital data obtained
by ENRK54 method with ϕ3(h) and h = 4.
Finally, for the purpose of comparison, we apply the method proposed by Wood
and Kojouharov [27] to the model. This is an NSFD scheme preserving the positivity
and elementary stability of the dynamical system based on non-local approximation.
For this method, the positivity step size threshold is H = ∞, and its elementary
stability threshold is the same as Euler method. Following the results in [27], we choose
the denominator function ϕ(h) = (1− e−1.0005h)/1.0005. The error and convergence
rate of the method are reported in Table 3.9. From the table, we see that the method
has the order 1 but the Euler method is somewhat better (see Table 3.4).
154
Table 3.9. The errors and rates of the Wood and Kojouharov methods.
h ϕ(h) = h ϕ1 rate1 ϕ2 rate2 ϕ3 rate3
0.2 0.5390 0.8143 0.5391 0.5621
0.1 0.2575 0.4041 1.0107 0.2576 1.0658 0.2918 0.9460
0.05 0.1257 0.2008 1.0093 0.1257 1.0348 0.1538 0.9237
0.01 0.0247 0.0399 1.0038 0.0247 1.0121 0.0326 0.9644
0.005 0.0123 0.0199 1.0012 0.0123 1.0036 0.0164 0.9889
0.001 0.0025 0.0040 1.0004 0.0025 1.0012 0.0033 0.9962
3.3.2. ENRK methods for a vaccination model with multiple endemic states
We consider the following vaccination model with multiple endemic states
[109], which was considered in [29]
dS
dt
= µN − βSI/N − (µ+ φ)S + cI + δV,
dI
dt
= βSI/N − (µ+ c)I,
dV
dt
= φS − (µ+ δ)V,
where the constants β = 0.7, c = 0.1, µ = 0.8, δ = 0.8 and φ = 0.8. In the above
model, the total (constant) population size N = 100 is divided into three classes
susceptibles (S), infectives (I) and vaccinated (V ), and it is assumed that the vaccine
is completely effective in preventing infection. A mathematical analysis of this system
shows that the disease free equilibrium (S∗, I∗, V ∗) =
( (µ+ δ)N
µ+ δ + φ
, 0,
φN
µ+ δ + φ
)
=(200
3
, 0,
100
3
)
is GAS [109]. The eigenvalues of J(S∗, I∗, V ∗) are given by λ1 =
−0.8, λ2 = −2.4, λ3 = −13/30. Therefore, it is easy to determine the elementary
stability threshold for ENRK methods. Moreover, the right-hand side of the system
belongs to the set Pα where α = max{β + µ + φ, µ + c, µ + δ} = 2.5, and hence,
we can determine the positivity thresholds for ENRK methods. Table 3.10 gives the
positivity and elementary stability thresholds for ENRK methods.
155
Table 3.10. Positivity and elementary stability thresholds for ENRK
Methods ENRK1 ENRK2 ENRK43 ENRK54 ENRK4
ϕ∗ 0.8333 0.8333 2.1499 2.2068 1.1605
H = R(A, b)/α 0.4 0.4 0.8 0.6631 *
τ ∗ = min{ϕ∗, H} 0.4 0.4 0.8 0.6631 *
Based on these results, we choose the denominator function for ENRK methods.
The numerical solutions obtained by ENRK54 are depicted in Figure 3.5, where each
color curve represents a trajectory corresponding a specific initial data. We see that
the GAS of the model is preserved, while the numerical simulations in [29] show that
the RK2 and Euler methods do not preserve this property. In addition, the advantage
of the ENRK methods is that they preserve the essential properties of the model for all
h > 0 and have high order of accuarcy when h is small.
30
40
50
60
70
0
10
20
30
40
50
0
10
20
30
40
50
SI
V
Figure 3.5. Phase portrait for the vaccination model with some different initial data
obtained by ENRK54 method for ϕ3(h) = e−h
6
he−0.5h
4
+ (1− e−h6)(1− e−1.6h)/1.6
and h = 2.
156
3.4. Conclusions
In this chapter, we have constructed EFD schemes and high order NSFD
schemes for a class of general dynamical systems based on the general standard
Runge-Kutta methods. The obtained results can be considered as a generalization of
the results formulated in [28, 29, 42]
Firstly, implicit and explicit EDS schemes for systems of three linear ODEs with
constant coefficients are constructed. Importantly, the obtained results not only answer
the open question posted by Roeger [40] but also can be extended to design EFD
schemes for general n-dimensional systems of linear ODEs with constant coefficients.
Secondly, we have constructed and analyzed high order ENRKmethods preserv-
ing two important properties of general autonomous dynamical systems, namely, the
positivity and LAS. The main result resolved the contradiction between the dynamics
consistency and high order of accuracy of NSFD schemes. Additionally, two important
applications of the constructed ENRK methods to the predator-prey model and the
vaccination model are also presented.
Lastly, the validity of the theoretical results and the superiority of the proposed
NSFD schemes are supported by many numerical examples. The results also indicate
that there is a good agreement between the numerical examples and the theoretical
results.
157
GENERAL CONCLUSIONS
In this thesis, we have successfully developed the Mickens’ methodology to con-
struct nonstandard finite difference (NSFD) methods for solving some important
classes of differential equations arising in fields of science and technology. The pro-
posed NSFD schemes are not only dynamically consistent with the differential equation
models, but also easy to be implemented; furthermore, they can be used to solve a
large class of mathematical problems in both theory and practice. The validity of the
theoretical results and the superiority of the NSFD schemes have been confirmed and
supported by many numerical simulations. The results have indicated that there is a
good agreement between the theoretical aspect and experimental one.
In the first part, we have successfully constructed NSFD schemes for some
mathematical models described by systems of ODEs including two metapopulation
models, one predator-prey model and two computer virus propagation models. It is
worth noting that all of the models possess at least one of the following characteristics:
(i) having large dimensions.
(ii) having non-hyberbolic equilibrium points.
(iii) having the GAS.
Firstly, we have investigated the GAS of the constructed NSFD schemes for the
metapopulation model formulated in [97] by using the standard techniques of math-
ematical analysis. Secondly, we have used the Lyapunov stability theorem to study
the GAS of proposed NSFD schemes for the computer virus propagation model and
the general predator-prey model constructed in [12] and [105], respectively. Lastly,
we have proposed two novel approaches to establish the stability properties of the
NSFD schemes for the metapopulation model and the propagation model of computer
viruses formulated in [98] and [16], respectively. The first approach is based on the
extension of the classical Lyapunov stability theorem, and the second one is based
on the Lyapunov stability theorem and its extensions in combination with a theorem
on the GAS of discrete-time nonlinear cascade systems. Both approaches lead the
study of the stability of the proposed NSFD schemes to the study of the stability of
158
discrete models with smaller dimension, and therefore, complicated calculations and
transforms are limited significantly.
In the second part, we have constructed EFD schemes and high order NSFD
schemes for a class of general dynamical systems based on the general Runge-Kutta
methods. The obtained results can be considered as a generalization of the results
formulated in [28, 29, 42]. Firstly, implicit and explicit EDS schemes for systems
of three linear ODEs with constant coefficients are constructed. Importantly, the
obtained results not only answer the open question posted by Roeger [40] but also
can be extended to design EFD schemes for general n-dimensional systems of linear
ODEs with constant coefficients. Next, we have constructed and analyzed high order
ENRK methods preserving two important properties of general autonomous dynamical
systems, namely, the positivity and LAS. The main result resolved the contradiction
between the dynamics consistency and high order of accuracy of NSFD schemes.
Additionally, two important applications of the constructed ENRK methods to the
predator-prey model and the vaccination model are also presented.
In the near future, the established results in this thesis will be developed to
construct highly effective NSFD schemes for PDEs, DDEs, FDEs and stochastic
differential equations. Also, we intent to study the combination of the Mickens’
methodology and other existing approaches to create new numerical methods with
high performance for both differential equations and integro-differential equations.
159
THE LIST OF THEWORKS OF THE AUTHOR
RELATED TO THE THESIS
[A1] Quang A Dang,Manh Tuan Hoang, Dynamically consistent discrete
metapopulation model, Journal of Difference Equations and Applications, 2016, 22,
1325-1349, (SCI-E).
[A2] Quang A Dang,Manh Tuan Hoang, Lyapunov direct method for investigating
stability of nonstandard finite difference schemes for metapopulation models, Journal
of Difference Equations and Applications, 2018, 24, 15-47, (SCI-E).
[A3] Quang A Dang,Manh Tuan Hoang, Complete global stability of a
metapopulation model and its dynamically consistent discrete models, Qualitative
Theory of Dynamical Systems, 2019, 18, 461-475, (SCI-E)
[A4] Quang A Dang,Manh Tuan Hoang, Numerical dynamics of nonstandard finite
difference schemes for a computer virus propagation model, International Journal of
Dynamics and Control, 2020, 8, 772-778, (SCOPUS).
[A5] Quang A Dang,Manh Tuan Hoang, Nonstandard finite difference schemes for
a general predator-prey system, Journal of Computational Science, 2019, 36, 101015,
(SCI-E).
[A6] Quang A Dang,Manh Tuan Hoang, Positivity and global stability preserving
NSFD schemes for a mixing propagation model of computer viruses, Journal of
Computational and Applied Mathematics 2020, 374, 112753, (SCI).
[A7]Manh Tuan Hoang, On the global asymptotic stability of a predator-prey model
with Crowley-Martin function and stage structure for prey, Journal of Applied Mathe-
matics and Computing, 2020, 64, 765-780, (SCI-E).
[A8] Quang A Dang,Manh Tuan Hoang, Exact finite difference schemes for
three-dimensional linear systems with constant coefficients, Vietnam Journal of
Mathematics, 2018, 46, 471-492, (ESCI, SCOPUS).
160
[A9] Quang A Dang,Manh Tuan Hoang, Positive and elementary stable explicit
nonstandard Runge-Kutta methods for a class of autonomous dynamical systems,
International Journal of Computer Mathematics, 2020, 97, 2036-2054, (SCI-E).
161
Bibliography
1. R. P. Agarwal, An Introduction to Ordinary Differential Equations, Springer, 2000.
2. L.J.S. Allen, An Introduction to Mathematical Biology, Prentice Hall, 2007, New
Jersey.
3. U. M. Ascher, L.R. Petzold, Computer Methods for Ordinary Differential Equa-
tions and Differential-Algebraic Equations, Society for Industrial and Applied
Mathematics, 1998, Philadelphia.
4. F. Brauer, C. Castillo-Chavez,Mathematical Models in Population Biology and
Epidemiology, Springer, 2001, New York.
5. L. Edelstein-Keshet,Mathematical Models in Biology, Society for Industrial and
Applied Mathematics, 1998, Philadelphia.
6. H. K. Khalil, Nonlinear systems, Prentice Hall, 2002.
7. M. Martcheva, An Introduction to Mathematical Epidemiology, Springer Sci-
ence+Business Media New York, 2015.
8. R.Mattheij, J. Molenaar, Classics in Applied Mathematics: Ordinary Differential
Equations in Theory and Practice, Society for Industrial and Applied Mathematics,
2002, New York.
9. L. Perko, Differential Equations and Dynamical Systems, Third Edition, Texts in
Applied Mathematics 7, Springer, 2001.
10. H. L. Smith, P. Waltman, The Theory of the Chemostat, Cambridge University
Press, 1995, Cambridge.
11. W. H. Murray, The application of epidemiology to computer viruses, Computers
& Security, 1988, 7, 130-50.
12. L. X. Yang, X. Yang, Q. Zhu, L. Wen, A computer virus model with graded cure
rates, Nonlinear Analysis: Real World Applications, 2013, 14, 414 -422.
162
13. L-X. Yang, X. Yang, A new epidemic model of computer viruses, Communications
in Nonlinear Science and Numerical Simulation, 2014, 19, 1935-1944.
14. X. Yang, B. K. Mishra, Y. Liu, Computer virus: theory, model, and methods,
Discrete Dynamics in Nature and Society, 2012, Article ID 473508.
15. L. X. Yang, X. Yang, L. Wen, J. Liu, A novel computer virus propagation model
and its dynamics, International Journal of Computer Mathematics, 2012, 89,
2307-2314.
16. Q. Zhu, X. Yang, L. Yang, X. Zhang, A mixing propagation model of computer
viruses and countermeasures, Nonlinear Dynamics, 2013, 73, 1433-1441.
17. A. M. Stuart, A. R. Humphries, Dynamical Systems and Numerical Analysis,
Cambridge monographs on applied and computational mathematics 2, Cambridge
University Pres, 1998, New York.
18. R. L. Burden, J. D. Faires, Numerical analysis, Brooks/Cole, Cengage Learning,
2011.
19. E. Hairer, P. S. Norsett, G. Wanner, Solving Ordinary Differential Equations I:
Nonstiff Problems, Springer Series in Computational Mathematics 8, Springer-
Verlag Berlin Heidelberg, 1993.
20. E. Hairer, G. Wanner, Solving Ordinary Differential Equations II: Stiff and
Differential-Algebraic Problems, Springer Series in Computational Mathematics
14, Springer-Verlag Berlin Heidelberg, 1996.
21. R. E. Mickens, Nonstandard Finite Difference Models of Differential Equations,
World Scientific, 1994.
22. R. E. Mickens, Applications of Nonstandard Finite Difference Schemes, World
Scientific, 2000.
23. R. E. Mickens, Advances in the Applications of Nonstandard Finite Difference
Schemes, World Scientific, 2005.
163
24. R. E. Mickens, Nonstandard Finite Difference Schemes for Differential Equations,
Journal of Difference Equations and Applications, 2012, 8, 823-847.
25. R. Anguelov, J. M. S. Lubuma, Nonstandard finite difference method by nonlocal
approximations, Mathematics and Computers in Simulation, 2003, 61, 465-475.
26. J. Cresson, F. Pierret, Non standard finite difference scheme preserving dynamical
properties, Journal of Computational and Applied Mathematics, 2016, 303, 15-30.
27. D. T. Wood, H. V. Kojouharov, A class of nonstandard numerical methods for
autonomous dynamical systems, Applied Mathematics Letters, 2015, 50, 78-82.
28. D. T. Dimitrov, H. V. Kojouharov, Nonstandard finite-difference schemes for
general two-dimensional autonomous dynamical systems, Applied Mathematics
Letters, 2005, 18, 769-774.
29. D. T. Dimitrov, H. V. Kojouharov, Stability-preserving finite-difference methods for
general multi-dimensional autonomous dynamical systems, International journal
of numerical analysis and modeling, 2007, 4, 280-290.
30. P. Liu and S. N. Elaydi, Discrete competitive and cooperative models of Lotka-
Volterra type, Journal of Computational Analysis and Applications, 2001, 3, 53-73.
31. R. E. Mickens, K. Oyedeji, S. Rucker, Exact finite difference scheme for second-
order, linear ODEs having constant coefficients, Journal of Sound and Vibration,
2005, 287, 1052-1056.
32. R. E. Mickens, Exact finite difference schemes for two-dimensional advection
equations, Journal of Sound and Vibration, 1997, 207, 426-428.
33. L.-I. W. Roeger, R. E. Mickens, Exact finite-difference schemes for first order
differential equations having three distinct fixed-points, Journal of Difference
Equations and Applications, 2007, 13, 1179-1185.
34. L. -I. W. Roeger, R. E. Mickens, Exact finite difference and non-standard finite
difference schemes for dx/dt = −λyα, Journal of Difference Equations and
Applications, 2012, 18, 1511-1517.
164
35. L. -I. W. Roeger, Exact nonstandard finite-difference methods for a linear system-
the case of centers, Journal of Difference Equations and Applications, 2008, 14,
381-389.
36. L. -I. W. Roeger, R. E. Mickens, Exact finite difference scheme for linear differ-
ential equation with constant coefficients, Journal of Difference Equations and
Applications, 2013, 19, 1663-1670.
37. L. -I. W. Roeger, Exact finite-difference schemes for two-dimensional linear sys-
tems with constant coefficients, Journal of Computational and Applied Mathemat-
ics, 2008, 219, 102-109.
38. K. C. Patidar, On the use of nonstandard finite difference methods, Journal of
Difference Equations and Applications, 2005, 11,11 735-758.
39. K. C. Patidar, Nonstandard finite difference methods: recent trends and further
developments, Journal of Difference Equations and Applications, 2016, 22, 817-
849.
40. L.-I. W. Roeger, General nonstandard finite-difference schemes for differential
equations with three fixed-points, Computers and Mathematics with Applications,
2009, 57, 379-383.
41. L. -I. W. Roeger, Nonstandard finite difference schemes for differential equations
with n+ 1 distinct fixed-points, Journal of Difference Equations and Applications,
2009, 15, 133-151.
42. L. -I. W. Roeger, R. E. Mickens, Non-standard finite difference schemes for the
differential equation y′(t) = bnyn+bn−1yn−1+. . .+b1y+b0, Journal of Difference
Equations and Applications, 2002, 18, 305-312.
43. B. M. Chen-Charpentier, D. T. Dimitrov, H. V. Kojouharov, Combined nonstandard
numerical methods for ODEs with polynomial right-hand sides, Mathematics and
Computers in Simulation, 2006, 73, 105-113.
165
44. F. J. Solis, B. Chen-Charpentier, Nonstandard discrete approximations preserv-
ing stability properties of continuous mathematical models, Mathematical and
Computer Modelling, 2004, 40, 481-490.
45. R. E. Mickens, Discretizations of nonlinear differential equations using explicit
nonstandard methods, Journal of Computational and Applied Mathematics, 1999,
110, 1810-185.
46. R. E. Mickens, I. Ramadhani, Finite-difference schemes having the correct linear
stability properties for all finite step-sizes III, Computers & Mathematics with
Applications, 1994, 27, 77-84.
47. U. Erdogan, T. Ozis, A smart nonstandard finite difference scheme for second
order nonlinear boundary value problems, Journal of Computational Physics,
2011, 230, 6464-6474.
48. H. Kojouharov, B. Welfert, A nonstandard Euler scheme for y′′+g(y)y′+f(y) = 0,
Journal of Computational and Applied Mathematics, 2003, 151, 335-353.
49. P. A. Zegeling, S. Iqbal, Nonstandard finite differences for a truncated Bratu-
Picard model, Applied Mathematics and Computation, 2018, 324, 266-284.
50. R. E. Mickens, A nonstandard finite-difference scheme for the Lotka-Volterra
system, Applied Numerical Mathematics, 2003, 45, 309-314.
51. L. -I. W. Roeger, A nonstandard discretization method for Lotka-Volterra models
that preserves periodic solutions, Journal of Difference Equations and Applica-
tions, 2005, 11, 721-733.
52. L. -I. W. Roeger, Periodic solutions preserved by nonstandard finite-difference
schemes for the Lotka–Volterra system: a different approach, Journal of Difference
Equations and Applications, 2008, 14, 481-493.
53. L. -I. W. Roeger, Nonstandard finite-difference schemes for the Lotka-Volterra
systems: generalization of Mickens’s method, Journal of Difference Equations and
Applications, 2006, 12, 937-948.
166
54. L. -I. W. Roeger, Dynamically consistent discrete Lotka-Volterra competition
models derived from nonstandard finite-difference schemes, Discrete & Continuous
Dynamical Systems B, 2008, 9, 415-429
55. L. -I. W. Roeger, G. L. Jr, Dynamically consistent discrete Lotka-Volterra com-
petition systems, Journal of Difference Equations and Applications, 2003, 19,
191-200.
56. D. T. Dimitrov, H. V. Kojouharov, Positive and elementary stable nonstandard
numerical methods with applications to predator-prey models, Journal of Compu-
tational and Applied Mathematics, 2006, 189, 98-108.
57. D. T. Dimitrov, H. V. Kojouharov, Nonstandard finite-difference methods for
predator-prey models with general functional response, Mathematics and Com-
puters in Simulation, 2008, 78, 1-11.
58. R. Anguelov, Y. Dumont, J. M. S. Lubuma, M. Shillor, Dynamically consistent
nonstandard finite difference schemes for epidemiological models, Journal of
Computational and Applied Mathematics, 2014, 255, 161-182.
59. A. J. Arenas, G. González-Parra, B. M. Chen-Charpentier, A nonstandard nu-
merical scheme of predictor-corrector type for epidemic models, Computers and
Mathematics with Applications, 2010, 59, 3740-3749.
60. O. F. Egbelowo, Nonstandard finite difference approach for solving 3-compartment
pharmacokinetic models, International Journal for Numerical Methods in Biomed-
ical Engineering, 2018, https://doi.org/10.1002/cnm.3114.
61. S.M.Garba, A.B.Gumel, J. M. S. Lubuma, Dynamically-consistent non-standard
finite difference method for an epidemic model, Mathematical and Computer
Modelling, 2011, 53, 131-150.
62. K. F. Gurski, A simple construction of nonstandard finite-difference schemes for
small nonlinear systems applied to SIR models, Computers and Mathematics with
Applications, 2013, 66, 2165-2177.
167
63. L. Jodar, R. J.Villanueva, A. J. Arenas, G. C.González, Nonstandard numeri-
cal methods for a mathematical model for influenza disease, Mathematics and
Computers in Simulation, 2008, 79, 622-633.
64. A. Korpusik, A nonstandard finite difference scheme for a basic model of cellular
immune response to viral infection, Communications in Nonlinear Science and
Numerical Simulation, 2017, 43, 369-384.
65. D. T. Wood, H. V. Kojouharov, D. T. Dimitrov, Universal approaches to approxi-
mate biological systems with nonstandard finite difference methods, Mathematics
and Computers in Simulation, 2017, 133, 337-350.
66. R. E. Mickens, A. B. Gumel, Numerical study of a non-standard finite-difference
scheme for the Van der Pol equation, Journal of Sound and Vabration, 2000, 250,
955-963.
67. R. E. Mickens, Step-Size Dependence Of The Period For A Forward-Euler Scheme
Of The Van Der Pol Equation, Journal of Sound and Vabration, 2002, 258, 199-
202.
68. R. E. Mickens, A numerical integration technique for conservative oscillators
combining nonstandard finite-difference methods with a Hamilton’ s principle,
Journal of Sound and Vabration, 2005, 285, 477-482.
69. D. T. Wood, Advancements and applications of nonstandard finite difference
methods, The University of Texas at Arlington, https://rc.library.uta.edu/uta-
ir/handle/10106/25406.
70. O. F. Egbelowo, The Nonstandard Finite Difference Method Applied to Phar-
macokinetic Models, University of the Witwatersrand, Johannesburg, 2018,
https://hdl.handle.net/10539/27306.
71. G. Gonzalez-Parra, A. J. Arenas, B. M. Chen-Charpentier, Combination of non-
standard schemes and Richardson’s extrapolation to improve the numerical so-
lution of population models, Mathematical and Computer Modelling, 2010, 52,
1030-1036.
168
72. J. Martin-Vaquero, A. Martin del Rey, A. H. Encinas, J. D. Hernandez Guillen, A.
Queiruga-Dios, G. Rodriguez Sanchez, Higher-order nonstandard finite difference
schemes for a MSEIR model for a malware propagation, Journal of Computational
and Applied Mathematics, 2017, 317, 146-156.
73. J. Martin-Vaquero, A. Queiruga-Dios, A. Martin del Rey, A. H. Encinas, J. D.
Hernandez Guillen, G. Rodriguez Sanchez, Variable step length algorithms with
high-order extrapolated non-standard finite difference schemes for a SEIR model,
Journal of Computational and Applied Mathematics, 2018, 330, 848-854.
74. R. E. Mickens, A nonstandard finite difference scheme for the diffusionless Burgers
equation with logistic reaction, Mathematics and Computers in Simulation, 2003,
62, 117-124.
75. R. E. Mickens, A nonstandard finite difference scheme for a Fisher PDE having
nonlinear diffusion, Computers and Mathematics with Applications, 2003, 45,
429-436.
76. R. E. Mickens, A nonstandard finite difference scheme for a PDE modeling
combustion with nonlinear advection and diffusion, Mathematics and Computers
in Simulation, 2005, 69, 439-446.
77. R. E. Mickens, A nonstandard finite difference scheme for a nonlinear PDE having
diffusive shock wave solutions, Mathematics and Computers in Simulation, 2001,
55, 549-555.
78. M. Chapwanya, J. M.-S. Lubuma, R. E. Mickens, Positivity-preserving non-
standard finite difference schemes for cross-diffusion equations in biosciences,
Computers and Mathematics with Applications, 2014, 68, 1071-1082.
79. A. J. Arenas, G. González-Parra, B. Melendez Caraballo, A nonstandard finite
difference scheme for a nonlinear Black-Scholes equation, Mathematical and
Computer Modeling, 2013, 57, 1663-1670.
80. M. Ehrhardt, R. E.Mickens, A nonstandard finite difference scheme for convection-
diffusion equations having constant coefficients, Applied Mathematics and Com-
putation, 2013, 219, 6591-6604.
169
81. Y. Yang, J. Zhou, X. Ma, T. Zhang, Nonstandard finite difference scheme for a
diffusive within-host virus dynamics model with both virus-to-cell and cell-to-
cell transmissions, Computers and Mathematics with Applications, 2016, 72,
1013-1020.
82. A.J. Arenas, G. Gonzalez-Parra, B. M. Chen-Charpentier, Construction of non-
standard finite difference schemes for the SI and SIR epidemic models of fractional
order, Mathematics and Computers in Simulation, 2016, 121, 48-63.
83. J. Cresson, A. Szafran´ska, Discrete and continuous fractional persistence
problems-the positivity property and applications, Communications in Nonlinear
Science and Numerical Simulation, 2017, 44, 424-448.
84. K. Moaddy, I. Hashim, S. Momani, Non-standard finite difference schemes for
solving fractional order Rossler chaotic and hyperchaotic systems, Computers and
Mathematics with Applications, 2011, 62, 1068-1074.
85. A. M. Nagy, N. H. Sweilam, An efficient method for solving fractional Hodgkin-
Huxley model, Physics Letters A, 2014, 378, 1980-1984.
86. H. Su, X. Ding, Dynamics of a nonstandard finite-difference scheme for Mackey-
Glass system, Journal of Mathematical Analysis and Applications, 2008, 344,
932-941.
87. Y. Wang, Dynamics of a nonstandard finite-difference scheme for delay differential
equations with unimodal feedback, Communications in Nonlinear Science and
Numerical Simulation, 2012, 17, 3967-3978.
88. S. M. Garba, A. B. Gumel, A. S. Hassan, J .M. S. Lubuma, Switching from
exact scheme to nonstandard finite difference scheme for linear delay differential
equation, Applied Mathematics and Computation, 2015, 258, 388-403.
89. S. Elaydi, An introduction to Difference Equations, Springer-Verlag, 2005, New
York.
90. V. Sundarapandian, Global asymptotic stability of nonlinear cascade systems,
Applied Mathematics Letters, 2002, 15, 275-277.
170
91. A. Iggidr, M. Bensoubaya, New Results on the Stability of Discrete-Time Systems
and Applications to Control Problems, Journal of Mathematical Analysis and
Applications, 1998, 219, 392-414.
92. A. Gerisch and R. Weiner, The Positivity of Low-Order Explicit Runge-Kutta
Schemes Applied in Splitting Methods, Computers and Mathematics with Applica-
tions, 2003, 45, 53-67.
93. Z. Horváth, Positivity of Runge-Kutta methods and diagonally split Runge-Kutta
methods, Applied Numerical Mathematics, 1998, 28, 306-326.
94. Z. Horváth, On the positivity step size threshold of Runge-Kutta methods, Applied
Numerical Mathematics, 2005, 53, 341-356.
95. J. F. B. M. Kraaijevanger, Contractivity of Runge-Kutta methods, BIT, 1991, 31,
482-528.
96. P. M. Manning, G. F. Margrave, Introduction to non-standard finite-difference
modelling, CREWES Research Report, 2006, 16.
97. J.E. Keymer, P.A. Marquet, J.X. Velasco-Hernandez, S.A. Levin, Extinction thresh-
olds and metapopulation persistence in dynamic landscapes, The American Natu-
ralist, 2000, 156, 478-494.
98. P. Amarasekare, H. Possingham, Patch Dynamics and Metapopulation Theory: the
Case of Successional Species, Journal of Theoretical Biology 2001, 209, 333-344.
99. H.R. Thieme, Convergence results and a Poincar Bendixson trichotomy for asymp-
totically autonomous differential equations, Journal of Mathematical Biology,
1992, 30, 755-763.
100. J. La Salle, S. Lefschetz, Stability by Liapunov’s Direct Method, Academic Press,
1961, New York.
101. J. C. Piqueira, V. O. Araujo, A modified epidemiological model for computer
viruses, Applied Mathematics and Computation, 2009, 213, 355-360.
171
102. P. Szor, The art of computer virus research and defense, Addison-Wesley Educa-
tion Publishers Inc, 2005.
103. A. A. Berryman, The origins and evolution of predator-prey theory, Ecology,
1992, 73, 1530-1535.
104. L. Edelstein-Keshet, Mathematical Models in Biology, Society for Industrial and
Applied Mathematics, 1998, Philadelphia.
105. L. M. Ladino, E. I. Sabogal, J. C. Valverde, General functional response and
recruitment in a predator-prey system with capture on both species, Mathematical
Methods in the Applied Sciences, 2015, 38, 2876-2887.
106. X. Y. Meng, H. F. Hou, H. Xiang, Q. Y. Yin, Stability in a predator-prey model
with Crowley-Martin function and stage structure for prey, Applied Mathematics
and Computation, 2014, 232, 810-819.
107. D. T. Dimitrov and H. V. Kojouharov, Complete mathematical analysis of
predator-prey models with linear prey growth and Beddington-DeAngelis func-
tional response, Applied Mathematics and Computation, 2005, 162, 523-538.
108. G. J. Cooper, J. H. Verner, Some Explicit Runge-Kutta Methods of High Order,
SIAM Journal on Numerical Analysis, 1972, 9, 389-405.
109. C. M. Kribs-Zaleta and J. X. Velasco-Hernández, A simple vaccination model
with multiple endemic states, Mathematical Biosciences, 2000, 164, 183-201.
110. R. E. Mickens, Exact solutions to a finite-difference model of a nonlinear
reaction-advection equation: Implications for numerical analysis, Numerical
Methods for Partial Differential Equations, 1989, 5, 313-325.
111. R. E. Mickens, Dynamic consistency: a fundamental principle for construct-
ing nonstandard finite difference schemes for differential equations, Journal of
Difference Equations and Applications, 2005, 11, 645-653.
112. P. Seibert, R. Suarez, Global stabilization of nonlinear cascade systems, Systems
and Control Letters, 1990, 14, 347-352.
172