Analytical approximate solution of fractional order smoking epidemic model

Abstract

In this paper, the fractional smoking epidemic model is presented. The model is presented in terms of Caputo’s fractional derivation. The fractional differential transformation method (FDTM) is presented to find an approximate analytical solution to the model. The method is tested on the model and the solution is compared with the homotopy transform method. The method shows the form of fast converging series and the results prove the applicability of the proposed technique, which gives accurate results.

Introduction

Smoking is considered one of the biggest health problems in the world. The fact that more than 5 million deaths in the world are due to smoking makes it a priority for health organizations to pay more attention to this issue. The risk of suffering a heart attack is 70% higher in smokers than in non-smokers, and smokers are 10% more likely to develop lung cancer than others. This is because smoking can affect various organs, causing them to break down in the short and long term. In the short term, high blood pressure, bad breath, cough, and diseased teeth are some of the consequences of smoking, while various cancers are the long-term consequences of smoking. According to the report by WHO, smokers live 10–13 years shorter than non-smokers. In addition, the World Health Organization has counted that more than 250 million children have died as a result of smoking and more than 10 million people will perish from smoking-related diseases by 2030. ¹ For this reason, smoking can be considered a leading global health problem that requires further action to make people aware of the dangers and effects of smoking. Mathematics can play an important role in this process by providing different models of smoking that examine the effects of smoking on different individuals. As recently as 1997, Castillo-Garsow et al. ² proposed a numerical model to simulate the smoking cessation model. They approached this problem by considering three classes of smokers: smokers ( $S$ ), potential smokers ( $P$ ), and quitters ( $Q$ ). They used the Routh-Huwartiz method to simulate the local stability of the model, while the global stability was only addressed by some numerical simulations. In recent years, many other researchers have proposed other methods to address the smoking phenomenon and better understand its effect and dynamics. For example, Huo and Zhu studied the effects of relapse in smokers who quit smoking using a new mathematical model. ³ Lahrouz et al. ⁴ studied the stochastic stability of these models using the Lyapunov function. Guerrero et al. ⁵ presented an approximate analytical solution using the homotopy analysis method for a dynamic model for smoking habits in Spain with different mortality rates. Singh et al. ⁶ discussed a new fractional order dynamic smoking cessation model with a non-singular kernel and used numerical based iterative technique to simulate it. Sikander et al. ⁷ developed a parameter variation combined with an auxiliary parameter to find an approximate solution of the epidemic model for the evolution of smoking habits in a constant population. Zaman ⁸ studied the qualitative response to smoking cessation dynamics.

Delgado et al. ⁹ also discussed the dynamics of the model using the homotopy perturbation method along with the Laplace transform method to obtain analytical solutions and studied the uniqueness and existence of the solution using the fixed point theorem and the Picard-Lindelof approach. Some studies on these models can be found as follows: Erturk et al. ¹⁰ applied two numerical methods for the fractional order smoking model, and Zeb et al. ¹¹ also applied two methods for fractional order smoking model with fallback class. Alkhudhari et al. ¹² performed a stability analysis of the smoking model. Zeb et al. ¹³ performed a dynamic analysis of the smoker model. For more examples of this model, we refer to Zeb et al., ¹⁴ where the proposed smoker model is successfully solved using two different numerical methods.

Veeresha et al. ¹⁵ have applied the q-homotopy analysis transformation method for the fractional order smoking epidemic model.

The fractional analysis is very useful in simulating phenomena in physics, mathematics, biology, fluid mechanics, chemistry, signal processing, and many other disciplines. Fractal derivatives and integrals are important aspects in all branches of fractional calculus. The basic definitions for fractional order derivatives were introduced by Caputo, ¹⁶ which allows for conventional initial and boundary conditions in real problems. After Caputo’s definition of fractional order derivative, there were a number of other definitions that can also describe fractional order, including Caputo’s definition of fractional derivative. The fractional derivative of Caputo was first derived by Khalil et al. ¹⁷ Numerous researchers use this definition to solve real-world problems. For example, Eslami and Rezazadeh ¹⁸ use the first integral method together with Caputo’s fractional derivative definition to solve the Wu-Zhang system. Moreover, Ekici et al. ¹⁹ adapted the same method for solving the optical soliton perturbation with fractional-temporal evolution. Aminikhah et al. ²⁰ use the sub-equation method to solve the regularized long wave equations with fractional derivatives of Caputo. El-Ajou et al. ²¹ employed the same definition for solving solitary solutions for time-fractional nonlinear dispersive PDEs. Singh et al.^22,23 gave a brief explanation of the fractional population model and the fractional coupled Burger equation, respectively. Area et al. ²⁴ investigated a classical and fractional model that simulates the Ebola epidemic model. Carvalho and Pinto ²⁵ examined a fractional-order model of the co-infection of Malaria and HIV/AIDS with delay and also studied the stability of the disease-free equilibrium. Further details regarding the fractional models and the methods for solving them can be found in Refs.^26–30 and references therein. Recently, researchers have designed several smoking models under various.^11,31–35

The differential transform method (DTM) was first created by Zhou to obtain approximate-analytical solutions to ordinary differential equations using Taylor series formulation. ³⁶

Then, Arikoglu and Ozkol ³⁷ developed this approach by using the power series formulation for fractional-order differential equations (see Günerhan et al. ³⁸ for more background information about the applicability of this technique for solving the HIV model). These include techniques described in Ullah et al., ³⁹ Alkresheh and Ismail, ⁴⁰ Patel and Tandel, ⁴¹ and Chiranjeevi and Biswas. ⁴²

DTM has been applied for solving fractional differential-algebraic equations in Birol et al. ⁴³ In addition, a modified version of the method has been used for solving a multi-term time-fractional diffusion equation in Abuasad et al. ⁴⁴ Nazari and Shahmorad ⁴⁵ used the method for solving a fractional order integro differential equation acquiring good results.

In this paper, we are concerned with analyzing the epidemic system of fractional order for modeling the smoking dynamics. The fractional derivative is considered in the Caputo fractional form and the system is considered in the form.

$$\begin{array}{cc}\multicolumn{1}{c}{D_{*}^{\alpha }P}& =\lambda -\beta \mathrm{PS}-\eta P,\hfill \\ \multicolumn{1}{c}{D_{*}^{\alpha }O}& =\beta \mathrm{PS}-{\alpha }_{1}O-\eta O,\hfill \\ \multicolumn{1}{c}{D_{*}^{\alpha }S}& ={\alpha }_{1}O+{\alpha }_2\mathrm{SQ}-(\eta +\gamma )S,\hfill \\ \multicolumn{1}{c}{D_{*}^{\alpha }Q}& =-{\alpha }_2\mathrm{SQ}-\eta Q+\gamma (1-\delta )S,\hfill \\ \multicolumn{1}{c}{D_{*}^{\alpha }L}& =\gamma S\delta -\eta L,\hfill \end{array}$$

(1)

with initial conditions

$$P\left(0\right)=P_0,O\left(0\right)=O_0,S\left(0\right)=S_0,Q\left(0\right)=Q_0,L\left(0\right)=L_0,$$

(2)

$D_{*}^{\alpha }$ is the operator of the fractional-order derivative of order $\alpha ,(0<\alpha \le 1)$ in the system of equation (1). Let the total population size at time $t$ is represented by $N(t)$ . We separate the population $N(t)$ into five subgroups, namely potential smoker $P(t)$ , smoker $S(t)$ , occasional smoker $O(t)$ , permanently quit smoker $L(t)$ , and temporarily quit smoker $Q(t)$ . Each parameter in (1) represents the following ¹⁵ :

$\lambda$ represents recruitment rate; $\beta$ represents linkage index between potential and light smokers, $\eta$ represents natural mortality rate, $\gamma$ represents smoking cessation rate, $\delta$ symbolizes smokers who quit smoking forever, ${\alpha }_1$ represents rate between light and heavy smokers, and ${\alpha }_2$ represents linkage index between temporary quitters who may start smoking again.

The development of the paper is as follows: In Section 2, a brief overview of the Caputo fractional order derivation is given. Section 3 discusses permilnarios for the fractional differential transform method, which will be needed later. In Section 4, the proposed mathematical model of smoking is solved using FDTM and a comparison is made with the q-homotopy analysis transform method (q-HATM) from the literature. Section 4 is the last section where a conclusion for the study is given.

Basic definitions

In this section, we will introduce some basic definitions and properties of the theory of fractional calculus that will be later.

• Definition 1. A real function $f\left(x\right),x>0$ is said to be in the space $C_{\mu },\mu ϵR$ if there exists a real number $P>\mu$ such that $f\left(x\right)=x^p{f}_1\left(x\right)$ where $f_1\left(x\right)ϵC[0,\infty )$ . Clearly $C_{\mu }<C_{\beta }$ if $\mu <\beta$ .

• Definition 2. A function $f\left(x\right),x>0$ is said to be in the space $C_{\mu }^m$ , $mϵN\cup \left\{0\right\}$ if $f^{\left(m\right)}\in C_{\mu }$ .

• Definition 3. The Riemann-Liouville fractional integral operator of the order $\alpha >0$ of a function, $f\in C_{\mu },\mu \ge -1$ is defined as

$$\begin{array}{c}(J_a^{\alpha }f)(x)=\frac{1}{\mathrm{\Gamma }(\alpha )}\underset{a}{\overset{x}{\int }}(x-\tau {)}^{\alpha -1}f(\tau )d\tau ,x>a,\hfill \\ (J_a^{0}f)(x)=f(x).\hfill \end{array}$$

(3)

All the properties of the operator $J^{\alpha }$ can be found in Caputo ⁴⁶ which we mention only the following

For $f\in C_{\mu },\mu \ge -1,\alpha ,\beta \ge 0,$ and $\gamma >-1$

$$(J_a^{\alpha }{J}_a^{\beta }f)(x)=(J_a^{\alpha +\beta }f)(x),$$

(4)

$$\left(J_a^{\alpha }{J}_a^{\beta }f\right)\left(x\right)=\left(J_a^{\beta }{J}_a^{\alpha }f\right)\left(x\right),$$

(5)

$$J_a^{\alpha }{x}^{\gamma }=\frac{\mathrm{\Gamma }(\gamma +1)}{\mathrm{\Gamma }(\alpha +\gamma +1)}{x}^{\alpha +\gamma }.$$

(6)

Equilibrium point and stability of the fractional order smoking model:

To calculate endemic equilibrium points, we will consider infected nodes.

The Jacobian of system (1) is as

$$J=\left[\begin{array}{ccccc}-\beta S-\eta \hfill & 0\hfill & -\beta P\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{\beta S}& -{\alpha }_1-\eta \hfill & \beta P\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& {\alpha }_1\hfill & {\alpha }_{2}Q-(\eta +\gamma )\hfill & {\alpha }_{2}S\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & {\alpha }_{2}Q+\gamma (1-\delta )\hfill & -{\alpha }_{2}S-\eta \hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & \delta \gamma \hfill & 0\hfill & -\eta \hfill \end{array}\right]$$

The equilibrium points are given as $E_0=(0,0,0,0,0)$ . The Jacobean matrix at equilibrium point (0, 0, 0, 0, 0) is given as follows

$$J(E_0)=\left[\begin{array}{ccccc}-\eta \hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& -{\alpha }_1-\eta \hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& {\alpha }_1\hfill & -(\eta +\gamma )\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & \gamma (1-\delta )\hfill & -\eta \hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & \delta \gamma \hfill & 0\hfill & -\eta \hfill \end{array}\right]$$

Now by substituting the values of parameters represented in Table 1 into $J(E_0),$ we have

$$J(E_0)=\left[\begin{array}{ccccc}-0.05\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& -0.052\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0.002\hfill & -0.85\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0.72\hfill & -0.05\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0.08\hfill & 0\hfill & -0.05\hfill \end{array}\right]$$

Table 1.

Parameters values defined in equation (1).

Parameter	Value
$\lambda$	1.0
$\beta$	0.14
$\eta$	0.05
${\alpha }_1$	0.002
${\alpha }_2$	0.0025
$\gamma$	0.8
$\delta$	0.1

The Eigen values corresponding to matrix $J(E_0)$ are ${\lambda }_1=-0.85,{\lambda }_2=-0.052,{\lambda }_3=-0.05,$ ${\lambda }_4=-0.05\mathrm{and}{\lambda }_5=-0.05.$

Since $J$ has negative eigenvalues, the system (1) is stable at the equilibrium point (0, 0, 0, 0, 0) for 0 < ${\alpha }_i$ < 1, i = 1, 2, 3, 4, 5.

In system (1) that at least one of the infected components is non-zero for find the endemic equilibria, we consider (see, for details, previous study ⁴⁷ )

$E^{*}=(S^{*},P^{*},O^{*},Q^{*},L^{*})$ represents endemic equilibrium of the system (1). By solving system (1) in a steady state, we have

$$\begin{array}{c}{P}^{*}=\frac{\lambda }{\beta S^{*}+\eta },O^{*}=\frac{\lambda \beta S^{*}}{(\beta S^{*}+\eta )({\alpha }_1+\eta )},\hfill \\ Q^{*}=\frac{\gamma (1-\delta )S^{*}}{({\alpha }_2{S}^{*}+\eta )},L^{*}=\frac{\delta \gamma S^{*}}{\eta }.\hfill \end{array}$$

Fractional differential transform method

The FDTM is semi-numerical and analytical approach, and it is considered as a traditional DTM reform. The differential transformation, $s(x)$ , can be expressed as follows:

$$S_{\phi }\left(k\right)=\frac{1}{\mathrm{\Gamma }(\phi k+1)}{\left[{\left(D_{x_0}^{\phi }\right)}^{k}s\left(x\right)\right]}_{x=x_0},$$

(7)

The differential transform (DT) inverse of $S_{\phi }\left(k\right)$ is provided as

$$s\left(x\right)=\sum _{k=0}^{\infty }{S}_{\phi }\left(K\right){(x-x_0)}^{\phi k}\equiv D^{-1}{S}_{\phi }\left(k\right),$$

(8)

By substituting equation (7) into equation (8), the following is obtained:

$$\sum _{k=0}^{\infty }\frac{1}{\mathrm{\Gamma }(\phi k+1)}{\left[{\left(D_{x_0}^{\phi }\right)}^{k}s\left(x\right)\right]}_{x=x_0}{(x-x_0)}^{\phi k}=s(x),$$

(9)

Assume that $S_{\phi }(k)$ is the DT of $s(x)$ . The approximate function $S\left(k\right)$ is written as:

$$s\left(x\right)=\sum _{k=0}^n{S}_{\phi }\left(K\right){(x-x_0)}^{\phi k}.s\left(x\right)=\sum _{k=0}^n{S}_{\phi }\left(K\right){(x-x_0)}^{\phi k}.$$

(10)

The above equation (equation (10)) represents the differential transformation resulting from the application of the Taylor series expansion, although the symbolic evaluation of the derivatives is not applicable to this technique. Moreover, comparative derivatives were obtained by applying the iterative procedure. The original function is represented by a lower case letter, while the transformed function is represented by an upper case letter. According to equations (9) and (10), it is obvious that the transformed functions have the basic mathematical values given in Table 2.

Table 2.

The operations for the fractional differential transform method.

Unique function	Transformed function
$w(x)=d(x)\pm b(x)$	$W(k)=D_{\phi }(k)\pm B_{\phi }(k)$
$w(x)=\mathrm{cb}(x)$	$W(k)=c{B}_{\phi }(k)$
$w(x)=q(\mathrm{cx})$	$W(k)=c^k{Q}_{\phi }(k)$
$w(x)=q(\frac{x}{c})$	$W(k)=\frac{Q_{\phi }(k)}{c^k}$
$w(x)=D_{x_0}^{m\phi }q(x)$	$W(k)=\frac{\mathrm{\Gamma }(\phi k+m\phi +1)}{\mathrm{\Gamma }(\phi k+1)}{Q}_{\phi }(k+m)$
$w(x)=x^m$	$W\left(k\right)=\delta (k-m)=\{\begin{array}{c}1\hfill \\ \multicolumn{1}{c}{0}\end{array}\begin{array}{c}k=m\hfill \\ \multicolumn{1}{c}{k\ne m}\end{array}$
$w(x)=e^{x+c}$	$W(k)=\frac{e^c}{k!}$

In this part, we will show converges of the differential transform method to the exact solution by proof of some theorems. ⁴⁷

• Theorem 1: Consider the first order ordinary differential equation in here $m,n$ are constants

$$m\frac{\mathrm{dw}(x)}{\mathrm{dx}}+\mathrm{nw}(x)=0,w\left(0\right)=\gamma ,$$

(a)

The differential transform method converges to the exact solution.

Proof. The general solution of differential equation is as $w\left(x\right)=\beta e^{-\mathrm{nx}/m}$ .

Now By taking the differential transform of equation (a), we have

$$mW[K+1]+\mathrm{nW}\left[K\right]=0,W\left[0\right]=\beta ,$$

We have $W[K+1]=\frac{-\mathrm{nW}\left[K\right]}{m}.$ By solve this recurrence relation, we get $W\left[K\right]=\frac{1}{k!}{\left(\frac{-n}{m}\right)}^k\beta$

By using equation (8), the solution obtained by DTM is as follows:

$$w(x)=\beta \sum _{k=0}^{\infty }\frac{1}{k!}{\left(\frac{-nx}{m}\right)}^k=\beta e^{-nx/m}.$$

• Theorem 2: Consider the homogeneous ordinary differential equation of order n, where $v_i,i=1,2,\dots ,n$ are variables

$$v_0\left(x\right)w+\sum _{r=1}^n{v}_r\left(x\right)\frac{d^{r}w}{d{x}^r}=0,$$

(b)

With initial conditions

$$w\left(x\right)={\beta }_0,{\frac{d^{i}w}{d{x}^i}|}_{x=0}={\beta }_i,i=1,2,\dots ,n-1,$$

(c)

The differential transform method converges to the exact solution.

Proof. By using with both Table 1 and the equation (b), we have

$$\sum _{m=0}^k(V_0[k]W[k-m]+\sum _{r=1}^n[V_r\left[k\right][\underset{j=1}{\overset{r}{\Pi }}(k-m+j)]W[k-m+r]])=0,$$

From the initial conditions, we have

$$W\left[0\right]={\beta }_0,W\left[i\right]=\frac{{\beta }_i}{i!},i=1,2,\dots ,n-1,$$

By substituting the differential transforms, we have

$$w\left(x\right)=\sum _{k=0}^{\infty }W\left[k\right]x^k,v_r\left(x\right)=\sum _{k=0}^{\infty }{A}_r\left[k\right]x^k,$$

Substituting the above formulas in equation (b), we get

$$\begin{array}{c}\multicolumn{1}{c}{0=\sum _{k=0}^{\infty }{A}_0\left[k\right]x^k\sum _{k=0}^{\infty }W\left[k\right]x^k+\sum _{r=1}^n(\sum _{k=0}^{\infty }{A}_r\left[k\right]x^k.\left(\frac{d^r}{d{x}^r}\sum _{k=0}^{\infty }W\left[k\right]x^k\right))}\\ =\sum _{k=0}^{\infty }\sum _{m=0}^{\infty }{A}_0\left[m\right]W[k-m]x^k+\sum _{r=1}^n(\sum _{k=0}^{\infty }{A}_r\left[k\right]x^k.\left(\sum _{m=0}^{\infty }W[m+r]\right)\left(\underset{j=0}{\overset{r-1}{\Pi }}(m+r-j)\right)x^m)\hfill \\ =\sum _{k=0}^{\infty }\sum _{m=0}^k(A_0\left[m\right]W[k-m]+\sum _{r=1}^n{A}_r\left[m\right]W[k-m+r](\underset{j=1}{\overset{r}{\Pi }}(k-m+j))x^k\hfill \end{array}$$

we comparing coefficients and get

$$\begin{array}{c}\sum _{m=0}^k{A}_0\left[m\right]W[k-m]+\sum _{r=1}^n{A}_r\left[m\right]W[k-m+r]\hfill \\ (\underset{j=1}{\overset{r}{\Pi }}(k-m+j)=0.\hfill \end{array}$$

Next, we will use the last-mentioned definitions for solving the model described in equation (1) through a numerical example to prove the efficiency of the proposed technique.

Numerical example

In this section, we will examine our method based on the FDTM for the system of fractional order smoking epidemic model represented in equation (1) for various values of the parameters.

By applying FDTM to the proposed smoking mathematical model, and by using with both Table 1 and the equation (1), a system of equations is obtained as follows:

$$\begin{array}{c}{P}_{\alpha }(k+1)=\frac{\mathrm{\Gamma }(1+\alpha k)}{\mathrm{\Gamma }(\alpha (k+1)+1)}\hfill \\ [\lambda \delta (k)-\beta \sum _{l=0}^k{P}_{\alpha }\left(l\right)S_{\alpha }(k-l)-\eta P_{\alpha }(k)],\hfill \\ O_{\alpha }(k+1)=\frac{\mathrm{\Gamma }(1+\alpha k)}{\mathrm{\Gamma }(\alpha (k+1)+1)}\hfill \\ [\beta \sum _{l=0}^k{P}_{\alpha }\left(l\right)S_{\alpha }(k-l)-\alpha O_{\alpha }(k)-\eta O_{\alpha }(k)],\hfill \\ S_{\alpha }(k+1)=\frac{\mathrm{\Gamma }(1+\alpha k)}{\mathrm{\Gamma }(\alpha (k+1)+1)}\hfill \\ [{\alpha }_1{O}_{\alpha }(k)+{\alpha }_2\sum _{l=0}^k{Q}_{\alpha }\left(l\right)S_{\alpha }(k-l)-(\eta +\gamma )S_{\alpha }(k)],\hfill \\ Q_{\alpha }(k+1)=\frac{\mathrm{\Gamma }(1+\alpha k)}{\mathrm{\Gamma }(\alpha (k+1)+1)}\hfill \\ [-{\alpha }_2\sum _{l=0}^k{Q}_{\alpha }\left(l\right)S_{\alpha }(k-l)-\eta Q_{\alpha }\left(l\right)+\gamma (1-\delta )S_{\alpha }\left(k\right)]\hfill \\ L_{\alpha }(k+1)=\frac{\mathrm{\Gamma }(1+\alpha k)}{\mathrm{\Gamma }(\alpha (k+1)+1)}[\gamma \delta S_{\alpha }\left(k\right)-\eta L_{\alpha }\left(k\right)].\hfill \end{array}$$

(11)

From the initial condition given by equation (2), we obtain:

$$\begin{array}{c}{P}_{\alpha }\left(0\right)=40,O_{\alpha }\left(0\right)=10,S_{\alpha }\left(0\right)=20,\hfill \\ Q_{\alpha }\left(0\right)=10,L_{\alpha }\left(0\right)=5,\hfill \end{array}$$

(12)

Now, by substituting the values of parameters represented in Table 1 into equation (11), we get

$$\begin{array}{c}{P}_{\alpha }(k+1)=\frac{\mathrm{\Gamma }(1+\alpha k)}{\mathrm{\Gamma }(\alpha (k+1)+1)}\hfill \\ [\delta (k)-0.14\sum _{l=0}^k{P}_{\alpha }\left(l\right)S_{\alpha }(k-l)-0.05P_{\alpha }(k)],\hfill \\ O_{\alpha }(k+1)=\frac{\mathrm{\Gamma }(1+\alpha k)}{\mathrm{\Gamma }(\alpha (k+1)+1)}\hfill \\ [0.14\sum _{l=0}^k{P}_{\alpha }\left(l\right)S_{\alpha }(k-l)-0.002O_{\alpha }(k)-0.05O_{\alpha }(k)],\hfill \\ S_{\alpha }(k+1)=\frac{\mathrm{\Gamma }(1+\alpha k)}{\mathrm{\Gamma }(\alpha (k+1)+1)}\hfill \\ [0.002O_{\alpha }(k)+0.0025\sum _{l=0}^k{Q}_{\alpha }\left(l\right)S_{\alpha }(k-l)-0.85S_{\alpha }(k)],\hfill \\ Q_{\alpha }(k+1)=\frac{\mathrm{\Gamma }(1+\alpha k)}{\mathrm{\Gamma }(\alpha (k+1)+1)}\hfill \\ [-0.0025\sum _{l=0}^k{Q}_{\alpha }\left(l\right)S_{\alpha }(k-l)-0.05Q_{\alpha }(k)+0.72S_{\alpha }(k)],\hfill \\ L_{\alpha }(k+1)=\frac{\mathrm{\Gamma }(1+\alpha k)}{\mathrm{\Gamma }(\alpha (k+1)+1)}[0.08S_{\alpha }\left(k\right)-0.05L_{\alpha }\left(k\right)].\hfill \end{array}$$

(13)

where $P\left(t\right),O\left(t\right),S\left(t\right),Q(t)$ and $L(t)$ can be expanded with the aid of equation (8) in the form

$$\begin{array}{c}\multicolumn{1}{c}{P\left(t\right)=\sum _{k=0}^{\infty }{P}_{\alpha }\left(\mathrm{k}\right)t^{k\alpha }=P_{\alpha }\left(0\right)+P_{\alpha }\left(1\right)t^{\alpha }}\\ +P_{\alpha }\left(2\right)t^{2\alpha }+P_{\alpha }\left(3\right)t^{3\alpha }+\dots \hfill \end{array}$$

(14)

$$\begin{array}{c}\multicolumn{1}{c}{O\left(t\right)=\sum _{k=0}^{\infty }{O}_{\alpha }\left(\mathrm{k}\right)t^{k\alpha }=O_{\alpha }\left(0\right)+O_{\alpha }\left(1\right)t^{\alpha }}\\ +O_{\alpha }\left(2\right)t^{2\alpha }+O_{\alpha }\left(3\right)t^{3\alpha }+\dots \hfill \end{array}$$

(15)

$$\begin{array}{c}\multicolumn{1}{c}{S\left(t\right)=\sum _{k=0}^{\infty }{S}_{\alpha }\left(\mathrm{k}\right)t^{k\alpha }=S_{\alpha }\left(0\right)+S_{\alpha }\left(1\right)t^{\alpha }+S_{\alpha }\left(2\right)t^{2\alpha }}\\ +S_{\alpha }\left(3\right)t^{3\alpha }+\dots \hfill \end{array}$$

(16)

$$\begin{array}{c}\multicolumn{1}{c}{Q\left(t\right)=\sum _{k=0}^{\infty }{Q}_{\alpha }\left(\mathrm{k}\right)t^{k\alpha }=Q_{\alpha }\left(0\right)+Q_{\alpha }\left(1\right)t^{\alpha }+Q_{\alpha }\left(2\right)t^{2\alpha }}\\ +Q_{\alpha }\left(3\right)t^{3\alpha }+\dots \hfill \end{array}$$

(17)

$$\begin{array}{c}\multicolumn{1}{c}{L\left(t\right)=\sum _{k=0}^{\infty }{L}_{\alpha }\left(\mathrm{k}\right)t^{k\alpha }=L_{\alpha }\left(0\right)+L_{\alpha }\left(1\right)t^{\alpha }+L_{\alpha }\left(2\right)t^{2\alpha }}\\ +L_{\alpha }\left(3\right)t^{3\alpha }+\dots \hfill \end{array}$$

(18)

Numerical simulation of $P\left(t\right),O\left(t\right),S\left(t\right),Q(t)$ and $L(t)$ are shown in Tables 3to 11 for different values of $\alpha$ . We compared the result of FDTM with the results of the modified homotopy analysis transformation method (q-HATM) ¹⁵ with the same values of the parameters listed in Table 1. In the α = 1 case, we compared the result of FDTM with the results of the modified homotopy analysis transformation method (q-HATM), and compared both methods with the ODE45(Runge-Kutta (4,5) method). The data for both cases are compared with the known Matlab code (ODE45) for numerical simulation of this system and it is found that the figures are in perfect agreement with the results of the proposed technique for α = 1. This ensures that our technique accurately reproduces the results for α = 0.8, 0.9. Figures 1 to 4 also shows the approximate solution profile for the variables and for different values of $\alpha$ . The results for the solution of model (1) are examined for different values of α to prove the effectiveness and validity of the proposed algorithm. The values of the parameters used in the numerical simulations are summarized in Table 2.

Table 3.

Approximate solution of $P\left(t\right)$ and $S(t)$ for $\alpha =1$ and the obtained results of FDTM with and q-HATM.

	$P\left(t\right)$			$S(t)$
t	ODE 45	FDTM	q-HATM	ODE 45	FDTM	q-HATM
0	40	40	40	20	20	20
0.1	30.5072	30.9766	30.7667	18.4221	18.5009	18.4244
0.2	23.7850	24.1740	25.6668	16.9756	17.0509	16.9938
0.3	18.9262	19.2495	24.7002	15.6484	15.7202	15.7080
0.4	15.3484	15.6184	27.8670	14.4298	14.4982	14.5671
0.5	12.6687	12.8957	35.1672	13.3102	13.3753	13.5711
0.6	10.6306	10.8229	46.6008	12.2810	12.3429	12.7201
0.7	9.05875	9.2230	62.1678	11.3344	11.3932	12.0138
0.8	7.8313	7.9727	81.8682	10.4633	10.5192	11.4525
0.9	6.8620	6.9847	105.702	9.6614	9.7144	11.0360
1.0	6.0889	6.1962	133.669	8.9230	8.9733	10.7645

Table 4.

Approximate solution of $Q\left(t\right)$ and $L(t)$ for $\alpha =1$ and the obtained results of FDTM with and q-HATM.

	$Q\left(t\right)$			$L(t)$
t	ODE 45	FDTM	q-HATM	ODE 45	FDTM	q-HATM
0	10	10	10	5	5	5
0.1	11.2780	11.2142	11.2760	5.1282	5.1218	5.1280
0.2	12.4398	12.3791	12.4241	5.2438	5.2377	5.2422
0.3	13.4959	13.4382	13.4443	5.3477	5.3420	5.3426
0.4	14.4556	14.401	14.3365	5.4410	5.4356	5.4291
0.5	15.3276	15.2759	15.1008	5.5245	5.5194	5.5017
0.6	16.1194	16.0705	15.7372	5.5990	5.5942	5.5605
0.7	16.8381	16.7919	16.2456	5.6652	5.6608	5.6054
0.8	17.4898	17.4462	16.6260	5.7239	5.7197	5.6365
0.9	18.0802	18.0392	16.8786	5.7756	5.7717	5.6537
1.0	18.6145	18.5759	17.0032	5.8209	5.8173	5.6570

Table 5.

Approximate solution of $O\left(t\right)$ for $\alpha =1$ and the obtained results of FDTM with and q-HATM.

t	ODE 45	FDTM	q-HATM
0	10	10	10
0.05	15.0909	14.5834	15.0559
0.1	19.3402	18.8789	19.0756
0.15	22.9048	22.4858	22.0590
0.20	25.9090	25.5286	24.0063
0.25	28.4522	28.1067	24.9173
0.30	30.6139	30.2999	24.7922
0.35	32.4586	32.1728	23.6308
0.40	34.0380	33.7777	21.4332
0.45	35.3949	35.1575	18.1994
0.50	36.5641	36.3473	13.9294

Table 6.

Approximate solution of $P\left(t\right)$ and $S(t)$ for $\alpha =0.9$ and the obtained results of FDTM with and q-HATM.

	$P\left(t\right)$		$S(t)$
t	FDTM	q-HATM	FDTM	q-HATM
0	40	40	20	20
0.1	26.2807	29.1162	17.5281	17.9798
0.2	17.8878	26.0052	15.2715	16.4515
0.3	12.8287	28.4733	13.3178	15.1913
0.4	9.6387	35.8752	11.6235	14.1491
0.5	7.5495	47.8402	10.1521	13.2994
0.6	6.1373	64.1152	8.8727	12.6260
0.7	5.1580	84.5111	7.7594	12.1178
0.8	4.4654	108.878	6.7897	11.7662
0.9	3.9687	137.095	5.9444	11.5645
1.0	3.6099	169.057	5.2073	11.5073

Table 7.

Approximate solution of $Q\left(t\right)$ and $L(t)$ for $\alpha =0.9$ and the obtained results of FDTM with and q-HATM.

	$Q\left(t\right)$		$L(t)$
t	FDTM	q-HATM	FDTM	q-HATM
0	10	10	5	5
0.1	11.9966	11.6331	5.1998	5.1636
0.2	13.7923	12.8520	5.3765	5.2841
0.3	15.3195	13.8408	5.5235	5.3803
0.4	16.6167	14.6413	5.6446	5.4564
0.5	17.7164	15.2749	5.7434	5.5148
0.6	18.6462	15.7550	5.8229	5.5567
0.7	19.4293	16.0913	5.8856	5.5832
0.8	20.0859	16.2908	5.9338	5.5950
0.9	20.6329	16.3594	5.9693	5.5928
1.0	21.0851	16.3015	5.9939	5.5770

Table 8.

Approximate solution of $O\left(t\right)$ for $\alpha =0.9$ and the obtained results of FDTM with and q-HATM.

t	FDTM	q-HATM
0	10	10
0.05	17.3313	16.6946
0.1	23.4740	20.6741
0.15	28.0928	22.8892
0.20	31.6072	23.5858
0.25	34.3096	22.8977
0.30	36.4067	20.9135
0.35	38.0470	17.6978
0.40	39.3381	13.3004
0.45	40.3594	7.76208
0.50	41.1701	1.11636

Table 9.

Approximate solution of $P\left(t\right)$ and $S(t)$ for $\alpha =0.8$ and the obtained results of FDTM with and q-HATM.

	$P\left(t\right)$		$S(t)$
t	FDTM	q-HATM	FDTM	q-HATM
0	40	40	20	20
0.1	20.3649	28.0338	16.0031	17.4502
0.2	11.6303	28.5365	12.6902	15.8891
0.3	7.5441	35.8119	10.0857	14.7230
0.4	5.4384	48.4484	8.0307	13.8383
0.5	4.2779	65.6923	6.4053	13.1806
0.6	3.6132	87.0555	5.1169	12.7172
0.7	3.2304	112.189	4.0942	12.4259
0.8	3.0192	140.825	3.2815	12.2905
0.9	2.9188	172.753	2.6349	12.2987
1.0	2.8938	207.799	2.1203	12.4407

Table 10.

Approximate solution of $Q\left(t\right)$ and $L(t)$ for $\alpha =0.8$ and the obtained results of FDTM with and q-HATM.

	$Q\left(t\right)$		$L(t)$
t	FDTM	q-HATM	FDTM	q-HATM
0	10	10	5	5
0.1	13.2126	12.0554	5.3198	5.20537
0.2	15.8016	13.2887	5.5687	5.32615
0.3	17.7628	14.1876	5.7471	5.41201
0.4	19.2377	14.8466	5.8702	5.47262
0.5	20.3338	15.3112	5.9498	5.51271
0.6	21.1337	15.6089	5.9953	5.53511
0.7	21.7012	15.7583	6.0139	5.54178
0.8	22.086	15.7730	6.0114	5.53414
0.9	22.3271	15.6636	5.9925	5.51326
1.0	22.4551	15.4383	5.9606	5.48004

Table 11.

Approximate solution of $O\left(t\right)$ for $\alpha =0.8$ and the obtained results of FDTM with and q-HATM.

T	FDTM	q-HATM
0	10	10
0.05	21.3023	18.4931
0.1	29.2192	21.6874
0.15	34.2409	22.3060
0.20	37.5148	20.9520
0.25	39.6940	17.9345
0.30	41.1651	13.4483
0.35	42.1640	7.63003
0.40	42.8397	3.58228
0.45	43.2886	1.61433
0.50	43.5748	0.58228

Figure 1.

The solution of $P\left(t\right),S\left(t\right),L\left(t\right),Q\left(t\right),O\left(t\right)$ obtained by ODE 45( Runge-Kutta (4,5)) for $0<t<1$ .

Figure 2.

The solution of $P\left(t\right),S\left(t\right),L\left(t\right),Q\left(t\right),O\left(t\right)$ obtained by ODE 45, LADM and q-homotopy analysis transform method for $\alpha =1$ .

Figure 3.

The solution of $P\left(t\right),S\left(t\right),L\left(t\right),Q\left(t\right),O\left(t\right)$ obtained by ODE 45, LADM and q-homotopy analysis transform method for $\alpha =0.9$ .

Figure 4.

The solution of $P\left(t\right),S\left(t\right),L\left(t\right),Q\left(t\right),O\left(t\right)$ obtained by LADM and q-homotopy analysis transform method for $\alpha =0.8$ .

It is noticed that the FDTM method is effective in generating approximate solutions for the proposed mathematical model. Numerical simulations of $P\left(t\right),O\left(t\right),S\left(t\right),Q(t)$ and $L(t)$ are shown in Figures 1 to 4 over an interval of $0<t<1$ for various values of $\alpha =0.8,0.9,1$ and all results were compared using the q-homotopy method for the problem under study. The values of the parameters used in the numerical simulations are summarized in Table 2.

Figure 1 shows the result of numerical solution of P(t), S(t), L(t), Q(t), O(t) that obtained by Runge-Kutta (4,5). In Figure 2, P(t) shows the number of potential smokers decreases with time. S(t) shows the number of smokers decreases with time, and O(t) shows the number of occasional smokers increases with time. L(t) shows the number of permanently quit smokers increases with time and Q(t) shows the group of temporarily quit smokers increases with time This shows the accuracy of the proposed model. In this figure, we compare the results of the ODE 45, FDTM, and q-homotopy analysis transformation method in the ALpha = 1 case. We can see that Clearly, we can see that the results of FDTM converge with the results of the Runge-Kutta (4,5) (ODE 45). The results data are ensured to be accurate and investigated for different values of the fractional-order α including α = 1 which is the basic model. In addition, the results are compared with some numerical verification tools such as the Matlab code (ODE45) which ensures that these results are accurate and can describe the true dynamics of these models. The behavior of each variable is presented through figures and the method is shown to be a straightforward, efficient, and simple approach and this may help to adapt to several other related models. Our used proposed method is reliable and efficient. Numerical experiments have been performed with the applied methods for different values of α, which have successfully provided good results for the studied problem.

Conclusion

In this research work, the proposed mathematical model of nonlinear fractional smoking associated with a modified version of Caputo’s fractional derivative was successfully solved using the fractional differential transform method. All obtained results were analyzed and compared for different cases. By using ode45 which is based on an explicit Runge-Kutta (4,5) formula, we got the numerical solution of the model, and by using the Fractional differential transform method we got approximate analytical solutions. We compare both results of FDTM and the modified homotopy analysis transformation method (q-HATM) in the ALpha = 1 case. Clearly, the results of FDTM converge to the exact solution as shown in Tables 3 to 6. The results obtained by the FDTM method agree well with the given exact solutions for α = 1. This ensures that our technique is accurate in the results for $\alpha =0.8,0.9.$ Moreover, the results for different values of α are presented to prove the ability of the method to provide accurate results. Finally, all results prove the validity and efficiency of these methods in solving nonlinear fractional differential equations.