EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS Vol. 17, No. 4, 2024, 2676-2691 ISSN 1307-5543 – ejpam.com Published by New York Business Global Using Modified Fractional Euler Formula for Solving the Fractional Smoking Model Mohamed M. Khader1,2 1 Department of Mathematics and Statistics, College of Science, Imam Mohammad Ibn Saud Islamic University (IMSIU), Riyadh, Saudi Arabia 2 Department of Mathematics, Faculty of Science, Benha University, Benha, Egypt Abstract. Worldwide, smoking is a common social practice, especially in places like schools and on important occasions. World Health Organization (WHO) states that smoking is the third leading cause of death worldwide and the most significant avoidable reason for disease. So, this work is devoted to giving a numerical simulation of the smoking model in its fractional form (Liouville-Caputo sense). The novel numerical scheme used is the modified fractional Euler method (MFEM). We compare the given technique with the traditional fractional Euler method. In the hope of providing some recommendations to reduce the risks of this bad behavior, the effect of some parameters and external factors affecting the solution behavior of this proposed mathematical model was studied, including the recruitment rate (due to immigration or birth) and the smoking cessation rate. The results confirm the implemented scheme is a straightforward and efficient tool for obtaining solutions to these problems. 2020 Mathematics Subject Classifications: 26A33, 34K37, 65M60, 65N12 Key Words and Phrases: Fractional smoking model, Liouville-Caputo fractional derivative, Modified fractional Euler method 1. Introduction According to a report by the WHO, many smokers pass away in their prime years. More than 5 million deaths globally occur each year as a result of smoking’s effects on various body processes; by 2030, this number may increase to 8 million ([14], [17]). Smokers are 70% more likely to experience a heart attack compared to non-smokers. Lung cancer occurs 10% more frequently in smokers than in non-smokers. Smokers often live 10 to 13 years less than others do. Researchers work to extend people’s lives to reduce smoking. Many researchers have attempted to investigate a variety of efficient smoking models to give the finest description of the phenomena of cigarette smoking. By categorizing the entire population into three groups-potential smokers, chain smokers, and permanently abstinent smokers. The smoking model was expressed in mathematical equations and was DOI: https://doi.org/10.29020/nybg.ejpam.v17i4.5366 Email address: mmkhader@imamu.edu.sa (M. M. Khader) https://www.ejpam.com 2676 Copyright: © 2024 The Author(s). (CC BY-NC 4.0) M. M. Khader / Eur. J. Pure Appl. Math, 17 (4) (2024), 2676-2691 2677 originally developed in 1997 [7]. More recently, in 2007, a survey conducted in Korea by Ham [8] recorded the various stages and practices of smoking among students. In this work, the creation of an integer-order model with a dynamic interaction and a new class of infrequent smokers was described as an expansion of the model. The smoking models have also been provided by several other authors in integer and fractional orders ([3]-[13]). The numerical evaluation of fractional-order delay differential equations is done us- ing a spectral collocation-based approximate technique [11]. The authors used fractional operators namely, Caputo-Fabrizio and Liouville-Caputo derivatives. The collocation ap- proach uses Narayana polynomials and their generalization forms. Using convergence analysis in a weighted L2 norm, they found an upper bound on the Narayana polynomi- als’ series expansion form. The present matrix collocation solves three fractional operator and fractional order test situations to demonstrate its performance. In [10], the authors solved fractional-order differential equations with singularity and strong nonlinearity for electro-hydrodynamic flow in a circular cylindrical conduit numerically. By using the quasi-linearization technique, they obtained a set of linearized equations from the nonlin- ear model. Using the generalized shifted airfoil polynomials of the second kind (SAPSK) and appropriate collocation points as SAPSK roots, they created an iterative algebraic system of linear equations. The error analysis and convergence are established in L2, and L∞ norms. The novelty in the current research is the attempt to reach numerical simulations (through a good numerical method) to study the important system under consideration which is of interest to many researchers, so we presented by shedding some light on the convergence and calculating the resulting error as well as comparisons with the same method but in a less accurate case and finally the effect of the parameters in the system on the behavior of the solution to provide recommendations that can be used by those interested in studying this model medically or industrially. This work is devoted to giving a numerical solution for the fractional smoking model. This is by using the MFEM see [9]. The outline of the paper is given as follows: Section 2, presents some basic concepts of fractional calculus. Section 3, describes the smoking system in its fractional form. Section 4, introduces the modified fractional Euler method. Section 5, solves the fractional-order smoking model. Section 6, presents a numerical simulation for the model under study. Section 7, gives the conclusions and discussions. 2. Basic concepts of fractional calculus To improve the system’s ability to describe memory and global correlation, fractional derivative models have employed power-law memory kernels to give the fractional deriva- tives. The most important definition that is utilized in the creation of fractional calculus theory is this one. While modeling certain real-world problems, the Riemann-Liouville formulation has some restrictions [11]. On the other hand, these issues were addressed in the creation of the Liouville-Caputo definition. Here, we will present some concepts of the fractional calculus ([2]-[12]). Definition 1. The fractional Riemann-Liouville integral of a function Ψ(τ) of order β ∈ M. M. Khader / Eur. J. Pure Appl. Math, 17 (4) (2024), 2676-2691 2678 (0, 1) is defined as: IβΨ(τ) = 1 Γ(β) ∫ τ 0 (τ − s)β−1Ψ(s)ds, t > 0. Definition 2. The Liouville-Caputo fractional derivative Dβ of a function Ψ(τ) of order β ∈ (m− 1,m], m ∈ N is formulated by: DβΨ(τ) = 1 Γ(m− β) ∫ τ 0 (τ − s)m−β−1Ψ(m)(s)ds. The relation between the fractional derivatives and fractional integral for their signif- icance is stated as: IβDβΨ(t) = Ψ(t)− m∑ j=1 Ψ(j)(0+) tj j! , t > 0, β ∈ (m− 1,m]. 3. Description of the smoking system Given the importance of mathematical modeling as a good tool for pandemic grasp in recent decades, we may use it to stop tobacco smoking from spreading. A generic model is the susceptible-exposed-infected-recovered (SEIR). In this section, to gain a deeper understanding of the qualitative analysis and the numerical iterative analysis utilized to solve it, the updated version of the smoking model will be considered as follows ([1], [4]): Ṡ1(t) = κ− ε1 S1(t)S2(t) + ω S4(t)− µ S1(t), Ṡ2(t) = ε1 S1(t)S2(t)− ε2 S2(t)S3(t)− (κ1 + µ) S2(t), Ṡ3(t) = ε2 S2(t)S3(t)− (φ+ κ2 + µ) S3(t), Ṡ4(t) = φS3(t)− (λ+ µ+ ω) S4(t), Ṡ5(t) = λS4(t)− µS5(t), (1) where Ṡi(t) = dS dt , i = 1, 2, 3, 4, 5. With the initial conditions Si(0) = S0i . The whole population is split into five classes in this model: At a time t, the susceptible smokers (S1, S2, S3, S4, and S5) represent respectively, the snuffing class, irregular smokers, regular smokers, and quitters. The parameters (κ, φ, λ, µ, ω, ε1, ε2, κ1, κ2) utilized in the system (1) are explained and defined as in the following Table 1: M. M. Khader / Eur. J. Pure Appl. Math, 17 (4) (2024), 2676-2691 2679 Table 1: Description of the included parameters of the system (1). Symbol Description S1 The susceptible smokers S2 The snuffing (ingestion) class S3 The irregular smokers S4 The regular smokers S5 The quit smokers κ The rate of recruitment (due to migration or birth) φ The percentage of infrequent smokers who start smoking regularly λ, µ, ω The rate of natural death, the rate of departure, and the pace of recuperation ε1 The proportion of the population that moves from being vulnerable to being snuffers ε2 Snuffing rates start to fluctuate among smokers κ1 The percentage of smoking-related deaths class κ2 The percentage of smoking-related deaths among the snuffing class 3.1. Liouville-Caputo-fractional smoking system The fractional derivative increases the opportunity to study and evaluate the real phenomenon in terms of fractional, which is why a lot of academics have focused on this subject ([15], [16]). We reinterpret the smoking system in terms of fractional derivatives of order α ∈ (0, 1] as follows: DαS1(t) = κ− ε1 S1(t)S2(t) + ω S4(t)− µ S1(t), DαS2(t) = ε1 S1(t)S2(t)− ε2 S2(t)S3(t)− (κ1 + µ) S2(t), DαS3(t) = ε2 S2(t)S3(t)− (φ+ κ2 + µ) S3(t)), DαS4(t) = φS3(t)− (λ+ µ+ ω) S4(t), DαS5(t) = λS4(t)− µS5(t)). (2) Throughout the examination of the model presented in (2), the subsequent approaches will be employed. We will provide the necessary definitions in the next sections. The advanced portions will then contain numerical results, simulations, and a conclusion. This study’s primary goal is to examine the model under the fractional derivative. 4. Modified fractional Euler method Theorem 1. For the function Ψ(t) which satisfies the condition DkαΨ(t) ∈ C(0, a], for k = 0, 1, ..., n+ 1, where 0 < α ≤ 1. Then the generalization of Taylor’s formula that involves Liouville-Caputo fractional deriva- tives can be defined as follows [18]: Ψ(t) = n∑ i=0 tiα Γ(iα+ 1) DiαΨ(0+) + D(n+1)αΨ(ξ) Γ((n+ 1)α+ 1) t(n+1)α, 0 ≤ ξ ≤ t, ∀ t ∈ (0, a]. (3) M. M. Khader / Eur. J. Pure Appl. Math, 17 (4) (2024), 2676-2691 2680 The classical Taylor’s formula [5] can be obtained from this formula (3) at α = 1. Now, we introduce the generalized Euler’s method (GEM) which was derived by Zaid and Momani [18]. Consider the following general form of the initial value problem (IVP) in the interval [0, T ] with Liouville-Caputo derivatives of order 0 < α ≤ 1 [18]: Dαz(t) = g(t, z(t)), z(0) = z0. (4) To get the required scheme of the GEM, we follow the following steps: (i) Divide the interval [0, T ] into n subintervals [tj , tj+1] of equal width (step size) h = T/n with nodes tj+1 = tj + h, for j = 0, 1, ..., n. (ii) Assume that Dℓ αz(t), ℓ = 0, 1, 2 are continuous on [0, T ]. (iii) Use the generalized Taylor’s formula (3) to expand z(t) about t = t0 = 0 as follows: z(t) = z(t0) + Dαz(t0) Γ(α+ 1) tα + D2αz(c1) Γ(2α+ 1) t2α, 0 < c1 < T. (5) (iv) Substitute by Dαz(t0) = g(t0, z(t0)) and h = t1 − t0 in equation (5), to get the following expression for z(t1): z(t1) = z(t0) + g(t0, z(t0)) hα Γ(α+ 1) +D2αz(c1) h2α Γ(2α+ 1) . (v) If we choose h small enough, then we may neglect the second-order term (involving h2α) in the previous expansion and get the following approximation: z(t1) = z(t0) + g(t0, z(t0)) hα Γ(α+ 1) . (vi) Repeat the above process to generate a sequence of points zj = z(tj) that approxi- mates the solution z(t) at the nodes tj , j = 0, 1, 2, ... . The general formula for the GEM is given as follows: z(tj+1) = z(tj) + g(tj , z(tj)) hα Γ(α+ 1) , j = 0, 1, ..., n− 1. (6) (vii) Apply the developed a new modification in [6], to obtain the general form of the scheme of the MFEM for solving fractional IVP (4) as follows: z(tj+1) = z(tj)+ hα Γ(α+ 1) g ( tj + hα 2Γ(α+ 1) , z(tj) + hα 2Γ(α+ 1) g(tj , z(tj) ) . j = 0, 1, ..., n−1. (7) M. M. Khader / Eur. J. Pure Appl. Math, 17 (4) (2024), 2676-2691 2681 The upper bound of the error generated by the given numerical scheme (7) will be esti- mated through the next theorem. Theorem 2. [6] Suppose that g is a continuous real-valued function and satisfies the following Lipschitz condition: |g(t, z1)− g(t, z2)| ≤ κ|z1 − z2|, κ ∈ Ω = [a, b]× R. Suppose that a constant δ exists with |Dnαz(t)| ≤ δ, ∀ a ≤ t ≤ b. Then, we have: |z(tk)− zk| ≤ γ (ek θ − 1), k = 0, 1, 2, ..., n, (8) where γ and θ are defined as follows: γ = 2Γ2(α+ 1)h2α δ Γ(2α+ 1)(2Γ(α+ 1)hα + h2ακ) , θ = 2Γ(α+ 1)hα + h2ακ 2Γ2(α+ 1) . Proof. The detail of the proof of this theorem can be found in [6]. 5. Solving fractional-order smoking model In this section, we obtain numerical solutions for the smoking model in its fractional form using the method proposed in the previous section (MFEM). For this purpose, we can revise this model (2) again as follows: DαSr(t) = gr(S1,S2, S3,S4,S5, t), r = 1, 2, 3, 4, 5, (9) where g1(S1, S2,S3,S4, S5, t) = κ− ε1 S1(t)S2(t) + ω S4(t)− µS1(t), g2(S1, S2,S3,S4, S5, t) = ε1 S1(t)S2(t)− ε2 S2(t)S3(t)− (κ1 + µ) S2(t), g3(S1, S2, S3,S4, S5, t) = ε2 S2(t)S3(t)− (φ+ κ2 + µ)S3(t), g4(S1,S2, S3,S4, S5, t) = φS3(t)− (λ+ µ+ ω)S4(t), g5(S1,S2, S3,S4,S5, t) = λ S4(t)− µ S5(t). (10) We assume that DkαSr(t), k = 0, 1, 2 are continuous on (0, T ]. Then we can generate the numerical solution Sr(tk), r = 1, 2, 3, 4, 5 of the system (9) by applying the derived numerical scheme (7) with the following iteration form: Sr(tj+1) = Sr(tj) + hα Γ(α+ 1) gr ( tj + hα 2Γ(α+ 1) ,Sr(tj) + hα 2Γ(α+ 1) . gr(tj , S1(tj),S2(tj),S3(tj),S4(tj), S5(tj)) ) , r = 1, 2, 3, 4, 5, (11) where gr, r = 1, 2, 3, 4, 5 are defined in (10). M. M. Khader / Eur. J. Pure Appl. Math, 17 (4) (2024), 2676-2691 2682 6. Computational results In this section, we give a numerical simulation of the proposed model (2) in the interval [0, 30] with varying values of the initial solution and fractional order α, and the important two parameters µ, κ; with step size h = 0.3 and we will be able to confirm the accuracy and persistency of the given technique. However, we use the same values in all figures: κ = 0.1, ε1 = ω = κ1 = κ2 = 0.003, ε2 = µ = 0.002, φ = λ = 0.05. The following two cases of the initial solutions are considered: (i) Case 1: S01 = 25, S02 = 20, S03 = 15, S04 = 10, S05 = 5; (ii) Case 2: S01 = 65, S02 = 55, S03 = 45, S04 = 25, S05 = 15. To assess the accuracy and efficiency of the proposed scheme, we also provide a comparison between the results produced by the given scheme and those obtained by applying the generalized Euler method. Labels of Figures. No. Description Used values Initial conditions 1 Numerical solution α = 1.0, 0.9, 0.8, 0.7 Case 1 2 Numerical solution α = 1.0, 0.9, 0.8, 0.7 Case 2 3 Comparison between the method and GEM α = 0.95 Case 1 4 Comparison between the method and GEM α = 0.95 Case 2 5 The impact of the system’s natural death rate µ = (3, 4, 11, 14)× 10−4 Case 1 6 The impact of the system’s recruitment rate κ = 0.1, 0.3, 0.5, 0.7 Case 1 Figures 1 through 6 present the numerical findings for the model under study that were produced by using the suggested technique. Based on these findings, we can conclude that the given technique is appropriate for solving the current model in its Liouville-Caputo fractional form, as the behavior of the numerical solution obtained by applying the proposed method depends on the values of α, h, µ, and κ. Furthermore, the efficiency of the process and the outcomes are signif- icantly enhanced by the suggested methodology. In Figures 5 and 6, respectively, the effect of the recruitment rate (due to immigration or birth), and the smoking cessation rate on the behavior and dynamics of the solutions was studied. From Figure 5, we find that indeed, increasing the values of the departure rate leads to a decrease in the values of each of the four components S2(t) to S5(t) without affecting the component S1(t), and this is largely consistent with the natural meaning of the model under study, especially from a medical and therapeutic perspective. Likewise, from Figure 6 we found that in- deed, increasing the values of the recruitment rate (due to migration or birth) leads to an increase in the values of each of the four components S1(t) to S4(t) without affecting the component S5(t), and this is largely consistent with the natural meaning of the system under study, especially from a medical and therapeutic perspective. M. M. Khader / Eur. J. Pure Appl. Math, 17 (4) (2024), 2676-2691 2683 Fig. 1: The solution Si(t), i = 1, 2, 3, 4, 5 for some α with small initial values. M. M. Khader / Eur. J. Pure Appl. Math, 17 (4) (2024), 2676-2691 2684 Fig. 2: The solution Si(t), i = 1, 2, 3, 4, 5 for some α with large initial values. M. M. Khader / Eur. J. Pure Appl. Math, 17 (4) (2024), 2676-2691 2685 Fig. 3: The Si(t), i = 1, 2, 3, 4, 5 by the present and GEM methods α = 1 with small initial values. In addition, to strongly prove and confirm the effectiveness of the given technique, we present a comparison between the MFEM with the generalized fractional Euler method. This comparison is done through Table 2 at some different iterations of the procedure with α = 0.95, and h = 0.15 in each method, but with the same initial conditions as in Case 2 and the same parameters as in Figure 1. Where for each method, we computed the Relative Approximate Error (RAE) which is defined as below to show that our presented M. M. Khader / Eur. J. Pure Appl. Math, 17 (4) (2024), 2676-2691 2686 Fig. 4: The solution Si(t), i = 1, 2, 3, 4, 5 by the present and GEM methods α = 1 with large initial values. method is more accurate and computationally effective in solving the given system: RAE of Sp(tk) = ∣∣∣∣Sp(tk+1)− Sp(tk) Sp(tk+1) ∣∣∣∣ , p = 1, 2, 3, 4, 5. Also, we computed the allowed time t̄ for obtaining these results by applying the two methods (t̄=75 sec. for MFEM, t̄=160 sec. for GEM). M. M. Khader / Eur. J. Pure Appl. Math, 17 (4) (2024), 2676-2691 2687 Fig. 5: The effect of µ on the approximate solution Si(t), i = 1, 2, 3, 4, 5. M. M. Khader / Eur. J. Pure Appl. Math, 17 (4) (2024), 2676-2691 2688 Fig. 6: The effect of κ on the approximate solution Si(t), i = 1, 2, 3, 4, 5. M. M. Khader / Eur. J. Pure Appl. Math, 17 (4) (2024), 2676-2691 2689 Table 2: Values of the RAE for the present method (first row) and GEM (second row) with α = 0.95. k RAE − S1(tk) RAE − S2(tk) RAE − S3(tk) RAE − S4(tk) RAE − S5(tk) 10 0.2583× 10−5 5.2660× 10−4 8.2511× 10−3 8.6542× 10−3 6.3508× 10−4 7.5621× 10−3 6.3255× 10−2 7.5621× 10−1 0.0226× 10−2 9.2145× 10−3 20 4.2154× 10−4 1.9512× 10−5 7.2584× 10−5 6.2581× 10−5 1.2005× 10−4 8.5201× 10−2 8.2509× 10−2 7.2156× 10−4 9.2812× 10−3 9.2580× 10−2 30 8.0056× 10−6 9.2258× 10−5 2.0897× 10−5 9.2581× 10−6 6.2581× 10−5 8.2581× 10−4 9.2581× 10−3 8.2554× 10−3 6.0057× 10−4 6.0257× 10−4 40 6.2014× 10−5 9.8025× 10−6 5.6241× 10−5 8.2943× 10−5 6.2581× 10−5 5.2210× 10−3 4.0579× 10−4 9.3025× 10−2 8.2581× 10−3 9.2508× 10−3 50 8.2030× 10−7 6.0258× 10−7 5.2017× 10−6 8.2589× 10−6 7.5623× 10−6 8.2587× 10−5 4.2587× 10−5 9.3652× 10−4 8.5214× 10−5 9.0258× 10−4 60 9.2584× 10−6 5.0254× 10−7 8.2541× 10−7 6.2587× 10−7 6.0258× 10−6 7.2564× 10−4 7.5824× 10−5 9.3587× 10−5 9.2587× 10−4 5.3698× 10−4 70 6.2874× 10−8 9.2580× 10−7 7.2054× 10−7 6.3258× 10−8 7.6250× 10−8 4.6528× 10−7 7.9542× 10−5 0.1593× 10−5 5.2810× 10−7 6.2591× 10−6 80 5.2914× 10−7 7.5630× 10−8 4.6528× 10−7 2.2147× 10−8 8.9254× 10−7 6.2581× 10−6 7.5829× 10−7 5.2810× 10−6 5.9190× 10−7 6.2581× 10−5 90 6.2587× 10−9 0.2580× 10−8 1.2250× 10−9 2.5879× 10−9 6.2415× 10−8 5.2614× 10−7 2.2581× 10−6 3.9650× 10−7 4.2589× 10−7 3.2590× 10−6 100 5.6280× 10−8 7.2941× 10−9 1.0236× 10−9 1.5587× 10−8 9.0054× 10−8 8.2541× 10−6 5.2913× 10−7 7.2564× 10−7 6.2591× 10−6 2.9514× 10−7 7. Conclusions and discussions In this research work, we applied the modified fractional Euler method to provide a numerical study and simulation of the smoking model in its fractional form, which is rep- resented by a system of mathematical equations. Therefore, we paid attention to studying the behavior of the numerical solution of the problem under study for different fractional motions as well as for the standard motion α = 1. The resulting numerical solution was also compared with its counterpart and the result of the generalized FEM, and this com- parison helped us to fully confirm that the numerical solution resulting from this applied numerical method agrees well with the standard case of the same method. In addition, through this study, we found that this proposed approach can solve problems effectively with high accuracy. Finally, the present study may contribute to providing more robust physical explanations for future theoretical and computational studies on the same topic. All calculations in this research work were done using Mathematica 8.0. Conflict of interest: The author has no relevant financial or non-financial interests to disclose. Data Availability: The datasets generated during the current study are available from the corresponding author upon reasonable request. REFERENCES 2690 References [1] M. Adel, M.M. Khader, I. Alraddadi, N.H. Sweilam, M.B. Riaz, and H. Ahmad. Ef- fective numerical method for studying the fractal-fractional smoking model. Fractals, pages xx–yy, 2024. [2] M. Adel, H.M. Srivastava, and M.M. Khader. Implementation of an accurate method for the analysis and simulation of electrical R-L circuits. Math. Meth. Appl. Sci., 46:8362–8371, 2023. [3] Z. Alkhudhari, S. Al-Sheikh, and S. Al-Tuwairqi. Global dynamics of a mathematical model on smoking. Appl. Math., pages 1–16, 2014. [4] Y.N. Anjam, R. Shafqat, I.E. Sarris, M. Rahma, S. Touseef, and M. Arshad. A fractional order investigation of smoking model using Caputo-Fabrizio differential operator. Fractal Fract., 6:1–15, 2022. [5] A.A.M. Arafa, S.Z. Rida, and M. Khalil. Fractional modeling dynamics of HIV and CD4+ T-cells during primary infection. Nonl. Bio. Phy., 6:1–7, 2012. [6] I.M. Batiha, A. Bataihah, Abeer A. Al-Nana, S. Alshorm, I.H. Jebril, and A. Zraiqat. A numerical scheme for dealing with fractional initial value problem. Int. J. Innov. Comput. Inform. Cont., 19:763–774, 2023. [7] C.C. Garsow, G.J. Salivia, and A.R. Herrera. Mathematical models for dynamics of tobacco use, recovery and relapse. Technical report, Cornell University, London, UK, 2000. [8] O.K. Ham. Stages and processes of smoking cessation among adolescents. West. J. Nurs. Res., 29:301–315, 2007. [9] M.A. Hammad, I.H. Jebril, S. Alshorm, I.M. Batiha, and N.A. Hammad. Numerical solution for fractional-order mathematical model of immune-chemotherapeutic treat- ment for breast cancer using modified fractional formula. Int. J. Anal. Appl., 21:1–9, 2023. [10] M. Izadi and H.M. Srivastava. Generalized shifted Airfoil polynomials of the second kind to solve a class of singular electrohydrodynamic fluid model of fractional order. Fractal Fract., 7:93–105, 2023. [11] M. Izadi and H.M. Srivastava. An accurate numerical algorithm based on the gener- alized Narayana polynomials to solve a class of Caputo-Fabrizio and Liouville-Caputo fractional-order delay differential equations. Math. Commun., 29:61–81, 2024. [12] M.M. Khader and K.M. Saad. On the numerical evaluation for studying the fractional KdV, KdV-Burger’s, and Burger’s equations. The Eur. Ph. J. Plus, 33:1–13, 2018. REFERENCES 2691 [13] T. Liu, T.J. Deiss, M.W. Lippi, A. Jauregui, K. Vessel, S. Ke, A. Belzer, H. Zhuo, K.N. Kangelaris, and et al. A.D. Gomez. Alternative tobacco product use in critically ill patients. Int. J. Environ. Res. Public Health, 17:1–15, 2020. [14] World Health Organisation. Monitoring tobacco use and prevention policies. executive summary, 2017. [15] K.M. Owolabi, A. Atangana, and A. Akgul. Modelling and analysis of fractal- fractional partial differential equations: Application to reaction-diffusion model. Alex. Eng. J., 59:2477–2490, 2020. [16] M.P. Yadav and R. Agarwal. Numerical investigation of fractional-fractal Boussinesq equation. Chaos, 29:01–31, 2019. [17] A. Yildirim and Y. Cherruault. Analytical approximate solution of a SIR epidemic model with constant vaccination strategy by homotopy perturbation method. Emer- ald, 38:1566–1581, 2009. [18] M.O. Zaid and N.T. Shawagfeh. Generalized Taylor’s formula. App. Math. Com., 186:286–293, 2007.