
Monkeypox (MPX) is a zoonotic illness that is analogous to smallpox. Monkeypox infections have moved across the forests of Central Africa, where they were first discovered, to other parts of the world. It is transmitted by the monkeypox virus, which is a member of the Poxviridae species and belongs to the Orthopoxvirus genus. In this article, the monkeypox virus is investigated using a deterministic mathematical framework within the Atangana-Baleanu fractional derivative that depends on the generalized Mittag-Leffler (GML) kernel. The system's equilibrium conditions are investigated and examined for robustness. The global stability of the endemic equilibrium is addressed using Jacobian matrix techniques and the Routh-Hurwitz threshold. Furthermore, we also identify a criterion wherein the system's disease-free equilibrium is globally asymptotically stable. Also, we employ a new approach by combining the two-step Lagrange polynomial and the fundamental concept of fractional calculus. The numerical simulations for multiple fractional orders reveal that as the fractional order reduces from 1, the virus's transmission declines. The analysis results show that the proposed strategy is successful at reducing the number of occurrences in multiple groups. It is evident that the findings suggest that isolating affected people from the general community can assist in limiting the transmission of pathogens.
Citation: Maysaa Al Qurashi, Saima Rashid, Ahmed M. Alshehri, Fahd Jarad, Farhat Safdar. New numerical dynamics of the fractional monkeypox virus model transmission pertaining to nonsingular kernels[J]. Mathematical Biosciences and Engineering, 2023, 20(1): 402-436. doi: 10.3934/mbe.2023019
[1] | Saulo Orizaga, Maurice Fabien, Michael Millard . Efficient numerical approaches with accelerated graphics processing unit (GPU) computations for Poisson problems and Cahn-Hilliard equations. AIMS Mathematics, 2024, 9(10): 27471-27496. doi: 10.3934/math.20241334 |
[2] | Mingliang Liao, Danxia Wang, Chenhui Zhang, Hongen Jia . The error analysis for the Cahn-Hilliard phase field model of two-phase incompressible flows with variable density. AIMS Mathematics, 2023, 8(12): 31158-31185. doi: 10.3934/math.20231595 |
[3] | Joseph L. Shomberg . Well-posedness and global attractors for a non-isothermal viscous relaxationof nonlocal Cahn-Hilliard equations. AIMS Mathematics, 2016, 1(2): 102-136. doi: 10.3934/Math.2016.2.102 |
[4] | Martin Stoll, Hamdullah Yücel . Symmetric interior penalty Galerkin method for fractional-in-space phase-field equations. AIMS Mathematics, 2018, 3(1): 66-95. doi: 10.3934/Math.2018.1.66 |
[5] | Jean De Dieu Mangoubi, Mayeul Evrard Isseret Goyaud, Daniel Moukoko . Pullback attractor for a nonautonomous parabolic Cahn-Hilliard phase-field system. AIMS Mathematics, 2023, 8(9): 22037-22066. doi: 10.3934/math.20231123 |
[6] | Alain Miranville . The Cahn–Hilliard equation and some of its variants. AIMS Mathematics, 2017, 2(3): 479-544. doi: 10.3934/Math.2017.2.479 |
[7] | Pierluigi Colli, Gianni Gilardi, Jürgen Sprekels . Distributed optimal control of a nonstandard nonlocal phase field system. AIMS Mathematics, 2016, 1(3): 225-260. doi: 10.3934/Math.2016.3.225 |
[8] | Haifeng Zhang, Danxia Wang, Zhili Wang, Hongen Jia . A decoupled finite element method for a modified Cahn-Hilliard-Hele-Shaw system. AIMS Mathematics, 2021, 6(8): 8681-8704. doi: 10.3934/math.2021505 |
[9] | Aymard Christbert Nimi, Daniel Moukoko . Global attractor and exponential attractor for a Parabolic system of Cahn-Hilliard with a proliferation term. AIMS Mathematics, 2020, 5(2): 1383-1399. doi: 10.3934/math.2020095 |
[10] | Hyun Geun Lee . A mass conservative and energy stable scheme for the conservative Allen–Cahn type Ohta–Kawasaki model for diblock copolymers. AIMS Mathematics, 2025, 10(3): 6719-6731. doi: 10.3934/math.2025307 |
Monkeypox (MPX) is a zoonotic illness that is analogous to smallpox. Monkeypox infections have moved across the forests of Central Africa, where they were first discovered, to other parts of the world. It is transmitted by the monkeypox virus, which is a member of the Poxviridae species and belongs to the Orthopoxvirus genus. In this article, the monkeypox virus is investigated using a deterministic mathematical framework within the Atangana-Baleanu fractional derivative that depends on the generalized Mittag-Leffler (GML) kernel. The system's equilibrium conditions are investigated and examined for robustness. The global stability of the endemic equilibrium is addressed using Jacobian matrix techniques and the Routh-Hurwitz threshold. Furthermore, we also identify a criterion wherein the system's disease-free equilibrium is globally asymptotically stable. Also, we employ a new approach by combining the two-step Lagrange polynomial and the fundamental concept of fractional calculus. The numerical simulations for multiple fractional orders reveal that as the fractional order reduces from 1, the virus's transmission declines. The analysis results show that the proposed strategy is successful at reducing the number of occurrences in multiple groups. It is evident that the findings suggest that isolating affected people from the general community can assist in limiting the transmission of pathogens.
In the context of phase field modeling, the equations that are often used to describe phase-separation in micro and nano-structure evolution of materials are in the form of higher order parabolic nonlinear partial differential equations (PDEs). The Cahn Hilliard (CH) equation originally proposed in 1958 to model phase separation in binary alloys [1], has found many other applications in physical processes related to block-copolymers, crystal growth, drug delivery, tumor growth, and image processing [2,3,4,5,6]. For the extended applications of the CH equation, full material separation is prevented leading to micro-phase separation, which is a result of interfacial energy interactions of the materials. For such models, different pattern formations are possible, which are often related to key structural properties of the materials [7].
The classic CH equation with constant mobility is a fourth order nonlinear PDE that derives from an energy functional that penalizes sharp transitions and selects minima related to the two materials in the mixture. There have been a number of numerical methods which have been proposed to solve the CH equation using finite differences, finite elements, finite volume methods, and spectral methods [8,9,10,11,12]. An extended version of the CH equation equation is found in the block copolymer (BCP) equation, which was proposed by Ohta and Kawasaki [13] to model BCP microstructure evolution. The energy functional for the BCP equation is similar to the one used in the CH equation but there is an added term to account for polymer stretching. Another extension for the CH equation is the phase field crystal (PFC) equation proposed by Elder and Grant [14]. Computational approaches for CH, BCP, and PFC equations can be found in [9,10,13,15,16].
Another very important equation in the phase field modeling area is the Functionalized Cahn-Hilliard (FCH) equation [17], which has been used to model micro-phase separation in ternary systems (oil-water-surfactant). The FCH equation is a sixth order parabolic PDE for which explicit methods are not practical since they would require very small steps for numerical stability (h≤Odx6), where h is the time step and dx is the grid spacing in one dimension. On the other hand, fully implicit methods would ensure stability, but can be computationally expensive for problems in more than one dimension. The FCH equation is a more challenging problem compared to the previous mentioned extensions of the CH equation. The challenges are due to the wide variations in space and time scales that need to be captured. The high PDE order makes the time evolution of spatial discretizations very stiff. Some progress has been made in the numerical solutions to the FCH equation [18], but more work needs to be done. Recent numerical approaches have been proposed in [19] using the (scalar auxiliary variable) SAV method. In our approach, we extend the ideas of the bi-harmonic modified (BHM) [20] method and introduce a parameter M3 for the FCH equation that allows for all nonlinear terms in the equation to be computed explicitly. This in turn gives an efficient implementation for the numerical solution at each time step. Numerical examples using the tri-harmonic modified (THM) approach are presented and compared to the work presented in [19] to demonstrate the efficiency, convergence properties, and ease of implementation of the proposed methodology. The purpose of this paper is to provide a numerical approach for the FCH equation that is efficient, accurate, and that is easy to implement with complexity similar to the original convexity-splitting (CS) approach [21] and the BHM approach [20].
In this paper, we utilize the THM approach as the basis for the new algorithms, which are coupled with different linear extrapolation and higher order time discretization techniques. In Section 2, we introduce the main mathematical model for the FCH equation to be considered in the first part of the numerical experiments. In Section 3, we explain how the THM approach can be applied to this version of the FCH equation and we provide a stability analysis of the proposed scheme. In Section 4, we introduce a series of benchmark problems for the FCH equation and we directly compare the performance of the THM approach versus the SAV methodology in terms of accuracy and order of convergence. Section 5 introduces a more complicated version of the FCH equation and we explain the adaptability of the THM approach to such models. In addition, we also show the adaptability of the THM approach with a higher order time stepping technique using an implicit-explicit (IMEX) Runge-Kutta (RK) formulation.
The FCH free energy is given in the form [17],
E(u)=∫Ω12(−ϵ2Δu+W′(u))2−η(ϵ22|∇u|2+W(u))dx | (2.1) |
where W(u) is a symmetric double well function given by W(u)=14(u2−1)2. Here, η>0 is a parameter that represents the properties of amphiphilic phase at the interface (η<0 defines the Cahn-Hilliard-Willmore (CHW) energy [19], which is not considered in this paper), and ϵ is a parameter that characterizes the thickness of the diffusive interface.
The FCH equation is obtained as the H−1 gradient flow of the free energy (2.1)
ut=MΔμ | (2.2) |
μ=ϵ2Δω−W″(u)ω+ηω | (2.3) |
ω=ϵ2Δu−W′(u) | (2.4) |
where μ is the chemical potential, W′(u)=u3−u, W″(u)=3u2−1 and M is the mobility constant. Considering periodic boundary conditions, the scalar solution u(x,t) will evolve into configurations in a way that the FCH energy is always nonincreasing, that is dE(u)dt≤0. Numerical solutions of the FCH must retain this energy decreasing property. We note that this means that the scheme should produce numerical solutions in which the energy is either constant or decreasing but not increasing over time.
The approach that we will consider for this problem relies in the extension of the BHM approach [20]. One of the goals for this formulation is maximum simplicity. We rewrite the Eqs (2.2)–(2.4) as
ut=∇⋅[M∇(ϵ4Δ2u+S(u))], | (3.1) |
S(u)=−ϵ2Δ(W′(u))−W″(u)(ϵ2Δu−W′(u))+η(ϵ2Δu−W′(u)). | (3.2) |
Having S(u) defined by (3.2) allows for us to write a single equation in the form
ut=∇⋅[M∇(ϵ4Δ2u+S(u))]. | (3.3) |
We note that (3.3) is identical to the FCH equation given by Eqs (2.2)–(2.4). The advantage of writing the FCH in this way is that now we can consider a splitting approach by extending the ideas of the BHM method proposed by Bertozzi [20]. Equation (3.3) is expanded by distributing the divergence operator to give
ut=ϵ4∇⋅[M∇Δ2u]+∇⋅[M∇S(u)], | (3.4) |
which now is suitable for introducing a constant M3 and rewriting Eq (3.4) as
ut=M3ϵ4Δ3u+ϵ4∇⋅[(M−M3)∇Δ2u]+∇⋅[M∇S(u)], | (3.5) |
which leads to the form of the FCH equation with a tri-harmonic splitting parameter M3.
We denote the discretized approximations to the solution u by u(xi,tn)≈Un,i, where the time step h is given by tn+1=tn+h. Accounting for the first order approximation of the left-hand side of (3.5) and treating the linear part of right-hand side in (3.5) implicitly gives the following semi-implicit scheme
Un+1−Unh=Ψ(Un+1)+Φ(Un) | (3.6) |
where Ψ(u)=M3ϵ4Δ3u and Φ(u)=(M−M3)ϵ4Δ3u+MΔS(u). We note that Eq (3.5) is the general form of the splitting that accounts for a variable mobility case (M=M(u)), but this case is not considered in the present problem. For problems with variable mobility, the interested reader is referred to [22,23]. For the class of sixth order parabolic problems, the scheme given by (3.6) could also be referred to as the THM method. With regards to splitting approaches, a well-known CS approach was proposed by Eyre [21] and several others have been used in the variable mobility CH equation [20]. In addition and more recently, the PFC equation has been computed with a stabilizing parameter [19]. Equation (3.6) can be understood as a class of splitting methods for phase field models.
We will use Fourier pseudo-spectral methods [24] for our numerical simulations. The CH equation and related phase field models have been implemented with Fourier spectral methods in the past [25,26,27,28,29]. More recently, the FCH was implemented using Fourier spectral method in the work of Wise et al [19] which is one of the main reasons for us to consider the same type of implementation for our paper.
We include additional details about the spatial discretization for (3.6). The two-dimensional discrete Fourier transform of U can be written as
U≈N/2−1∑kx=−N/2N/2−1∑ky=−N/2ˆU(kx,ky,t)exp[i(kxx+kyy)]. |
Consequently, the Fourier transforms of the harmonic, bi-harmonic and tri-harmonic operators can be expressed as ^ΔU=−k2ˆU, ^Δ2U=k4ˆU, and ^Δ3U=−k6ˆU where k2=k2x+k2y. Applying the Fourier spectral method to (3.6), gives the the following expression:
ˆUn+1=ˆUn+h^Φ(Un)1+hM3ϵ4k6, | (3.7) |
which provides an efficient formula for obtaining ˆUn+1, and then the inverse transform is applied to obtain Un+1.
For the stability criteria of the THM method given by (3.6), we will use the approach presented in [30]. Equation (3.3) is of the form ut=G(u), where G(u)=Ψ(u)+Φ(u). Using the expression ^Φ(Un)=^G(Un)−^Ψ(Un) in (3.7) gives the following:
ˆUn+1=ˆUn+h^G(Un)1+hM3ϵ4k6. | (3.8) |
To analyze the conditions for linear stability, we consider highest order terms in G(u) to arrive at the approximation ^G(Un+en)≈^G(Un+en)−C0ϵ4k6^en, which is then used in (3.8) to obtain the amplification factor associated with the growth of the errors in Un. The amplification factor is given by
σ=1−hCoϵ4k61+hM3ϵ4k6. | (3.9) |
The conditions to guarantee stability of the THM scheme require that |σ|<1, so one gets
−1<1−hCoϵ4k61+hM3ϵ4k6<1,−2<−hCoϵ4k61+hM3ϵ4k6<0. |
Given the quantities considered in this formulation, the less-than-zero inequality is satisfied so one works with
−2(1+hM3ϵ4k6)<−hCoϵ4k6 |
which can be further reduced by rearranging the terms to arrive at the following inequality
ϵ4hk6(C0−2M3)<2. |
In order to guarantee stability (|σ|<1), we satisfy the above inequality with the requirement that (C0−2M3)<0 which gives M3>C0/2. In the context of a constant mobility case, we can set C0=M, which gives a minimum value for the splitting parameter in the THM method (3.6). We note that in the analysis performed in [20], a criteria of M1>M (M1 is the splitting parameter for fourth order equations with constant mobility) was obtained for the original BHM approach. Following the approach in [30], we are able to get an improved criteria for this splitting parameter threshold. We will use M3>M/2 for all the simulations presented in this paper.
Remark 1. We emphasize that the THM approach is different than the original CS approach since the THM approach does not require splitting of the energy into convex and concave parts. The similarities with our approach and the CS approach are only in the shared simplicity of the splitting which leads to fast computations. In addition, we also make a note that the THM approach is an extension of the original BHM approach and many benefits in terms of efficiency and ease of implementation are inherited. Finally, in the next section we will present numerical experiments in which we show that the THM approach leads to good numerical results when compared to the SAV formulation.
Our numerical experiments will be presented in two parts: The first part will be related to the simpler version of the FCH equation and second part of the numerical experiments will include the more general and complicated verson of the FCH equation.
Since we are proposing a new THM method, it is important for us to understand how the method compares to the recently developed schemes for the FCH equation using a SAV formulation. We are interested in order of convergence, accuracy and ease of implementation for the THM approach. For the order of convergence and computing errors, there are several different forms such tests can take:
(ⅰ) We can construct a reference solution [9] using a very small time step h for a corresponding simulation final time tf. Then simulation runs using different h values reaching tf can be compared against this reference solution in order to compute errors and order of convergence.
(ⅱ) We can define an exact solution [31] and insert it into the FCH equation which will produce a modified version of the equation. We can then apply our methods to the modified FCH and for different choices of h we are able to compute the error since we have defined our exact solution.
In order to directly compare the performance of our schemes with those schemes presented in [19], in particular their SAV-BDF1 (scalar auxiliary variable backward differentiation formula of order 1) schemes, we will use a numerical test of the type (ⅱ) for the first few test problems. We note that SAV-BDF1 and the THM method presented in this paper are implemented using the same Fourier pseudo-spectral approach so that we can present a fair comparison of their performance.
We can consider an exact solution, whenever available, or accurate reference solution Uref(x,y,t) which is constructed with a small time step size of h=10−7 and a final time configuration of tf. We then compute the numerical solution using the THM scheme with difference choices of h, and compare against the reference solution at tf, and compute the L2 error with the formula
Error(h)=√∫Ω|U(x,y,tf;h)−Uref(x,y,tf)|2dx. | (4.1) |
We first investigate the convergence properties for our numerical approach by solving Eq (3.3) on Ω=[0,L]2 subject to periodic boundary conditions. For simplicity, all test problems will be considered with periodic boundary conditions. We use Fourier pseudo-spectral discretizations in space using N Fourier modes in each space direction [24,32]. We note that other discretizations are certainly possible (finite differences, finite elements, finite volume). The numerical solution is obtained by solving (3.3) with initial condition u0(x,y,t=0)=sinxcosy, parameter values M=1, M3=5, ϵ=0.1, η=ϵ2, L=4π, and with 2562 Fourier modes. The solution is evolved to a final time tf and compared against the exact solution given by
u(x,y,t)=sinxcosycos(t). | (4.2) |
Equation (4.2) is the solution to the forced FCH equation which is obtained by inserting (4.2) into (3.3), resulting in a forcing term that is used to modify the original equation [19]. To be concrete, the form of the forced FCH equation is given by
ut=∇⋅[M∇(ϵ4Δ2u+S(u))]+F(x,y,t), | (4.3) |
where F(x,y,t)=cos(y)sin(x)(1.9404cos(t)+cos(t)3(5.775−5.955cos(2y)+cos(2x) (5.955+16.965cos(2y)))−sin(t)+cos(t)5(−60cos(x)2 cos(y)4sin(x)2+sin(x)4(30cos(y)4−15sin(2y)2))). We compute the errors in the L2 norm. Here, we consider 3 cases for the computation of the forcing term F(x,y,t), the implicit time level (n+1), explicit (n), and semi-implicit (n+1/2). In addition, we compare our errors to the first order scheme SAV-BDF1 introduced in [19] to solve the FCH equation. For a comprehensive review of SAV methods, which were introduced by Shen [12] to efficiently solve gradient flow problems, we recommend the review-paper on SAV methods [33]. Our numerical results using the THM method are presented in Figure 1. It can be seen that the new method produces small errors for all possible evaluations of the forcing term, and, in addition, the errors from THM are found to be smaller than those from the SAV-BDF1 method [19]. Table 1 shows that the THM method reaches first order convergence with no problems and it is capable to perform well when compared to the SAV formulaton. Figure 1 (right) contains the error curves with several choices of the splitting parameter M3. For problems using a stabilization or splitting parameter, it is possible to introduce instability if the parameter is chosen very small (M3≈0, which defines a fully explicit problem). For the current benchmark problem, we tested several cases and found our lowest value of M3 to be M3≈1. For this case, we can see errors increasing for the range of values centered at h=10−2. On the other hand, when M3≥2, solutions remain accurate, numerical instability is suppressed, and the scheme performs with the expected level of accuracy (see Figure 1, right).
h | SAV-BDF1 | Order | THM | Order |
0.02/20 | 7.93e-04 | – | 8.36e-05 | – |
0.02/21 | 4.38e-04 | 0.86 | 4.45e-05 | 0.90 |
0.02/22 | 2.48e-04 | 0.82 | 2.45e-05 | 0.86 |
0.02/23 | 1.40e-04 | 0.82 | 1.37e-05 | 0.83 |
0.02/24 | 7.72e-05 | 0.86 | 7.76e-06 | 0.83 |
0.02/25 | 4.14e-05 | 0.90 | 4.35e-06 | 0.84 |
0.02/26 | 2.16e-05 | 0.94 | 2.39e-06 | 0.86 |
0.02/27 | 1.11e-05 | 0.96 | 1.28e-06 | 0.89 |
0.02/28 | 5.61e-06 | 0.98 | 6.73e-07 | 0.93 |
0.02/29 | 2.82e-06 | 0.99 | 3.46e-07 | 0.96 |
0.02/210 | 1.42e-06 | 0.99 | 1.75e-07 | 0.98 |
0.02/211 | 7.10e-07 | 1.00 | 8.85 e-08 | 1.00 |
To test the performance of the proposed method for longer runs, we considered a benchmark problem used for the FCH equation which can be found in the work done by Wise et al [19] that also traces back to the problems tested in [34]. We solve the FCH equation given by (3.3) with the following initial condition
u(x,y,t=0)=2exp(sinx+siny−2)+2.2exp(−sinx−siny−2)−1, | (4.4) |
with parameters values M=1, M3=5, ϵ=0.18, η=ϵ2, L=2π and with 1282 Fourier modes. The simulation snapshots are taken at values for t=0;and t=4 (see Figure 2). We test the numerical convergence of the method by following the numerical procedure (ⅰ), that is we construct a reference solution for this problem using a small time step h=10−6 using the THM method with a simple linear extrapolation [9] and compute the L2 errors for different choices of h. Figure 3 (left) shows the error curves for several values of M3. Local truncation errors are noticeable in this simulation and they can be spotted for h>10−2. To improve accuracy, we coupled the THM scheme (3.6) with a simple linear extrapolation which is given by
Un+1−Unh=Ψ(Un+1)+Φ(Un),U0:=2Un−Un−1, | (4.5) |
where we have now introduced the need of knowing two consecutive steps initially. However, this is easily accomplished since this is only required to be done during the initialization of the numerical procedure which allows for the computation cost of the linear extrapolation technique to be very efficient. We note that a linear extrapolation technique allows for a better choice of the initial condition and it is known to improve accuracy of numerical solutions [9]. In particular, the work presented in [9] is related to the classic CH equation and the PFC equation. For a detailed formulation using linear extrapolation and nonlinear extrapolation (not considered in this paper) applied to phase field models, we refer the readers to the work presented in [9]. We label the scheme given by (4.5) as the THM scheme with linear extrapolation as THMLX. Computational results given by THMLX allow for errors to be reduced and remain accurate while allowing larger time steps (h≈10−1.5). The error plots for this problem suggests h=0.01 to be a time step that will offer sufficient accuracy for the THM method, which coincides with the step size proposed in [19] for their SAV-BDF1 method. Figure 3 (right) shows the energy decreasing property for THM and THMLX schemes using h=0.01 and M3=5. We also demonstrate the performance of our methods by running longer simulations. Figure 4 displays the extended dynamics for test problem 2 by using the THM method with h=0.001 and M3=5. The energy decreasing property is preserved for the entire duration of the simulation (see Figure 5). We note that the extended dynamics simulation results are in agreement with the results presented in [19].
Our third test problem is aimed at testing our numerical approach with regards to correctly capturing the meandering instability [35]. We set our initial condition to
u(x,y,t=0)={−1,x>sin(y)+L/2+0.34−1,x<sin(y)+L/2−0.340,otherwise. | (4.6) |
with the rest of the parameters the same as in the previous test problem except for η=0.2 and L=4π. We note that this type of initial condition has been used by different authors [19,36] and a smoothing procedure has been implemented in order to avoid the Gibbs phenomena. Details of the smoothing procedure can be found in [19]. It is also worth mentioning that smoothing procedures are standard in the context of spectral methods [24,32]. We treated our initial condition with a smoothing function found in MATLAB. In particular, we used a Gaussian filter on the initial condition with a standard deviation of 2 which provided a smooth enough initial configuration for our schemes to run. The function used was imgaussfilt(u0,σ), where sigma is the standard deviation. The main goal of the smoothing process is that the spectral computations do not initially blow-up due to numerical instability [24]. We implement the THM method using 256×256 Fourier modes with h=0.001 and M3=5. The computed dynamics capturing the meandering instability are presented in Figure 6 matching results in [19]. Corresponding energy evolution with the decreasing property is shown in Figure 7.
We set up our final test problem from this section by using a random initial state, which is connected to the modeling of a bilayer structure [19]. In this case, the solution to the FCH equaton represents an amphiphilic binary mixture undergoing phase separation [37]. We let our initial condition take the form of u(x,y,t=0)=0.5+0.001rand(x,y) and u(x,y,z,t=0)=0.5+0.001rand(x,y,z) for the 2D and 3D cases, respectively. Simulations are presented in Figures 8–11 using the THM method with ϵ=0.1, η=ϵ2, and Ω=[0,4π]d with d=2,3 for 2D and 3D, respectively. It can be seen that our schemes are able to correctly capture the dynamics for the FCH equation [19] in either 2D and 3D and they do so with an ease of implementation and efficiency provided by the THM approach. The energy decreasing property was retained for the entire duration of the simulations (see Figures 9 and 11).
Test problem 5 : pearling instability
Now that we have tested the performance of the THM method by solving the benchmark problems given by test problems 1–4, we introduce a more complicated version of the FCH equation which is given by the following energy [17,19,35]:
˜E(u)=∫Ω12(ϵ2Δu−˜W′(u))2−η1(ϵ22|∇u|2+η2˜W(u))dx | (5.1) |
where ϵ is a parameter associated with thickness of the interface, and η1>0 and η2 are parameters related to the properties of the amphiphilic materials. The H−1 gradient flow for the above functional gives the following:
ut=MΔμ | (5.2) |
μ=(ϵ2Δ−˜W″(u)+η1)ω+(η1−η2)˜W′(u) | (5.3) |
ω=ϵ2Δu−W′(u) | (5.4) |
where ˜W(u) is a double well-potential of unequal well depths given by
˜W(u)=12(u+1)2(12(u−1)2+23τ(u−2)) | (5.5) |
which implies that ˜W′(u)=(u+1)(u+τ)(u−1) and ˜W″(u)=(u+1)(u−1)+2u(u+τ).
The THM approach applied to this version of the FCH equation (5.2) gives the following scheme,
Un+1−Unh=˜Ψ(Un+1)+˜Φ(Un) | (5.6) |
where ˜Ψ(u)=M3ϵ4Δ3u, ˜Φ(u)=(M−M3)ϵ4Δ3u+MΔP(u) and similar to the previous formulation there is an associated term analogous to S(u) which is labeled as P(u) for this version of the FCH equation and it is given by
P(u)=−ϵ2Δ˜W′(u)−ϵ2˜W″(u)Δu+˜W″(u)˜W′(u)+η1ω+(η1−η2)˜W′(u). | (5.7) |
In addition to the THM approach, we also incorporate a higher order time stepping approach based on IMEX-RK schemes for the equation ut=˜Ψ(Un+1)+˜Φ(Un). The THM-IMEX2 schemes read as,
U(1)=Un+h(2−√22˜Ψ(U(1))+2−√22˜Φ(Un)) | (5.8) |
Un+1=Un+h(2−√22˜Ψ(U(n+1)+√22˜Ψ(U(1))+12−√2˜Φ(U(1))+1−√22−√2˜Φ(U(n))). | (5.9) |
We note that several phase field models have been adapted with the IMEX approach in the past in order to increase accuracy. The classic CH equation was implemented with an IMEX scheme in the work presented in [10]. Ceniceros presented a numerical procedure for the variable CH equation using IMEX time stepping [23] technique. More recently a higher order IMEX scheme was used for the PFC equation in [38]. Our current THM-IMEX2 scheme differs from all the ones mentioned since to the best of our knowledge the THM approach has not been implemented before for the FCH equation.
We will test the numerical convergence of the method by following procedure (ⅰ), that is, we construct a reference solution using a small time step h=10−7 and evolve the solution to tf=1.0. An initial condition that illustrates the pearling instability is used in the form of an elliptical annulus which is given by
u(x,y,t=0)={−1,g(x,y)>L/4+0.2−1,g(x,y)<L/4−0.20,otherwise. | (5.10) |
where g(x,y)=√(x−L/2)2+0.5(y−L/2)2. Similar to the case of meandering instability, we filter our initial condition with a Gaussian filter to obtain a smooth initial configuration. To observe the pearling instability, we evolve the filtered initial condition to a final configuration using tf=1 for which the instability is visible on the longer sides of the ellipse (see Figure 12). The rest of the parameters used in this simulation are given by L=4π,M=1,M3=5,ϵ=0.1,η1=1.45ϵ,η2=2ϵ, and τ=0.125.
The THM-BE (which refers to the THM approach coupled with a backward Euler (BE) approximation for the time derivative) and THM-IMEX2 computed errors are shown in Figure 13. The errors are able to reach 1st order accuracy and reach second order accuracy for small enough choices of the time step h. We note that the THM-IMEX2 with M3=5 scheme requires h<0.01 in order to produce stable results. For these range of values of the time step h, the THM-IMEX2 generates smaller errors than THM-BE, which well justifies the formulation. It is common for the splitting parameter in these type of formulations to influence the stability of the numerical results [9,12,19,33,36]. For this reason, we tried M3=10,15 and we noticed that the THM-IMEX2 scheme was able to generate stable results with larger time steps h, but this came at the expense of losing accuracy. Overall, we found that THM-IMEX2 increased the accuracy of our THM method which can be seen in Figure 13. This improvement in accuracy demonstrates that the THM approach can be adapted with a time stepping technique and generate accurate and stable results with the benefits of ease of implementation. Of course, here the method of choice will be based on the preference of the user in which a smaller time step h could be chosen to implement a slightly more efficient THM scheme to generate accurate results. The other option is to consider a few more computations to implement the THM-IMEX2 scheme with the benefits of added accuracy. We now test the performance of the new THM-IMEX2 method for longer simulation runs using the same initial condition. We note that both schemes are very efficient and are able to run in a modest laptop. Simulation plots are presented in Figure 14 where it can be seen that the pearling instability has fully dominated the ellipse by time t=10 but this is followed by a reconnecting of the pearls/dots to then evolve into a configuration that fully reconnects around t=15 and then follows a meandering instability for t≥20. The computed energy for this simulation is presented in Figure 15 where the energy decreasing property is displayed.
Finally, we demonstrate the performance of the new methods by solving the FCH in full 3D space variables by using parameters associated with the different solutions that the FCH equation is intended to model in mixtures of amphiphilic molecules in solvent [17,35,39]. We set our initial condition to u(x,y,z,t=0)=0.65+0.001rand(x,y,z) and run the simulations using THM-IMEX2 on Ω=[0,2π]3 with h=0.001,M=1,M3=5,ϵ=0.1,η1=5ϵ,τ=−0.4 and tf=4. We then vary η2=6ϵ,4ϵ,−ϵ,−3ϵ. Simulation results are shown in Figure 16 in which level sets of the solutions have been labeled for (u=0.4) and (u=0.45) with the colors yellow and purple, respectively. The solution structures are related to a highly branched pore network η2=6ϵ, a pore network η2=4ϵ, a pearled pore and micelle mixture η2=−ϵ, and isolated micelles η2=−3ϵ, which are given from left to right and top to bottom in Figure 16. The simulation results computed with our new proposed approach are in agreement with the results presented in [39].
A THM numerical approach for solving the FCH equation was proposed and the convergence properties of the scheme were investigated using relevant benchmark problems in the phase field community. Our methods retain the energy decreasing property and reach first and second order convergence while producing small errors. In particular, the errors produced using the THM method are found to be smaller than those provided by the SAV formulation [19], by way of a direct comparison. The main advantage of our approach is that the computations are very efficient since all nonlinearities in the equation are computed at the explicit time level. The computation cost for the THM method is similar to the original BHM method, which was originally developed for variable mobility problems. The THM shares a comparable ease of implementation to the original CS method, which was originally developed for the constant mobility CH equation. One additional advantage of the THM method is found on its simplicity. A simple splitting allows for accurate implementation of the FCH equation under the THM approach. We note that the SAV formulation [19,33] requires additional stabilizer parameters for their methods to produce accurate results. In fact, two different SAV formulations had to be derived in order to handle the two versions of the FCH equation in [19] which gives an edge to the THM method in ease of implementation. The THM method requires only one simple splitting on the tri-harmonic term that relies on a single parameter M3 and does not suffer for the required additional steps and reformulation of the energy as the SAV method [12,19,33]. In addition to efficiency, our method is easy to implement and can run in a modest laptop for 2D and 3D implementations in the course of a few hours. In addition to numerical convergence, we tested our method for longer time runs with and without linear extrapolation in which dynamics for the FCH equation were properly displayed by capturing the microstructure evolution of material and the more challenging meandering type of instability. To demonstrate the robustness of the THM method, we applied the method to a more complicated form of the FCH equation and we also coupled our method with an IMEX time stepping technique to improve errors and to capture the pearling instability dynamics. Additionally, we tested the performance of the THM-IMEX2 scheme in full 3D to capture more complex dynamics associated with the more complicated version of the FCH equation. The THM approach provides a simple and efficient implementation for the FCH equation and it produces small errors. We believe that the proposed THM approach could serve as a powerful tool, that is easy to implement, to study not only the FCH equation but also a general class of sixth order phase field models and a class of 4th order problems [30,40], for which the original BHM approach remains applicable. The work presented in this numerical investigation can be viewed as preliminary results and formulations of improved/extended numerical schemes that will lead to further studies and rigorous numerical analysis. Future research work will include theoretical analysis of the energy stability and the convergence analysis of the scheme with constant and variable mobility. We will also study the applicability of the new approach to systems of CH equations in 2D and 3D related to the modeling of biomedical applications.
All authors contributed equally for the completion of the manuscript. All authors have read and approved the final version of the manuscript for publication.
The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.
S. Orizaga is grateful for the support from New Mexico IDeA Networks of Biomedical Research Excellence (NM-INBRE) and from the 2024 NM-INBRE NIRF Faculty Fellowship award. S. Orizaga thanks Professors Snezna Rogelj and Sally Pias for summer funding opportunities. We thank the anonymous reviewers who provided valuable feedback resulting in an improved version of this paper.
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
[1] |
A. W. Rimoin, P. M. Mulembakani, S. C. Johnstonm, J. O. Lloyd Smith, N. K. Kisalu, T. L. Kinkela, et al., Major increase in human monkeypox incidence 30 years after smallpox vaccination campaigns cease in the Democratic Republic of Congo, Proc. Natl. Acad. Sci. USA, 107 (2010), 16262–16267. http://dx.doi.org/10.1073/pnas.1005769107 doi: 10.1073/pnas.1005769107
![]() |
[2] | N. Sklenovská, M. Van Ranst, Emergence of monkeypox as the most important orthopoxvirus infection in humans, Front. Public Health, 241 (2018). https://doi.org/10.3389/fpubh.2018.00241 |
[3] | Singapore Ministry of Health, Europe, US on alert over detection of Monkeypox cases: What is the virus, symptoms and its transmission across the globe. Available from: https://news.knowledia.com/IN/en/articles/europe-us-on-alert-over-detection-of-monkeypox-cases-what-is-the-virus-5. |
[4] | F. Fenner, D. A. Henderson, I. Arita, Z. Jezek, I. D. Ladnyi, Smallpox and its eradication, W. H. O., 1988. |
[5] | R. B. Kennedy, J. M. Lane, D. A. Henderson, G. A. Poland, Smallpox and vaccinia, Vaccines (chapter 32), Amsterdam: Elsevier, (2012), 718–727. |
[6] |
P. E. M. Fine, Z. Jezek, B. Grab, H. Dixon, The transmission potential of monkeypox virus in human populations, Int. J. Epidemiol. , 17 (1988), 643–650. http://dx.doi.org/10.1093/ije/17.3.643 doi: 10.1093/ije/17.3.643
![]() |
[7] |
K. D. Reed, J. W. Melski, M. B. Graham, R. L. Regnery, M. J. Sotir, M. V. Wegner, et al., The detection of monkeypox in humans in the Western Hemisphere, Engl. J. Med. , 350 (2004), 342–350. http://dx.doi.org/10.1056/ NEJMoa032299 doi: 10.1056/NEJMoa032299
![]() |
[8] | M. Roberts, Monkeypox to get a new name, says WHO. Available from: https://en.kataeb.org/articles/monkeypox-to-get-a-new-name-says-who. |
[9] |
L. A. Learned, M. G. Reynolds, D. W. Wassa, Y. Li, V. A. Olson, K. Karem, et al., Extended interhuman transmission of monkeypox in a hospital community in the Republic of the Congo, Am. J. Trop. Med. Hygiene, 73 (2005), 428–434. https://doi.org/10.4269/ajtmh.2005.73.428 doi: 10.4269/ajtmh.2005.73.428
![]() |
[10] |
R. A. Elderfield, S. J. Watson, A. Godlee, W. E. Adamson, C. I. Thompson, J. Dunning, M. Fernandez-Alonso, D. Blumenkrantz, T. Hussell, M. Zambon, Accumulation of human-adapting mutations during circulation of A (H1N1) pdm09 influenza virus in humans in the United Kingdom, J. Virol. , 88 (2014), 13269–13283. https://doi.org/10.1128/JVI.01636-14 doi: 10.1128/JVI.01636-14
![]() |
[11] |
N. C. Elde, S. J. Child, M. T. Eickbush, J. O. Kitzman, K. S. Rogers, J. Shendure, et al., Poxviruses deploy genomic accordions to adapt rapidly against host antiviral defenses, Cell, 150 (2012), 831–841. https://doi.org/10.1016/j.cell.2012.05.049 doi: 10.1016/j.cell.2012.05.049
![]() |
[12] |
R. J. Jackson, A. J. Ramsay, C. D. Christensen, S. Beaton, D. F. Hall, I. A. Ramshaw, Expression of mouse interleukin-4 by a recombinant ectromelia virus suppresses cytolytic lymphocyte responses and overcomes genetic resistance to mousepox. J. Virol. , 75 (2001), 1205–1210. https://doi.org/10.1128/JVI.75.3.1205-1210.2001 doi: 10.1128/JVI.75.3.1205-1210.2001
![]() |
[13] |
S. Bidari, E. E. Goldwyn, Stochastic models of influenza outbreaks on a college campus, Lett. Biomath. , 6 (2019), 1–14. https://doi.org/10.1080/23737867.2019.1618744 doi: 10.1080/23737867.2019.1618744
![]() |
[14] |
J. C. Blackwood, L. M. Childs, An introduction to compartmental modeling for the budding infectious disease modeler. Lett. Biomath. , 5 (2018), 195–221. https://doi.org/10.30707/LiB5.1Blackwood doi: 10.30707/LiB5.1Blackwood
![]() |
[15] | C. Bhunu, S. Mushayabasa, Modelling the transmission dynamics of Pox-like infections, IAENG Int. J. Appl. Math. , 41 (2011), 141–149. |
[16] |
S. Usman, I. I. Adamu, Modeling the transmission dynamics of the monkeypox virus infection with treatment and vaccination interventions, J. Appl. Math. Phy. , 5 (2017), 2335–2353. https://doi.org/10.4236/jamp.2017.512191 doi: 10.4236/jamp.2017.512191
![]() |
[17] |
A. Atangana, D. Baleanu, New fractional derivatives with nonlocal and non-singular kernel: Theory and application to heat transfer model, Therm. Sci. , 20 (2016), 763–769. https://doi.org/10.2298/TSCI160111018A doi: 10.2298/TSCI160111018A
![]() |
[18] |
S. A. Iqbal, M. G. Hafez, Y. M. Chu, C. Park, Dynamical Analysis of nonautonomous RLC circuit with the absence and presence of Atangana-Baleanu fractional derivative, J. Appl. Anal. Comput. 12 (2022), 770–789. https://doi.org/10.11948/20210324 doi: 10.11948/20210324
![]() |
[19] | Y. M. Chu, S. Bashir, M. Ramzan, M. Y. Malik, Model-based comparative study of magnetohydrodynamics unsteady hybrid nanofluid flow between two infinite parallel plates with particle shape effects, Math. Meth. Appl. Sci., 2022. https://doi.org/10.1002/mma.8234 |
[20] |
W. M. Qian, H. H. Chu, M. K. Wang, Y. M. Chu, Sharp inequalities for the Toader mean of order −1 in terms of other bivariate means, J. Math. Inequal. , 16 (2022), 127–141. https://doi.org/10.7153/jmi-2022-16-10 doi: 10.7153/jmi-2022-16-10
![]() |
[21] |
T. H. Zhao, H. H. Chu, Y. M. Chu, Optimal Lehmer mean bounds for the nth power-type Toader mean of n=−1,1,3, J. Math. Inequal. , 16 (2022), 157–168. https://doi.org/10.7153/jmi-2022-16-12 doi: 10.7153/jmi-2022-16-12
![]() |
[22] | M. Caputo, M. Fabrizio, A new definition of fractional derivative without singular kernel, Progr. Fract. Differ. Appl. , 1 (2015), 1–13. |
[23] |
J. F. Gómez-Aguilar, H. Yéppez-Marténez, C. Calderón-Ramón, I. Cruz-Orduña, R. F. Escobar-Jiménez, V. H. Olivares-Peregrino, Modeling of a mass-spring-damper system by fractional derivatives with and without a singular kernel, Entropy, 17 (2015), 6289–6303. https://doi.org/10.3390/e17096289 doi: 10.3390/e17096289
![]() |
[24] | F. Z. Wang, M. N. Khan, I. Ahmad, H. Ahmad, H. Abu-Zinadah, Y. M. Chu, Numerical solution of traveling waves in chemical kinetics: time-fractional fishers equations, Fractals, 30 (2022). https://doi.org/10.1142/S0218348X22400515 |
[25] | S. Rashid, E. I. Abouelmagd, S. Sultana, Y. M. Chu, New developments in weighted n-fold type inequalities via discrete generalized ˆh-proportional fractional operators, Fractals, 30 (2022). https://doi.org/10.1142/S0218348X22400564 |
[26] |
S. Rashid, B. Kanwal, A. G. Ahmad, E. Bonyah, S. K. Elagan, Novel numerical estimates of the pneumonia and meningitis epidemic model via the nonsingular kernel with optimal analysis, Complexity, 2022 (2022), 1–25. https://doi.org/10.1155/2022/4717663 doi: 10.1155/2022/4717663
![]() |
[27] | S. W. Yao, S. Rashid, M. Inc, E. Elattar, On fuzzy numerical model dealing with the control of glucose in insulin therapies for diabetes via nonsingular kernel in the fuzzy sense, AIMS Math., 7 (2022). https://doi.org/10.3934/math.2022987 |
[28] |
S. Rashid, F. Jarad, A. K. Alsharidi, Numerical investigation of fractional-order cholera epidemic model with transmission dynamics via fractal-fractional operator technique, Chaos Solitons Fractals, 162 (2022), 112477. https://doi.org/10.1016/j.chaos.2022.112477 doi: 10.1016/j.chaos.2022.112477
![]() |
[29] |
F. Jin, Z. S. Qian, Y. M. Chu, M. ur Rahman, On nonlinear evolution model for drinking behavior under Caputo-Fabrizio derivative, J. Appl. Anal. Comput. 12 (2022), 790–806. https://doi.org/10.11948/20210357 doi: 10.11948/20210357
![]() |
[30] |
S. Rashid, M. K. Iqbal, A. M. Alshehri, F. Jarad, R. Ashraf, A comprehensive analysis of the stochastic fractal-fractional tuberculosis model via Mittag-Leffler kernel and white noise, Results Phy. . 39 (2022), 105764. https://doi.org/10.1016/j.rinp.2022.105764 doi: 10.1016/j.rinp.2022.105764
![]() |
[31] |
M. Al-Qurashi, S. Rashid, F. Jarad, A computational study of a stochastic fractal-fractional hepatitis B virus infection incorporating delayed immune reactions via the exponential decay, Math. Biosci. Eng, 19 (2022), 12950–12980. https://doi.org/10.3934/mbe.2022605 doi: 10.3934/mbe.2022605
![]() |
[32] |
S. Rashid, A. Khalid, S. Sultana, F. Jarad, K. M. Abulanja, Y. S. Hamed, Novel numerical investigation of the fractional oncolytic effectiveness model with M1 virus via generalized fractional derivative with optimal criterion, Results Phy. , 37 (2022), 105553. https://doi.org/10.1016/j.rinp.2022.105553 doi: 10.1016/j.rinp.2022.105553
![]() |
[33] |
S. Rashid, B. Kanwal, F. Jarad, S. K. Elagan, A peculiar application of the fractal-fractional derivative in the dynamics of a nonlinear scabies model, Results Phys. , 38 (2022), 105634. https://doi.org/10.1016/j.rinp.2022.105634 doi: 10.1016/j.rinp.2022.105634
![]() |
[34] |
S. Rashid, Y. G. Sánchez, J. Singh, Kh. M. Abualnaja, Novel analysis of nonlinear dynamics of a fractional model for tuberculosis disease via the generalized Caputo fractional derivative operator (case study of Nigeria), AIMS Math. , 7 (2022), 10096–10121. https://doi.org/10.3934/math.2022562 doi: 10.3934/math.2022562
![]() |
[35] |
S. Rashid, F. Jarad, A. G. Ahmad, Kh. M. Abualnaja, New numerical dynamics of the heroin epidemic model using a fractional derivative with Mittag-Leffler kernel and consequences for control mechanisms, Results Phys. , 35 (2022), 105304. https://doi.org/10.1016/j.rinp.2022.105304 doi: 10.1016/j.rinp.2022.105304
![]() |
[36] |
S. N. Hajiseyedazizi, M. E. Samei, J. Alzabut, Y. M. Chu, On multi-step methods for singular fractional q-integro-differential equations, Open Math. 19 (2021), 1378–1405. https://doi.org/10.1515/math-2021-0093 doi: 10.1515/math-2021-0093
![]() |
[37] | S. Rashid, E. I. Abouelmagd, A. Khalid, F. B. Farooq, Y. M. Chu, Some recent developments on dynamical ℏ-discrete fractional type inequalities in the frame of nonsingular and nonlocal kernels, Fractals, 30 (2022). https://doi.org/10.1142/S0218348X22401107 |
[38] | Y. M. Chu, U. Nazir, M. Sohail, M. M. Selim, J. R. Lee, Enhancement in thermal energy and solute particles using hybrid nanoparticles by engaging activation energy and chemical reaction over a parabolic surface via finite element approach, Fractal Fract. 5 (2021), 119. https://doi.org/10.3390/fractalfract5030119 |
[39] |
T. Abdeljawad, D. Baleanu, Integration by parts and its applications of a new nonlocal fractional derivative with Mittag-Leffler nonsingular kernel, J. Nonlinear Sci. Appl. , 10 (2017), 1098–1107. http://dx.doi.org/10.22436/jnsa.010.03.20 doi: 10.22436/jnsa.010.03.20
![]() |
[40] |
J. Singh, D. Kumar, D. Baleanu, On the analysis of fractional diabetes model with exponential law, Adv. Differ. Equ. , 2018 (2018), 231. https://doi.org/10.1186/s13662-018-1680-1 doi: 10.1186/s13662-018-1680-1
![]() |
[41] |
Z. M. Odibat, N. T. Shawagfeh, Generalized Taylor's formula, Appl. Math. Comput. 186 (2007), 286–293. https://doi.org/10.1016/j.amc.2006.07.102 doi: 10.1016/j.amc.2006.07.102
![]() |
[42] |
S. V. Bankuru, S. Kossol, W. Hou, P. Mahmoudi, J. Rychtáˆr, D. Taylor, A game-theoretic model of Monkeypox to assess vaccination strategies, PeerJ, 8 (2020), 9272. http://doi.org/10.7717/peerj.927 doi: 10.7717/peerj.927
![]() |
[43] | O. J. Peter, S. Kumar, N. Kumari, F. A. Oguntolu, K. Oshinubi, R. Musa, Transmission dynamics of Monkeypox virus: A mathematical modelling approach, Model. Earth Sys. Envir., 2021. https://doi.org/10.1007/s40808-021-01313-2 |
[44] |
M. Toufik, A. Atangana, New numerical approximation of fractional derivative with non-local and non-singular kernel: Application to chaotic models, Eur. Phys. J. Plus, 132 (2017), 444. https://doi.org/10.1140/epjp/i2017-11717-0 doi: 10.1140/epjp/i2017-11717-0
![]() |
[45] |
C. Bhunu, W. Garira, G. Magombedze, Mathematical analysis of a two strain hiv/aids model with antiretroviral treatment. Acta Biotheor. , 57 (2009), 361–381. https://doi.org/10.1007/s10441-009-9080-2 doi: 10.1007/s10441-009-9080-2
![]() |
[46] | C. Bhunu, S. Mushayabasa, Modelling the transmission dynamics of pox-like infections, IAENG Int. J. Appl. Math. , 41 (2011), 141–149. |
[47] |
M. R. Odom, R. Curtis Hendrickson, E. J. Lefkowitz, Poxvirus protein evolution: Family wide assessment of possible horizontal gene transfer events, Virus Res. , 144 (2009), 233–249. https://doi.org/10.1016/j.virusres.2009.05.006 doi: 10.1016/j.virusres.2009.05.006
![]() |
[48] |
Y. Lia, Y. Chen, I. Podlubny, Mittag–Leffler stability of fractional order nonlinear dynamic systems, Automatica, 45 (2009), 1965–1969. https://doi.org/10.1016/j.automatica.2009.04.003 doi: 10.1016/j.automatica.2009.04.003
![]() |
[49] |
S. Somma, N. Akinwande, U. Chado, A mathematical model of monkey pox virus transmission dynamics, IFE J. Sci. , 21 (2019), 195–204. https://doi.org/10.4314/ijs.v21i1.17 doi: 10.4314/ijs.v21i1.17
![]() |
[50] |
O. Diekmann, J. Heesterbeek, M. G. Roberts, The construction of next-generation matrices for compartmental epidemic models, J. R. Soc. Interface, 7 (2010), 873–885. https://doi.org/10.1098/rsif.2009.0386 doi: 10.1098/rsif.2009.0386
![]() |
[51] | P. van den Driessche, J. Watmough, Further notes on the basic reproduction number, Springer, Berlin, Heidelberg, (2008), 159–178. |
[52] |
C. Castillo-Chavez, B. Song, Dynamical models of tuberculosis and their applications, Math. Biosci. Eng. , 1 (2004), 361. https://doi.org/10.3934/mbe.2004.1.361 doi: 10.3934/mbe.2004.1.361
![]() |
[53] |
E. Ahmed, A. M. A. El-Sayed, H. A. A. El-Saka, On some routh–hurwitz conditions for fractional order differential equations and their applications in lorenz, rössler, chua and chen systems, Phys. Lett. A, 358 (2006), 1–4. https://doi.org/10.1016/j.physleta.2006.04.087 doi: 10.1016/j.physleta.2006.04.087
![]() |
[54] |
C. T. Deressa, G. F. Duressa, Analysis of Atangana–Baleanu fractional-order SEAIR epidemic model with optimal control, Adv. Diff. Equ. , 2021 (2021), 174. https://doi.org/10.1186/s13662-021-03334-8 doi: 10.1186/s13662-021-03334-8
![]() |
[55] |
B. Wu, F. Fu, L. Wang, Imperfect vaccine aggravates the long-standing dilemma of voluntary vaccination, PLOS ONE, 6 (2011), 20577. https://doi.org/10.1371/journal.pone.0020577 doi: 10.1371/journal.pone.0020577
![]() |
[56] | L. Khodakevich, Z. Ježek, D. Messinger, Monkeypox virus: Ecology and public health significance, Bull. W. H. O. , 66 (1988), 747–752. |
[57] |
J. Kobe, N. Pritchard, Z. Short, I. V. Erovenko, J. Rychtář, J. T. Rowel, A game theoretic model of cholera with optimal personal protection strategies. Bull. Math. Bio., 80 (2018), 2580–2599. https://doi.org/10.1007/s11538-018-0476-5 doi: 10.1007/s11538-018-0476-5
![]() |
1. | Saulo Orizaga, Maurice Fabien, Michael Millard, Efficient numerical approaches with accelerated graphics processing unit (GPU) computations for Poisson problems and Cahn-Hilliard equations, 2024, 9, 2473-6988, 27471, 10.3934/math.20241334 |
h | SAV-BDF1 | Order | THM | Order |
0.02/20 | 7.93e-04 | – | 8.36e-05 | – |
0.02/21 | 4.38e-04 | 0.86 | 4.45e-05 | 0.90 |
0.02/22 | 2.48e-04 | 0.82 | 2.45e-05 | 0.86 |
0.02/23 | 1.40e-04 | 0.82 | 1.37e-05 | 0.83 |
0.02/24 | 7.72e-05 | 0.86 | 7.76e-06 | 0.83 |
0.02/25 | 4.14e-05 | 0.90 | 4.35e-06 | 0.84 |
0.02/26 | 2.16e-05 | 0.94 | 2.39e-06 | 0.86 |
0.02/27 | 1.11e-05 | 0.96 | 1.28e-06 | 0.89 |
0.02/28 | 5.61e-06 | 0.98 | 6.73e-07 | 0.93 |
0.02/29 | 2.82e-06 | 0.99 | 3.46e-07 | 0.96 |
0.02/210 | 1.42e-06 | 0.99 | 1.75e-07 | 0.98 |
0.02/211 | 7.10e-07 | 1.00 | 8.85 e-08 | 1.00 |
h | SAV-BDF1 | Order | THM | Order |
0.02/20 | 7.93e-04 | – | 8.36e-05 | – |
0.02/21 | 4.38e-04 | 0.86 | 4.45e-05 | 0.90 |
0.02/22 | 2.48e-04 | 0.82 | 2.45e-05 | 0.86 |
0.02/23 | 1.40e-04 | 0.82 | 1.37e-05 | 0.83 |
0.02/24 | 7.72e-05 | 0.86 | 7.76e-06 | 0.83 |
0.02/25 | 4.14e-05 | 0.90 | 4.35e-06 | 0.84 |
0.02/26 | 2.16e-05 | 0.94 | 2.39e-06 | 0.86 |
0.02/27 | 1.11e-05 | 0.96 | 1.28e-06 | 0.89 |
0.02/28 | 5.61e-06 | 0.98 | 6.73e-07 | 0.93 |
0.02/29 | 2.82e-06 | 0.99 | 3.46e-07 | 0.96 |
0.02/210 | 1.42e-06 | 0.99 | 1.75e-07 | 0.98 |
0.02/211 | 7.10e-07 | 1.00 | 8.85 e-08 | 1.00 |