Development of nonstandard finite difference methods for some classes of differential equations

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