Implementing Bernstein Operational Matrices to Solve a Fractional-Order Smoking Epidemic Model

Abstract

This paper leverages the Bernstein operational matrices method for the first time in order to resolve the nonlinear fractional smoking epidemic model presented in terms of Caputo’s fractional derivative. An approximate solution is derived using Bernstein’s operational matrices and strategically chosen collocation points. This is followed by the validation of the proposed method’s accuracy and reliability against the established Runge–Kutta fourth-order method. Furthermore, a comprehensive comparative analysis is conducted against two prominent techniques: the fractional differential transform method (FDTM) and the q-homotopy analysis transform method (q-HATM). The results show a superior and significant performance regarding accuracy as well as approximation. A residual corrected error technique is employed to enhance the precision of the presented method, thus effectively minimising absolute errors.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

The global response to infectious diseases highlights the importance of utilising mathematical epidemic models to comprehend and manage disease spread. As these models become more complex, the demand for advanced computing systems capable of handling this complexity with high accuracy and efficiency grows. Recently, numerous research papers have been published showcasing the latest advancements in computational techniques aimed at enhancing the reliability and performance of epidemic and biological modelling. These developments offer a crucial perspective on addressing current and future outbreaks [1–4].

Smoking denotes an example of an extremely hazardous behaviour that spreads through the population, like infectious diseases. It can cause cancer, chronic lung disease, cardiovascular disease and diabetes. According to the World Health Organization (WHO), tobacco use causes the deaths of over eight million people each year. In addition, 1.3 million nonsmokers are estimated to be exposed to secondhand smoke every year [5]. Approximately 80% of the world’s 1.3 billion tobacco users live in low- and middle-income countries. Data from 2020 show that 22.3% of the global population used tobacco, with usage being significantly higher among males (36.7%) compared with females (7.8%). As a result, numerical analysis modelling is often employed to study smoking behaviour.

In the year 1997, Garsow et al. [6] introduced a basic and simple model for smoking cessation so as to better understand the manner in which smoking spreads in a community. They divided the population into three distinct groups: potential smokers, represented by P; active smokers, signified by S; and quit smokers, denoted by Q, with the possibility of size changes over a period of time. Building upon this initial work, Zaman [7] expanded the smoking epidemic model through the inclusion of another group, occasional or light smokers. Zaman [8] further placed emphasis on controlling the smoking epidemic by considering optimal campaigns and proposing two control measures: treating smokers with medication and educating potential smokers about the hazards of smoking. The researcher formulated a control problem and exemplified the existence of an optimal control strategy, emphasising the prevalence of smoking can be reduced significantly by ensuring the proper implementation of both measures. Further investigations integrated supplementary parameters and smoking-related characteristics into the modelling approach.

A novel model for giving up smoking was proposed by Zeb et al. [9], who explored the local as well as global stability of the epidemic model, along with its general solutions. Subsequently, Erturk et al. [10] introduced fractional-order derivatives into the ordinary differential equation model by Zeb et al. [9] and employed generalised Euler as well as the generalised multistep differential transform. On the other hand, Khalid et al. [11] employed a perturbation–iteration algorithm in order to examine a smoking model consisting of five variables, incorporating the Caputo fractional derivative. Meanwhile, Haq et al. [12] implemented the Laplace transformation to determine numerical approximations for the smoking model with a Caputo fractional derivative. Furthermore, Singh et al. [13] established theorems with respect to the existence and uniqueness of a smoking cessation model equipped with the Caputo–Fabrizio derivative, employing the fixed-point postulate. Then, Uçar et al. [14] aimed to analyse the outcome of the Atangana–Baleanu derivative by addressing the smoking model introduced to delve into the existence and uniqueness of the fractionalised model. Recently, the q-homotopy analysis transform method (q-HATM) was utilised to approximate the fractional smoking model by Veeresha et al. [15] and then Günerhan et al. [16] used the fractional differential transform method (FDTM) to resolve the model, also comparing their results with q-HATM. Thereafter, Haidong et al. [17] investigated a new time noninteger-order smoking system of equations with a Caputo fractional-type derivative operator and examined its qualitative aspects. The model exhibited stability and convergence, with improved stability being showcased by smaller fractional orders displaying. Numerous other researchers have performed extensive investigations on the mathematical epidemic smoking model [18–26].

The model explored in this study is the fractional smoking epidemic model, designed to study and forecast how smoking affects individuals’ health and the population as a whole. The model divides the population into five distinct groups: potential smokers, occasional smokers, smokers, permanent quitters and temporary quitters. The investigation’s parameters include the recruiting rate, effective contact rate, natural death rate, smoking cessation rate and transition rates among smoking classes. Integrating fractional calculus into the model allows for the simulation of individual behaviour over time, incorporating memory effects and nonlocality. These variables are crucial for accurately representing physical activities such as smoking habits. The model provides a framework for studying and understanding how smoking habits evolve within a community, helping scholars address specific issues related to smoking trends. The results indicate that both the fractional order and the parameter system significantly influence the development of smoking. The solutions discovered for various categories of smokers are potential, occasional and permanent quitters, which offer meaningful insights into smoking habits over time and the impact of different factors on the model’s stability.

Numerous scientists have put forth various methodologies for evaluating smoking trends and the role of mathematical models in addressing this significant public health issue. The endeavour of this research is to provide an approximate solution for a fractional-order smoking pandemic model using the Bernstein operational matrices approach. This method assesses the impact of smoking on various population groups and compares the effectiveness of the Bernstein operational matrices technique to other well-established methods. This work lays the groundwork for future investigations into smoking behaviour and the development of strategies to reduce smoking-related illnesses and fatalities. Consequently, there will be comprehensive model improvements, an examination of additional factors influencing smoking behaviour and an expansion of the research to include diverse populations and geographical regions.

Fractional calculus has been an interest of many researchers and is extensively used in engineering and science, with numerous applications in electromagnetics, fluid mechanics and epidemiology [27–31]. Alshbool et al. [32] solved some mathematical systems comprising second-order BVPs, a Brusselator system, as well as a nonlinear stiff system using the Bernstein polynomials approach. This research extends the application of Bernstein polynomials to solve a fractional smoking epidemic system, offering several benefits. It reduces the time required to calculate approximate solutions, thereby accelerating numerical analysis and computation. At the same time, it decreases the computational workload compared to traditional methods whilst maintaining the accuracy and efficiency of the approximate solutions. This method is particularly powerful for analysing nonlinear fractional differential equations in complex problems across epidemiology, science and engineering.

This study provides unique insights into smoking cessation dynamics and presents a significant approximate solution. However, practical applications of these findings should consider other factors, such as the complexity of human behaviour and the social and environmental influences on smoking. The need for interdisciplinary collaboration is essential to translate the mathematical model into effective smoking-cessation programs. Therefore, while the study offers an intriguing framework for explaining smoking habit dynamics, its real-life application may require further research and collaboration with specialists in public health, psychiatry and behavioural sciences.

The structure of the paper is as follows. In Section 2, important definitions and preliminaries necessary for the study are presented. Section 3 presents the fractional order smoking epidemic mathematical model, the existence of the smoking-free equilibrium (SFE), the basic reproduction number $R_0$, and the stability of the endemic equilibrium of the smoking system. Section 4 illustrates the Bernstein operational matrices method. In Section 5, numerical approximation results of the fractional smoking epidemic model are presented, with comparisons made to FDTM and q-HATM from the literature, along with error estimates. Section 6 elaborates on the residual correction technique. Section 7 discusses the numerical findings of the smoking model. Finally, Section 8 provides a summary of the study’s findings.

2. Preliminaries

In this section, essential definitions and properties of fractional derivatives in Caputo’s sense that are needed throughout this paper to numerically solve the given smoking system of fractional order are introduced [33–36].

Definition 1.

The Caputo fractional derivative of a function $gt$ having an order of $\alpha >0$ is given by$$\begin{array}{}\text{(1)}& D^{\alpha }gt=\frac{1}{\Gamma s-\alpha }{\int }_0^x{t-u}^{s-\alpha -1}{g}^{s}tdu,\end{array}$$such that $s-1<\alpha <s,s\in N,t>0$, and $g\in C_{-1}^s$, where $C_{-1}^s$ is the space of functions with minimal smoothness, ensuring the Caputo derivative is well defined for $0\le \alpha \le s$.

For the Caputo fractional derivative, we have two essential properties defined as

a. For $g\in C_{-1}^s$, $s\in N$, we have ${\mathrm{D}}^{\alpha }g$ is well defined when $0\le \alpha \le s$ and ${\mathrm{D}}^{\alpha }g\in C_{-1}$.

b. For $s-1\le \alpha \le s$, $s\in N$, $g\in C_q^s$ for $q\ge -1$, where $C_q^s$ is a smoother function space ensuring continuity and differentiability for $s-1\le \alpha \le s$, we have

$$\begin{array}{}\text{(2)}& J^{\alpha }{\mathrm{D}}^{\alpha }gt=gt-{\displaystyle \sum }_{k=0}^{s-1}{g}^k{0}^{+}\frac{t^k}{k!}.\end{array}$$

In the case of the Caputo derivative, we have

a. ${\mathrm{D}}^{\alpha }k=0$, ($k$ constant),

b. ${\mathrm{D}}^{\alpha }{t}^{\psi }=\begin{array}{c}\begin{array}{cc}0,& \text{if\hspace{0.17em}}\psi \in N_0\text{\hspace{0.17em}and\hspace{0.17em}}\psi <\lceil \alpha \rceil ,\\ \Gamma \psi +1/\Gamma \psi +1-\alpha ,& \text{if\hspace{0.17em}}\psi \in N_0\text{\hspace{0.17em}and\hspace{0.17em}}\psi \ge \lceil \alpha \rceil \text{\hspace{0.17em}or\hspace{0.17em}}\psi >\lceil \alpha \rceil .\end{array}\end{array}$

Definition 2.

The Bernstein polynomial of order $n$ on the interval [0, 1] is$$\begin{array}{}\text{(3)}& B_{i,n}t=\begin{array}{c}n\\ i\end{array}{t}^i{1-t}^{n-i}.\end{array}$$

Employing the Binomial theorem in equation (3), we obtain$$\begin{array}{}\text{(4)}& B_{i,n}t={\displaystyle \sum }_{j=0}^{n-i}{-1}^j\begin{array}{c}n\\ i\end{array}\begin{array}{c}n-i\\ j\end{array}{t}^{i+j}.\end{array}$$

Definition 3.

The Bernstein polynomial of order $n$ on the interval $0,R$ is$$\begin{array}{}\text{(5)}& B_{i,n}t=\begin{array}{c}n\\ i\end{array}\frac{t^i{R-t}^{n-i}}{R^n}={\displaystyle \sum }_{j=0}^{n-i}{-1}^j\begin{array}{c}n\\ i\end{array}\begin{array}{c}n-i\\ j\end{array}\frac{t^{i+j}}{R^{i+j}},i=0,1,\dots ,n.\end{array}$$

Definition 4.

The Bernstein polynomial of order $n$ on the interval $a,b$ is$$\begin{array}{}\text{(6)}& B_{i,n}t=\begin{array}{c}n\\ i\end{array}\frac{{t-1}^i{b-t}^{n-i}}{{b-a}^n},i=0,1,\dots ,n.\end{array}$$

Substituting $t⟶t^{\alpha }$ in equation (3), we obtain the fractional Bernstein polynomial on $0,1$ as follows:$$\begin{array}{}\text{(7)}& B_{i,n}^{\alpha }t=\begin{array}{c}n\\ i\end{array}{t}^{i\alpha }{1-t^{\alpha }}^{n-i}={\displaystyle \sum }_{j=0}^{n-i}{-1}^j\begin{array}{c}n\\ i\end{array}\begin{array}{c}n-i\\ j\end{array}{t}^{i+j\alpha }.\end{array}$$

3. Analysis of the Fractional-Order Smoking Epidemic Model

In this section, we present the fractional-order smoking epidemic mathematical model, the existence of the SFE, the process of finding the basic reproduction number $R_0$ using the next-generation matrix method, and we assess the stability of the endemic equilibrium for the given fractional-order smoking model.

3.1. Fractional-Order Smoking Epidemic Model

This research investigates an existing fractional-order epidemic system of equations from literature [15, 16] used for simulating the dynamics of smoking. The Caputo fractional form incorporates the fractional derivative; the system is expressed as follows:$$\begin{array}{c}\text{(8)}& D^{\alpha }Pt=\omega -\sigma PtSt-\kappa Pt,D^{\alpha }Ot=\sigma PtSt-\beta Ot-\kappa Ot,\\ D^{\alpha }St=\beta Ot+\lambda StQt-\kappa +\eta St,\\ D^{\alpha }Qt=-\lambda StQt-\kappa Qt+\eta \cdot 1-\mu St,\\ D^{\alpha }Lt=\eta St\mu -\kappa Lt,\end{array}$$with initial conditions$$\begin{array}{c}\text{(9)}& P0={\delta }_1,O0={\delta }_2,\\ S0={\delta }_3,\\ Q0={\delta }_4,\\ L0={\delta }_5.\end{array}$$

In the system of equation (8), $D^{\alpha }$ denotes the operator for the fractional-order derivative of order $\alpha$ such that $0<\alpha \le 1$. At time $t$, the whole population size is represented as $Tt$ and splits into five distinct subgroups: potential smokers represented as $Pt$, occasional smokers as $Ot$, smokers as $St$, people who have quit smoking permanently as $Lt$ and those who have quit smoking temporarily as $Qt$.

The following parameters are accorded special significance in the given system of equation (8), where the recruitment rate is represented by $\omega$, the linkage index between light and potential smokers is signified by $\sigma$, the natural mortality rate is denoted by $\kappa$, the smoking cessation rate is represented by $\eta$, the smokers who never smoke again are signified by $\mu$, the rate between light and heavy smokers are denoted by $\beta$ and $\lambda$ showcases the linkage index between temporary quitters who might get back to smoking.

3.2. Existence of the SFE

The existence of the SFE is crucial, as it represents a state where the smoking prevalence reaches zero, indicating no further progression. It can be found by setting the time derivatives of the system (8) to zero at equilibrium.$$\begin{array}{c}\text{(10)}& D^{\alpha }Pt=0,D^{\alpha }Ot=0,\\ D^{\alpha }St=0,\\ D^{\alpha }Qt=0,\\ D^{\alpha }Lt=0.\end{array}$$

The SFE point for system (8) may be calculated directly as follows:$$\begin{array}{}\text{(11)}& E_0=\frac{\omega }{\kappa },0,0,0,0,\end{array}$$where the population at this point consists entirely of potential smokers, with no active smokers, occasional smokers, temporary quitters or permanent quitters [37].

3.3. Basic Reproduction Number $R_0$

In this section, the process of finding the basic reproduction number $R_0$ for the fractional-order smoking model using the next-generation matrix method is illustrated. The focus is on the infected compartments of the system, particularly the smokers $St$ and occasional smokers $Ot$, and linearise the system around the SFE to calculate $R_0$ [38–41]. The governing equations for $St$ and $Ot$ are$$\begin{array}{}\text{(12)}& D^{\alpha }Ot=\sigma PtSt-\beta +\kappa Ot,\text{(13)}& D^{\alpha }St=\beta Ot-\kappa +\eta St.\end{array}$$

At the SFE, $S_0=0$, $O_0=0$, and $P_0=\omega /k$. We linearise the system around this point to form the Jacobian matrix $J$. It is constructed by taking the partial derivatives of the system with respect to $St$ and $Ot$ as follows:$$\begin{array}{}\text{(14)}& J=\begin{array}{cc}-\beta +\kappa & \sigma \frac{\omega }{\kappa }\\ \beta & -\kappa +\eta \end{array}.\end{array}$$

To find the basic reproduction number $R_0$, we use the next-generation matrix method. The matrix $F$ represents new infections, and $V$ represents the transitions and removals. They are given by$$\begin{array}{c}\text{(15)}& F=\begin{array}{cc}0& \sigma \frac{\omega }{\kappa }\\ 0& 0\end{array},V=\begin{array}{cc}\beta +\kappa & 0\\ -\beta & \kappa +\eta \end{array}.\end{array}$$

The inverse of $V$ is calculated as$$\begin{array}{}\text{(16)}& V^{-1}=\frac{1}{\beta +\kappa \kappa +\eta }\begin{array}{cc}\kappa +\eta & 0\\ \beta & \beta +\kappa \end{array}.\end{array}$$

Multiply $F$ and $V^{-1}$ to get the next-generation matrix $G$ as follows:$$\begin{array}{}\text{(17)}& G=F{V}^{-1}=\frac{1}{\beta +\kappa \kappa +\eta }\begin{array}{cc}0& \sigma \frac{\omega }{\kappa }\\ 0& 0\end{array}\times \begin{array}{cc}\kappa +\eta & 0\\ \beta & \beta +\kappa \end{array}.\end{array}$$

This results in$$\begin{array}{}\text{(18)}& G=\frac{\sigma \omega }{\kappa \beta +\kappa \kappa +\eta }\begin{array}{cc}\beta & \beta +\kappa \\ 0& 0\end{array}.\end{array}$$

The basic reproduction number $R_0$ is the spectral radius of the next-generation matrix $G$ which is found as follows:$$\begin{array}{}\text{(19)}& R_0=\frac{\sigma \omega \beta }{\kappa \beta +\kappa \kappa +\eta }.\end{array}$$

If $R_0<1$, then $E_0$ is locally stable, indicating a decrease in smoking behaviour within the population. If $R_0>1$, then $E_0$ is unstable, indicating that smoking behaviour will spread among the population. Figure 1 illustrates a graphical phase diagram of $R_0$ in relation to $\sigma$ and $\omega$ for system (8), where $\kappa$, $\eta$, and $\beta$ are constant parameters from Table 1, highlighting the stable and unstable regions.

[figure(s) omitted; refer to PDF]

Table 1

Specified parameter values for equation (8).

Parameter	Description	Value (unit: 1/time)Source [15]
$\omega$	Recruitment rate in $Pt$	1
$\sigma$	Linkage index between $St$ and $Pt$	0.14
$\kappa$	Natural mortality rate	0.05
$\beta$	Rate at which $Ot$ become $St$	0.002
$\lambda$	Linkage index between $Qt$ who might get back to smoking	0.0025
$\eta$	Smoking cessation rate	0.8
$\mu$	Remaining fraction of smokers who permanently quit smoking	0.1

3.4. Endemic Equilibrium

In this section, we assess the stability of the endemic equilibrium for the fractional smoking model. Using the system of equations provided in the model, we first compute the endemic equilibrium points and then analyse the stability using the Jacobian matrix and eigenvalues.

The endemic equilibrium points for the smoking fractional model in steady state can be summarised as follows:$$\begin{array}{}\text{(20)}& E^{\mathrm{\ast }}=P^{\mathrm{\ast }},O^{\mathrm{\ast }},S^{\mathrm{\ast }},Q^{\mathrm{\ast }},L^{\mathrm{\ast }},\end{array}$$where$$\begin{array}{c}\text{(21)}& P^{\mathrm{\ast }}=\frac{\omega }{\sigma S^{\mathrm{\ast }}+\kappa },O^{\mathrm{\ast }}=\frac{\sigma \omega S^{\mathrm{\ast }}}{\sigma S^{\mathrm{\ast }}+\kappa \beta +\kappa },\\ Q^{\mathrm{\ast }}=\frac{\eta 1-\mu S^{\mathrm{\ast }}}{\lambda S^{\mathrm{\ast }}+\kappa },\\ L^{\mathrm{\ast }}=\frac{\eta \mu S^{\mathrm{\ast }}}{\kappa }.\end{array}$$

Using the parameter values from Table 1 and assuming $S^{\mathrm{\ast }}=1$ for simplicity, we compute the approximate endemic equilibrium values as follows.

To assess the stability of the endemic equilibrium, we compute the Jacobian matrix $J{E}^{\mathrm{\ast }}$ by taking the partial derivatives of the system’s equations with respect to $P$, $O$, $S$, $Q$, and $L$. The Jacobian at endemic equilibrium is given by$$\begin{array}{}\text{(22)}& J{E}^{\mathrm{\ast }}=\begin{array}{ccccc}-\sigma S^{\mathrm{\ast }}-\kappa & 0& -\sigma P^{\mathrm{\ast }}& 0& 0\\ \sigma S^{\mathrm{\ast }}& -\beta +\kappa & \sigma P^{\mathrm{\ast }}& 0& 0\\ 0& \beta & \lambda Q^{\mathrm{\ast }}-\kappa +\eta & \lambda S^{\mathrm{\ast }}& 0\\ 0& 0& \eta 1-\mu -\lambda Q^{\mathrm{\ast }}& -\lambda S^{\mathrm{\ast }}+\kappa & 0\\ 0& 0& \eta \mu & 0& -\kappa \end{array}.\end{array}$$

Using the parameter values from Table 1 and assuming $S^{\mathrm{\ast }}=1$, the Jacobian at the endemic equilibrium is given by$$\begin{array}{}\text{(23)}& J{E}^{\mathrm{\ast }}=\begin{array}{ccccc}-0.19& 0& -0.7368& 0& 0\\ 0.14& -0.052& 0.7368& 0& 0\\ 0& 0.002& -0.8157& 0.0343& 0\\ 0& 0& 0.6857& -0.0525& 0\\ 0& 0& 0.08& 0& -0.05\end{array}.\end{array}$$

The eigenvalues of the Jacobian matrix are$$\begin{array}{c}\text{(24)}& {\lambda }_1=-0.05,{\lambda }_2=-0.8203,\\ {\lambda }_3=-0.1877,\\ {\lambda }_4=-0.0503,\\ {\lambda }_5=-0.0520.\end{array}$$

For fractional-order systems, the stability criterion differs slightly from the classical integer-order case. The equilibrium is locally asymptotically stable if the real parts of all eigenvalues ${\lambda }_i$ satisfy$$\begin{array}{}\text{(25)}& \alpha \text{Re}{\lambda }_i<\frac{\pi }{2}.\end{array}$$

Given that all the real parts of the eigenvalues are negative, the endemic equilibrium $E^{\mathrm{\ast }}$ is locally asymptotically stable for any fractional order $0<\alpha \le 1$ [42].

4. Numerical Method

In this section, the Bernstein operational matrices approximation method [43] is presented using the Caputo fractional derivative properties and the collocation points to solve numerically the nonlinear system of fractional order of the epidemic smoking model presented in equation (8).

We start by considering the system of fractional differential equations as follows:$$\begin{array}{}\text{(26)}& {\mathrm{D}}^{\alpha }{g}_{j}t=y_{j}t,g_1,g_{2,}\dots ,g_k,\end{array}$$subject to the following initial conditions:$$\begin{array}{}\text{(27)}& g_j{t}_0={\theta }_j,\end{array}$$where ${\theta }_j$ are constants and $j=1,2,\dots ,k$. Our initial approach is to rewrite equation (26) in the following form:$$\begin{array}{}\text{(28)}& {\mathrm{D}}^{\alpha }{g}_{j}t-y_{j}t,g_1,g_{2,}\dots ,g_k=0,\end{array}$$subject to the initial conditions (27).

The Bernstein polynomials of the $n$ th degree are defined by$$\begin{array}{}\text{(29)}& B_{i,n}t=\begin{array}{c}n\\ i\end{array}\frac{t^i{R-t}^{n-i}}{R^n},i=0,1,2,\dots ,n,t\in 0,R.\end{array}$$

To change equation (29) to fractional Bernstein polynomial, we substitute $t⟶t^{\alpha }$ into $B_{i,n}t$ to obtain the following:$$\begin{array}{}\text{(30)}& B_{i,n}^{\alpha }t=\begin{array}{c}n\\ i\end{array}\frac{t^{i\alpha }{R-t^{\alpha }}^{n-i}}{R^n},0<\alpha <1,\hspace{1em}i=0,1,2,\dots ,n,t\in 0,R.\end{array}$$

For mathematical convenience, we set the following $B_{i,n}^{\alpha }t=0$, if $i<0$ or $i>n$. In the present method, we approximate any function $g_{j}t$ with the first $n+1$ Bernstein polynomials as follows:$$\begin{array}{}\text{(31)}& g_{j}t={\displaystyle \sum }_{i=0}^n{c}_{j,i}{B}_{i,n}^{\alpha }t=C_j^T\varphi t,j=1,2,\dots ,k,\end{array}$$where$$\begin{array}{}\text{(32)}& \varphi t=B_{0,n}^{\alpha }t,B_{1,n}^{\alpha }t,\dots ,B_{n,n}^{\alpha }t.\end{array}$$

We express the vector $\varphi t$ as follows:$$\begin{array}{}\text{(33)}& \varphi t=AT,\end{array}$$where$$\begin{array}{}\text{(34)}& A=\begin{array}{ccccc}{a}_{0,0}& a_{0,1}& a_{0,2}& \dots & a_{0,k}\\ a_{1,0}& a_{1,1}& a_{1,2}& \dots & a_{1,k}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ a_{k0}& a_{k,1}& a_{k,2}& \dots & a_{k,k}\end{array},\end{array}$$such that $a_{i,j}=\begin{array}{cc}{-1}^{j-i}/R^j\begin{array}{c}k\\ i\end{array}\begin{array}{c}k-i\\ j-i\end{array},& i\le j,\\ 0,& i>j,\end{array}$$$\begin{array}{}\text{(35)}& T=\begin{array}{c}1\\ t\\ t^2\\ ⋮\\ t^k\end{array}.\end{array}$$

For $j=1$, we have$$\begin{array}{}\text{(36)}& g_{1}t={\displaystyle \sum }_{i=0}^n{c}_{1,i}{B}_{i,n}^{\alpha }t=C_1^T\varphi t,\end{array}$$where$$\begin{array}{}\text{(37)}& C_1^T=\begin{array}{cccc}{C}_{1,0}& C_{1,1}& \dots & C_{1,n}\end{array},\end{array}$$etc. For $j=k$, we have$$\begin{array}{}\text{(38)}& g_{k}t={\displaystyle \sum }_{i=0}^n{c}_{n,i}{B}_{i,n}^{\alpha }t=C_{\mathrm{k}}^T\varphi t,\end{array}$$where$$\begin{array}{}\text{(39)}& C_k^T=\begin{array}{cccc}{C}_{k,0}& C_{k,1}& \dots & C_{k,n}\end{array}.\end{array}$$

The vector $\varphi t$ can be expressed as a fractional derivative as follows:$$\begin{array}{}\text{(40)}& \frac{d^{\alpha }\varphi t}{dt}=D^{\alpha }\varphi t,\end{array}$$where $D^{\alpha }$ is $n+1\times n+1$ operational matrix of derivative defined as$$\begin{array}{}\text{(41)}& {\mathrm{D}}^{\alpha }=AMN{A}^{-1},\end{array}$$given that$$\begin{array}{c}\text{(42)}& A=\begin{array}{ccccc}{a}_{0,0}& a_{0,1}& a_{0,2}& \dots & a_{0,k}\\ a_{1,0}& a_{1,1}& a_{1,2}& \dots & a_{1,k}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ a_{k0}& a_{k,1}& a_{k,2}& \dots & a_{k,k}\end{array},M=\begin{array}{cccccc}0& 0& 0& 0& 0& 0\\ 0& \frac{\Gamma \alpha +1}{\Gamma 1}& 0& 0& 0& 0\\ 0& 0& \frac{\Gamma 2\alpha +1}{\Gamma \alpha +1}& 0& 0& 0\\ 0& 0& 0& \frac{\Gamma 3\alpha +1}{\Gamma 2\alpha +1}& 0& 0\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& 0& 0& \dots & \frac{\Gamma k\alpha +1}{\Gamma k-1\alpha +1}\end{array},\\ N=\begin{array}{cccccc}0& 0& 0& 0& 0& 0\\ 1& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& 0& 0& 1& 0\end{array}.\end{array}$$

By applying equations (31) and (40) in system (28), we obtain the residual $\mathfrak{R}t$ which measures the error in approximating the fractional differential equation as follows:$$\begin{array}{c}\text{(43)}& C_1^T{D}^{\alpha }\varphi t-y_{1}t,C_1^T\varphi t,C_2^T\varphi t,\dots ,C_k^T\varphi t=0,C_2^T{D}^{\alpha }\varphi t-y_{2}t,C_1^T\varphi t,C_2^T\varphi t,\dots ,C_k^T\varphi t=0,\\ ⋮\\ C_k^T{D}^{\alpha }\varphi t-y_{k}t,C_1^T\varphi t,C_2^T\varphi t,\dots ,C_k^T\varphi t=0,\end{array}$$with initial conditions$$\begin{array}{c}\text{(44)}& C_1^T\varphi t_0={\theta }_1,C_2^T\varphi t_0={\theta }_2\\ \dots \\ C_k^T\varphi t_0={\theta }_k.\end{array}$$

To avoid integration difficulty, the collocation points $t_0,t_1,\dots ,t_k$ in system (43) can be determined by implementing the Chebyshev roots as follows:$$\begin{array}{}\text{(45)}& t_i=0.5+\frac{1}{2\cos 2i+1\pi /2n},i=0,1,\dots ,n-1.\end{array}$$

Consequently, solving the system of equation (43) can be done by finding the collocation points using the equation (45) and using the given initial conditions of (44) to calculate all the $c_{j,i}$. Applying all this process will enable us to find approximate solutions $g_{1}t,g_{2}t,\dots ,g_{k}t$ of the system of fractional order (26).

5. Numerical Results and Error Estimates

In this section, the application of the illustrated method in Section 4 on the given system of equations of fractional order for the smoking epidemic model presented in equation (8) with different parameters is demonstrated. Moreover, a comparison is made between the method in Section 4 with two recent methods applied to the same equation (8), i.e., q-HATM [15] and FDTM [16].

For $k=5$ and $n=7$, we start implementing our presented method by letting$$\begin{array}{c}\text{(46)}& Pt=g_{1}t=C_1^T\varphi t,Ot=g_{2}t=C_2^T\varphi t,\\ St=g_{3}t=C_3^T\varphi t,\\ Qt=g_{4}t=C_4^T\varphi t,\\ Lt=g_{5}t=C_5^T\varphi t.\end{array}$$

We substitute all the variables in (46) in the system of smoking model (8), then we find residual $\mathfrak{R}\mathrm{t}$ as follows:$$\begin{array}{c}\text{(47)}& C_1^T{D}^{\alpha }\varphi t-\omega +\sigma C_1^T\varphi t{C}_3^T\varphi t+\kappa C_1^T\varphi t=0,C_2^T{D}^{\alpha }\varphi t-\sigma C_1^T\varphi t{C}_3^T\varphi t+\beta C_2^T\varphi t+\kappa C_2^T\varphi t=0,\\ C_3^T{D}^{\alpha }\varphi t-\beta C_2^T\varphi t-\lambda C_3^T\varphi t{C}_4^T\varphi t-\kappa +\eta C_3^T\varphi t=0,\\ C_4^T{D}^{\alpha }\varphi t+\lambda C_3^T\varphi t{C}_4^T\varphi t+\kappa C_4^T\varphi t-\eta \cdot 1-\mu C_3^T\varphi t=0,\\ C_5^T{D}^{\alpha }\varphi t-\eta C_3^T\varphi t\mu +\kappa C_5^T\varphi t=0,\end{array}$$subject to the following initial conditions:$$\begin{array}{c}\text{(48)}& C_1^T\varphi 0={\delta }_1=40,C_2^T\varphi 0={\delta }_2=10,\\ C_3^T\varphi 0={\delta }_3=20,\\ C_4^T\varphi 0={\delta }_4=10,\\ C_5^T\varphi 0={\delta }_5=5.\end{array}$$

Specified values for the parameters in the system of equation (8) are given in Table 1 to verify the effectiveness and precision of the method used in this research.

By finding the collocation points using equation (45) and using the given initial conditions in (48), we calculate all the $c_{j,i}$ with the help of Maple Software. Then, we find approximate solutions of $Pt$, $Ot$, $St$, $Qt$ and $Lt$. The efficiency and accuracy of the presented method are verified when comparing it with the well-known Runge–Kutta fourth-order method (RK4).

The following graphs show the approximate solutions of $Pt$, $Ot$, $St$, $Qt$ and $Lt$ with different values of $\alpha$, as shown in Figure 2. The accuracy and efficiency of the presented method are clearly seen in all the graphs. Moreover, Tables 2, 3, 4, 5 and 6 show the comparisons of the approximated solutions of the Bernstein operational matrices method with RK4, FDTM, and q-HATM for the given smoking epidemic model 8 at $\alpha =1$ and $\alpha =0.9$. The results show strong evidence for the validity of the presented method.

[figure(s) omitted; refer to PDF]

Table 2

Approximate solutions of $Pt$ at $n=7$ when $\alpha =1$ and $\alpha =0.9$ for Bernstein operational matrices, FDTM and q-HATM.

$Pt$
$t$	RK4	$\alpha =1$			$\alpha =0.9$
		Present method	FDTM	q-HATM	Present method	FDTM	q-HATM
0	40	40	40	40	40	40	40
0.1	30.5072	30.5063	30.9766	30.7667	30.3581	26.2807	29.1162
0.2	23.7850	23.7853	24.174	25.6668	23.8256	17.8878	26.0052
0.3	18.9262	18.9271	19.2495	24.7002	19.2628	12.8287	28.4733
0.4	15.3484	15.3480	15.6184	27.867	15.9889	9.6387	35.8752
0.5	12.6687	12.6675	12.8957	35.1672	13.5890	7.5495	47.8402
0.6	10.6306	10.6304	10.8229	46.6008	11.7965	6.1373	64.1152
0.7	9.0588	9.0596	9.223	62.1678	10.4314	5.158	84.5111
0.8	7.8313	7.8315	7.9727	81.8682	9.3728	4.4654	108.878
0.9	6.8620	6.8614	6.9847	105.702	8.5440	3.9687	137.095
1	6.0889	6.0890	6.1962	133.669	7.8904	3.6099	169.057

Table 3

Approximate solutions of $Ot$ at $n=7$ when $\alpha =1$ and $\alpha =0.9$ for Bernstein operational matrices, FDTM and q-HATM.

$Ot$
$t$	RK4	$\alpha =1$			$\alpha =0.9$
		Present method	FDTM	q-HATM	Present method	FDTM	q-HATM
0	10	10	10	10	10	10	10
0.05	15.0909	15.0913	14.5834	15.0559	15.2268	17.3313	16.6946
0.1	19.3402	19.3411	18.8789	19.0756	19.4832	23.474	20.6741
0.15	22.9048	22.9053	22.4858	22.059	22.9735	28.0928	22.8892
0.2	25.9090	25.9087	25.5286	24.0063	25.8561	31.6072	23.5858
0.25	28.4522	28.4513	28.1067	24.9173	28.2532	34.3096	22.8977
0.3	30.6140	30.6131	30.2999	24.7922	30.2590	36.4067	20.9135
0.35	32.4585	32.4582	32.1728	23.6308	31.9467	38.047	17.6978
0.4	34.0380	34.0385	33.7777	21.4332	33.3732	39.3381	13.3004
0.45	35.3950	35.3959	35.1575	18.1994	34.5835	40.3594	7.76208
0.5	36.5641	36.5653	36.3473	13.9294	35.6138	41.1701	1.11636

Table 4

Approximate solutions of $St$ at $n=7$ when $\alpha =1$ and $\alpha =0.9$ for Bernstein operational matrices, FDTM and q-HATM.

$St$
$t$	RK4	$\alpha =1$			$\alpha =0.9$
		Present method	FDTM	q-HATM	Present method	FDTM	q-HATM
0	20	20	20	20	20	20	20
0.1	18.4221	18.4221	18.5009	18.4244	18.3696	17.5281	17.9798
0.2	16.9756	16.9756	17.0509	16.9938	16.8930	15.2715	16.4515
0.3	15.6484	15.6484	15.7202	15.708	15.5537	13.3178	15.1913
0.4	14.4298	14.4298	14.4982	14.5671	14.3372	11.6235	14.1491
0.5	13.3102	13.3102	13.3753	13.5711	13.2309	10.1521	13.2994
0.6	12.2810	12.2810	12.3429	12.7201	12.2235	8.8727	12.626
0.7	11.3344	11.3344	11.3932	12.0138	11.3053	7.7594	12.1178
0.8	10.4633	10.4633	10.5192	11.4525	10.4675	6.7897	11.7662
0.9	9.6615	9.6615	9.7144	11.036	9.7023	5.9444	11.5645
1	8.9231	8.9231	8.9733	10.7645	9.0028	5.2073	11.5073

Table 5

Approximate solutions of $Qt$ at $n=7$ when $\alpha =1$ and $\alpha =0.9$ for Bernstein operational matrices, FDTM and q-HATM.

$Qt$
$t$	RK4	$\alpha =1$			$\alpha =0.9$
		Present method	FDTM	q-HATM	Present method	FDTM	q-HATM
0	10	10	10	10	10	10	10
0.1	11.2780	11.2780	11.2142	11.276	11.3198	11.9966	11.6331
0.2	12.4398	12.4398	12.3791	12.4241	12.5034	13.7923	12.852
0.3	13.4959	13.4959	13.4382	13.4443	13.5654	15.3195	13.8408
0.4	14.4556	14.4556	14.401	14.3365	14.5184	16.6167	14.6413
0.5	15.3276	15.3276	15.2759	15.1008	15.3737	17.7164	15.2749
0.6	16.1194	16.1194	16.0705	15.7372	15.3737	18.6462	15.755
0.7	16.8381	16.8381	16.7919	16.2456	16.8297	19.4293	16.0913
0.8	17.4898	17.4898	17.4462	16.626	17.4471	20.0859	16.2908
0.9	18.0802	18.0802	18.0392	16.8786	18.0003	20.6329	16.3594
1	18.6145	18.6145	18.5759	17.0032	18.4954	21.0851	16.3015

Table 6

Approximate solutions of $Lt$ at $n=7$ when $\alpha =1$ and $\alpha =0.9$ for Bernstein operational matrices, FDTM and q-HATM.

$Lt$
$t$	RK4	$\alpha =1$			$\alpha =0.9$
		Present method	FDTM	q-HATM	Present method	FDTM	q-HATM
0	5	5	5	5	5	5	5
0.1	5.1283	5.1283	5.1218	5.128	5.1324	5.1998	5.1636
0.2	5.2438	5.2438	5.2377	5.2422	5.2499	5.3765	5.2841
0.3	5.3478	5.3478	5.342	5.3426	5.3540	5.5235	5.3803
0.4	5.4410	5.4410	5.4356	5.4291	5.4460	5.6446	5.4564
0.5	5.5245	5.5245	5.5194	5.5017	5.5271	5.7434	5.5148
0.6	5.5990	5.5990	5.5942	5.5605	5.5983	5.8229	5.5567
0.7	5.6653	5.6653	5.6608	5.6054	5.6606	5.8856	5.5832
0.8	5.7239	5.7239	5.7197	5.6365	5.7149	5.9338	5.595
0.9	5.7756	5.7756	5.7717	5.6537	5.7618	5.9693	5.5928
1	5.8209	5.8209	5.8173	5.657	5.8021	5.9939	5.577

Tables 7, 8, and 9 compare the absolute error of $Pt$, $Ot$, $St$, $Qt$ and $Lt$ when $\alpha =1$ and $n=7$ for Bernstein operational method, FDTM, q-HATM and RK4. The calculations show that the absolute error for all the functions in our presented method is less than the other two methods (FDTM and q-HATM).

Table 7

Absolute error of P(t) and S(t) at $n=7$ when $\alpha =1$ for Bernstein operational, FDTM, q-HATM and RK4.

$t$	Absolute error of $Pt$			Absolute error of $St$
	Present method	FDTM	q-HATM	Present method	FDTM	q-HATM
0	0	0	0	0	0	0
0.1	$9.0\times {10}^{-4}$	$4.7\times {10}^{-1}$	$2.6\times {10}^{-1}$	$3.3\times {10}^{-7}$	$7.9\times {10}^{-2}$	$2.3\times {10}^{-3}$
0.2	$2.5\times {10}^{-4}$	$3.9\times {10}^{-1}$	$1.9\times {10}^0$	$2.4\times {10}^{-7}$	$7.5\times {10}^{-2}$	$1.8\times {10}^{-2}$
0.3	$8.9\times {10}^{-4}$	$3.2\times {10}^{-1}$	$5.8\times {10}^0$	$4.0\times {10}^{-7}$	$7.2\times {10}^{-2}$	$6.0\times {10}^{-2}$
0.4	$4.1\times {10}^{-4}$	$2.7\times {10}^{-1}$	$1.2\times {10}^1$	$2.4\times {10}^{-7}$	$6.8\times {10}^{-2}$	$1.4\times {10}^{-1}$
0.5	$1.1\times {10}^{-3}$	$2.3\times {10}^{-1}$	$2.3\times {10}^1$	$4.4\times {10}^{-7}$	$6.5\times {10}^{-2}$	$2.6\times {10}^{-1}$
0.6	$1.4\times {10}^{-4}$	$1.9\times {10}^{-1}$	$3.6\times {10}^1$	$2.0\times {10}^{-7}$	$6.2\times {10}^{-2}$	$4.4\times {10}^{-1}$
0.7	$8.4\times {10}^{-4}$	$1.6\times {10}^{-1}$	$5.3\times {10}^1$	$7.1\times {10}^{-7}$	$5.9\times {10}^{-2}$	$6.8\times {10}^{-2}$
0.8	$1.6\times {10}^{-4}$	$1.4\times {10}^{-1}$	$7.4\times {10}^{-1}$	$3.2\times {10}^{-7}$	$5.6\times {10}^{-2}$	$9.9\times {10}^{-1}$
0.9	$5.8\times {10}^{-4}$	$1.2\times {10}^{-1}$	$9.9\times {10}^1$	$1.2\times {10}^{-7}$	$5.3\times {10}^{-2}$	$1.4\times {10}^0$
1	$9.8\times {10}^{-5}$	$1.1\times {10}^{-1}$	$1.3\times {10}^2$	$8.4\times {10}^{-7}$	$5.0\times {10}^{-2}$	$1.8\times {10}^0$

Table 8

Absolute error of $Qt$ and $Lt$ at $n=7$ when $\alpha =1$ for Bernstein operational, FDTM, q-HATM and RK4.

$t$	Absolute error of $Qt$			Absolute error of $Lt$
	Present method	FDTM	q-HATM	Present method	FDTM	q-HATM
0	0	0	0	0	0	0
0.1	$1.1\times {10}^{-7}$	$6.4\times {10}^{-2}$	$2.0\times {10}^{-3}$	$6.9\times {10}^{-9}$	$6.5\times {10}^{-3}$	$3\times {10}^{-4}$
0.2	$4.4\times {10}^{-8}$	$6.0\times {10}^{-2}$	$1.6\times {10}^{-2}$	$3.8\times {10}^{-9}$	$6.1\times {10}^{-3}$	$1.6\times {10}^{-3}$
0.3	$7.8\times {10}^{-7}$	$5.8\times {10}^{-2}$	$5.2\times {10}^{-2}$	$5.0\times {10}^{-9}$	$5.8\times {10}^{-3}$	$5.2\times {10}^{-3}$
0.4	$8.6\times {10}^{-8}$	$5.5\times {10}^{-2}$	$1.2\times {10}^{-1}$	$1.1\times {10}^{-8}$	$5.4\times {10}^{-3}$	$1.2\times {10}^{-2}$
0.5	$1.8\times {10}^{-7}$	$5.2\times {10}^{-2}$	$2.3\times {10}^{-1}$	$2.5\times {10}^{-8}$	$5.1\times {10}^{-3}$	$2.3\times {10}^{-2}$
0.6	$4.5\times {10}^{-8}$	$4.9\times {10}^{-2}$	$3.8\times {10}^{-2}$	$2.9\times {10}^{-9}$	$4.8\times {10}^{-3}$	$3.9\times {10}^{-2}$
0.7	$1.2\times {10}^{-7}$	$4.6\times {10}^{-2}$	$5.9\times {10}^{-1}$	$6.9\times {10}^{-8}$	$4.5\times {10}^{-3}$	$6.0\times {10}^{-2}$
0.8	$1.3\times {10}^{-7}$	$4.4\times {10}^{-2}$	$8.6\times {10}^{-1}$	$8.0\times {10}^{-8}$	$4.2\times {10}^{-3}$	$8.7\times {10}^{-2}$
0.9	$2.7\times {10}^{-7}$	$4.1\times {10}^{-2}$	$1.2\times {10}^0$	$1.3\times {10}^{-7}$	$3.9\times {10}^{-3}$	$1.2\times {10}^{-1}$
1	$8.5\times {10}^{-7}$	$3.9\times {10}^{-2}$	$1.6\times {10}^0$	$1.6\times {10}^{-7}$	$3.6\times {10}^{-3}$	$1.6\times {10}^{-1}$

Table 9

Absolute error of $Ot$ at $n=7$ when $\alpha =1$ for Bernstein operational, FDTM, q-HATM and RK4.

Absolute error of $Ot$
$t$	Present method	FDTM	q-HATM
0	0	0	0
0.05	$4.0\times {10}^{-4}$	$5.1\times {10}^{-1}$	$3.5\times {10}^{-2}$
0.1	$9.2\times {10}^{-4}$	$4.6\times {10}^{-1}$	$2.6\times {10}^{-1}$
0.15	$5.4\times {10}^{-4}$	$4.2\times {10}^{-1}$	$8.5\times {10}^{-1}$
0.2	$8.9\times {10}^{-4}$	$3.8\times {10}^{-1}$	$1.9\times {10}^0$
0.25	$8.9\times {10}^{-4}$	$3.5\times {10}^{-1}$	$3.5\times {10}^0$
0.3	$4.0\times {10}^{-4}$	$3.1\times {10}^{-1}$	$5.8\times {10}^0$
0.35	$3.5\times {10}^{-4}$	$2.9\times {10}^{-1}$	$8.8\times {10}^0$
0.4	$4.1\times {10}^{-4}$	$2.6\times {10}^{-1}$	$1.3\times {10}^1$
0.45	$9.9\times {10}^{-4}$	$2.4\times {10}^{-1}$	$1.7\times {10}^1$
0.5	$1.1\times {10}^{-3}$	$2.2\times {10}^{-1}$	$2.3\times {10}^1$

6. Residual Correction Technique

In this section, we introduce a residual correction technique to approximate the solutions for the system of differential equations by minimising the absolute error to the greatest extent possible.

Let $g_{j}t$ and $g_{j_n}t$ be the exact and the approximate solutions, respectively, of equation (26). In the following technique, a residual correction can be applied to improve the estimations of the absolute error.

We start the residual correction technique by writing equation (26) as follows:$$\begin{array}{}\text{(49)}& g_{j_n}^{\alpha }t=y_{j}t,g_1,g_2,\dots ,g_k.\end{array}$$

Let $R_j$ be defined as follows:$$\begin{array}{}\text{(50)}& R_j=g_{j_n}^{\alpha }t-y_{j}t,g_{1_n},g_{2_n},\dots ,g_{k_n}.\end{array}$$

Adding and subtracting $R_j$ from equation (26) produces the following equation for the absolute error.$$\begin{array}{}\text{(51)}& e_{j_n}^{\alpha }t=y_{j}t,e_{1_n},e_{2_n},\dots ,e_{k_n}.\end{array}$$

We solve equation (51) using the presented method in Section 4 with degree $m$ such that $m>n$, where $e_{j_n}t=g_{j}t-g_{j_n}t$ having the following initial condition $e_{j_n}^{\alpha }{t}_0=0$.

Corollary 1.

Let $g_{j_n}t$ and $E_{j_m}t$ be the approximate solutions of equations (26) and (51), respectively. Then, $g_{j_n}t+E_{j_m}t$ is considered the corrected approximate solution of equation (26) having the following error function $e_{j_n}t-E_{j_m}t$. Taking into consideration that if ${e_{j_n}t+E_{j_m}t}_{\infty }<ϵ$, then $E_{j_m}t$ is the estimation of the absolute error. Furthermore, if ${e_{j_n}t+E_{j_m}t}_{\infty }<g_{j}t-g_{j_n}t$, we can conclude that $g_{j_n}t+E_{j_m}t$ is more precise solution compared with $g_{j_n}t$ in any given norm. The optimal $n$ in the presented method that finds the minimal absolute error might be found by $E_{j_m}t$ at any given norm.

Corollary 2.

To find the approximate solutions for different values on $n$, we start first by using the triangle inequality that implies the following.$$\begin{array}{}\text{(52)}& g_{j}t-g_{n_1}t-g_{j}t-g_{n_2}t⩽g_{n_1}t-g_{n_2}t.\end{array}$$

If the prior errors are sufficiently spread out, we may estimate a rough upper bound for the resultant error. To test the upper bound, we perform the following steps: if the error sequence is increasing (or decreasing), then$$\begin{array}{}\text{(53)}& gt-g_{n+1}t-gt-g_{n}t=1-Kgt-g_{n}t⩽g_{n+1}t-g_{n}t,\end{array}$$or$$\begin{array}{}\text{(54)}& gt-g_{n+1}t<gt-g_{n}t⩽\frac{1}{1-K}{g}_{n+1}t-g_{n}t,\end{array}$$where $gt-g_{n+1}t=Kgt-g_{n}t$$0⩽t<1$.

If $K>1$, then both absolute errors can be considered very good approximations.

The efficiency residual correction technique is shown in Figure 3 where the absolute error function at $n=3$ is much higher compared with the corrected absolute error function at $n=3$ and $m=8$ of the Bernstein operational matrices method at $\alpha =1$ for $St$, $Qt$ and $Lt$, which shows the applicability, efficiency and accuracy of the corrected absolute error technique.

[figure(s) omitted; refer to PDF]

7. Discussion

The numerical results reveal distinct behaviours among the subgroups $Pt$, $Ot$, $St$, $Qt$ and $Lt$ as the fractional order $\alpha$ varies, providing meaningful insights into smoking behaviour modelling. For $\alpha =1$, where memory effects are minimal, the decline in potential smokers $Pt$ and smokers $St$ is rapid, while permanent quitters $Lt$ increase quickly, indicating that interventions yield faster results. In this case, occasional smokers $Ot$ and temporary quitters $Qt$ stabilise quickly, reflecting faster relapse and cessation dynamics. As $\alpha$ decreases (i.e., $\alpha <1$), the model incorporates stronger memory effects, leading to slower declines in $Pt$ and $St$ with more gradual transitions between states. This suggests that smoking behaviour becomes more persistent and deeply rooted, making cessation efforts more challenging. Prolonged fluctuations in $Ot$ and $Qt$ further highlight the potential for relapse over time, underscoring the importance of sustained public health interventions. The reproduction number $R_0$ plays a crucial role in determining the system’s stability, with $R_0<1$ predicting a decline in smoking prevalence and $R_0>1$ indicating sustained or growing smoking behaviour. This implies that policies must aim to maintain $R_0<1$ by promoting smoking cessation and reducing initiation rates. In addition, the comparison between the Bernstein operational matrices method and other approaches—such as Runge–Kutta, FDTM and q-HATM—demonstrates the Bernstein method’s superior accuracy, efficiency and reduced computational error. These results validate the utility of Bernstein polynomials in solving complex nonlinear fractional models, reinforcing their value for studying public health dynamics.

8. Conclusions

This paper utilised the Bernstein operational matrices method to numerically solve a nonlinear system of fractional differential equations. By applying this method to the smoking epidemic model, we demonstrated its effectiveness and significance in this context. Our results underscored the method’s efficiency when compared with traditional approaches such as Runge–Kutta, emphasising the importance of integrating advanced techniques such as the Bernstein operational matrices method in modelling real-world problems. Furthermore, our findings revealed superior approximation outcomes compared with FDTM and q-HATM. In addition, the residual corrected error technique introduced in this research efficaciously minimised absolute error functions by refining previous approximation solutions for the given smoking model. The reliability of the model’s long-term predictions, derived from almost a decade of numerical analysis, is generally strong, particularly in capturing key trends in smoking behaviour over time. However, several factors could significantly alter these predictions, such as changes in public health policies, technological advances in cessation methods, the introduction of new tobacco products like e-cigarettes or shifts in societal attitudes towards smoking. These potential influences highlight the importance of continuously updating the model to maintain its accuracy and relevance. The presented method effectively captures the memory effect in smoking dynamics; however, future research could address constraints such as external influences affecting smoking behaviour and use sensitivity analysis to understand how parameter changes impact results. The application of machine learning for parameter estimation could further enhance model accuracy by adapting to real-time data. In addition, studying the social and economic aspects that influence smoking would improve the mathematical model and its use in public health planning. Finally, the successful implementation of the presented method in addressing epidemic models generates new research opportunities for precise simulation and analysis of complex systems across many disciplines of study.

Funding

The research was supported entirely by personal funds.

References

[1] W. Alhejili, M. M. Khader, K. Lotfy, A. A. El-Bary, M. Adel, "Studying of the Covid-19 Model by Using the Finite Element Method: Theoretical and Numerical Simulation," Soft Computing, vol. 28 no. 6, pp. 5263-5273, DOI: 10.1007/s00500-023-09310-6, 2024.

[2] M. Adel, M. M. Khader, T. A. Assiri, W. Kallel, "Numerical Simulation for COVID-19 Model Using a Multidomain Spectral Relaxation Technique," Symmetry, vol. 15 no. 4,DOI: 10.3390/sym15040931, 2023.

[3] Y. F. Ibrahim, S. E. Abd El-Bar, M. M. Khader, M. Adel, "Studying and Simulating the Fractional Covid-19 Model Using an Efficient Spectral Collocation Approach," Fractal and Fractional, vol. 7 no. 4,DOI: 10.3390/fractalfract7040307, 2023.

[4] A. F. S. Aboubakr, G. M. Ismail, M. M. Khader, M. A. E. Abdelrahman, A. M. T. AbdEl-Bar, M. Adel, "Derivation of an Approximate Formula of the Rabotnov Fractional-Exponential Kernel Fractional Derivative and Applied for Numerically Solving the Blood Ethanol Concentration System," AIMS Mathematics, vol. 8 no. 12, pp. 30704-30716, DOI: 10.3934/math.20231569, 2023.

[5] World Health Organization, World Health Organization: Tobacco, 2023.

[6] C. C. Garsow, G. J. Salivia, A. R. Herrera, "Mathematical Models for the Dynamics of Tobacco Use, Recovery and Relapse," Technical Report Series BU-1505-M, 1997.

[7] G. Zaman, "Qualitative Behavior of Giving up Smoking Models," Bulletin of the Malaysian Mathematical Sciences Society, vol. 34 no. 2, pp. 403-415, 2011.

[8] G. Zaman, "Optimal Campaign in the Smoking Dynamics," Computational and Mathematical Methods in Medicine, vol. 2011, pp. 163834-163839, DOI: 10.1155/2011/163834, 2011.

[9] A. Zeb, G. Zaman, S. Momani, "Square-root Dynamics of a Giving up Smoking Model," Applied Mathematical Modelling, vol. 37 no. 7, pp. 5326-5334, DOI: 10.1016/j.apm.2012.10.005, 2013.

[10] V. S. Erturk, G. Zaman, B. Alzalg, A. Zeb, S. Momani, "Comparing Two Numerical Methods for Approximating a New Giving up Smoking Model Involving Fractional Order Derivatives," Iranian Journal of Science and Technology, Transactions A: Science, vol. 41 no. 3, pp. 569-575, DOI: 10.1007/s40995-017-0278-x, 2017.

[11] M. Khalid, F. Sami, A. Iqbal, "Perturbation-Iteration Algorithm to Solve Fractional Giving Up Smoking Mathematical Model," International Journal of Computer Applications, vol. 142 no. 9,DOI: 10.5120/ijca2016909891, 2016.

[12] F. Haq, K. Shah, G. ur Rahman, M. Shahzad, M. Shahzad, "Numerical Solution of Fractional Order Smoking Model via Laplace Adomian Decomposition Method," Alexandria Engineering Journal, vol. 57 no. 2, pp. 1061-1069, DOI: 10.1016/j.aej.2017.02.015, 2018.

[13] J. Singh, D. Kumar, M. A. Qurashi, D. Baleanu, "A New Fractional Model for Giving Up Smoking Dynamics," Advances in Difference Equations, vol. 2017 no. 1, pp. 88-16, DOI: 10.1186/s13662-017-1139-9, 2017.

[14] S. Uçar, E. Uçar, N. Özdemir, Z. Hammouch, "Mathematical Analysis and Numerical Simulation for a Smoking Model With Atangana-Baleanu Derivative," Chaos, Solitons & Fractals, vol. 118, pp. 300-306, DOI: 10.1016/j.chaos.2018.12.003, 2019.

[15] P. Veeresha, D. G. Prakasha, H. M. Baskonus, "Solving Smoking Epidemic Model of Fractional Order Using a Modified Homotopy Analysis Transform Method," Mathematical Sciences, vol. 13 no. 2, pp. 115-128, DOI: 10.1007/s40096-019-0284-6, 2019.

[16] H. Günerhan, H. Rezazadeh, W. Adel, "Analytical Approximate Solution of Fractional Order Smoking Epidemic Model," Advances in Mechanical Engineering, vol. 14 no. 9,DOI: 10.1177/16878132221123888, 2022.

[17] Q. Haidong, M. U. Rahman, M. Arfan, "Fractional Model of Smoking with Relapse and Harmonic Mean Type Incidence Rate Under Caputo Operator," Journal of Applied Mathematics and Computing, vol. 69 no. 1, pp. 403-420, DOI: 10.1007/s12190-022-01747-6, 2023.

[18] N. Ilmayasinta, H. Purnawan, "Optimal Control in a Mathematical Model of Smoking," Journal of Mathematical and Fundamental Sciences, vol. 53 no. 3, pp. 380-394, DOI: 10.5614/j.math.fund.sci.2021.53.3.4, 2021.

[19] O. Khyar, J. Danane, K. Allali, "Mathematical Analysis and Optimal Control of Giving up the Smoking Model," International Journal of Differential Equations, vol. 2021,DOI: 10.1155/2021/8673020, 2021.

[20] H. Khan, J. Alzabut, A. Shah, S. Etemad, S. Rezapour, C. Park, "A Study on the Fractal-Fractional Tobacco Smoking Model," AIMS Mathematics, vol. 7 no. 8, pp. 13887-13909, DOI: 10.3934/math.2022767, 2022.

[21] P. Liu, T. Munir, T. Cui, A. Din, P. Wu, "Mathematical Assessment of the Dynamics of the Tobacco Smoking Model: An Application of Fractional Theory," AIMS Mathematics, vol. 7 no. 4, pp. 7143-7165, DOI: 10.3934/math.2022398, 2022.

[22] H. Günerhan, M. K. A. Kaabar, E. Çelik, "Novel Analytical and Approximate-Analytical Methods for Solving the Nonlinear Fractional Smoking Mathematical Model," Sigma: Linguistique Anglaise, Linguistique Generale, vol. 41 no. 2, pp. 331-343, 2023.

[23] I. R. Sofia, S. R. Bandekar, M. Ghosh, "Mathematical Modeling of Smoking Dynamics in Society with Impact of Media Information and Awareness," Results in Control and Optimization, vol. 11,DOI: 10.1016/j.rico.2023.100233, 2023.

[24] K. S. Nisar, R. Tabassum, M. A. Z. Raja, M. Shoaib, "Advanced Bio-Inspired Computing Paradigm for Nonlinear Smoking Model," Alexandria Engineering Journal, vol. 76, pp. 411-427, DOI: 10.1016/j.aej.2023.06.032, 2023.

[25] A. Ullah, H. Sakidin, S. Gul, "Sensitivity Analysis-Based Control Strategies of a Mathematical Model for Reducing Marijuana Smoking," AIMS Bioengineering, vol. 10 no. 4, pp. 491-510, DOI: 10.3934/bioeng.2023028, 2023.

[26] A. N. Chatterjee, P. K. Roy, "Smoking Habit: a Bio Mathematical Study," Commun. Math. Biol. Neurosci., vol. 2020, 2020.

[27] S. S. Ray, A. Atangana, S. C. O. Noutchie, M. Kurulay, N. Bildik, A. Kilicman, "Fractional Calculus and its Applications in Applied Mathematics and Other Sciences," Mathematical Problems in Engineering, vol. 2014,DOI: 10.1155/2014/849395, 2014.

[28] O. A. Arqub, A. El-Ajou, "Solution of the Fractional Epidemic Model by Homotopy Analysis Method," Journal of King Saud University Science, vol. 25 no. 1, pp. 73-81, DOI: 10.1016/j.jksus.2012.01.003, 2013.

[29] M. Farman, M. U. Saleem, A. Ahmad, M. O. Ahmad, "Analysis and Numerical Solution of SEIR Epidemic Model of Measles With Non-Integer Time Fractional Derivatives by Using Laplace Adomian Decomposition Method," Ain Shams Engineering Journal, vol. 9 no. 4, pp. 3391-3397, DOI: 10.1016/j.asej.2017.11.010, 2018.

[30] M. Adel, "Fractional Order System of Integro-Differential Equations: Theoretical and Numerical Simulations Using Appell-Type Changhee Polynomials," The Journal of Mathematics and Computer Science, vol. 34 no. 01, pp. 74-84, DOI: 10.22436/jmcs.034.01.07, 2024.

[31] M. H. T. Alshbool, "Bernstein Polynomials Method for Solving Multi-Order Fractional Neutral Pantograph Equations with Error and Stability Analysis," Results in Applied Mathematics, vol. 22,DOI: 10.1016/j.rinam.2024.100451, 2024.

[32] M. H. T. Alshbool, W. Shatanawi, I. Hashim, M. Sarr, "Residual Correction Procedure with Bernstein Polynomials for Solving Important Systems of Ordinary Differential Equations," Computers, Materials & Continua, vol. 64 no. 1, pp. 63-80, DOI: 10.32604/cmc.2020.09431, 2020.

[33] M. H. Alshbool, O. Isik, I. Hashim, "Fractional Bernstein Series Solution of Fractional Diffusion Equations With Error Estimate," Axioms, vol. 10 no. 1,DOI: 10.3390/axioms10010006, 2021.

[34] M. H. T. Alshbool, M. Mohammad, O. Isik, I. Hashim, "Fractional Bernstein Operational Matrices for Solving Integro-Differential Equations Involved by Caputo Fractional Derivative," Results in Applied Mathematics, vol. 14 no. 6,DOI: 10.1016/j.rinam.2022.100258, 2022.

[35] W. Li, L. Bai, Y. Chen, S. D. Santos, B. Li, "Solution of Linear Fractional Partial Differential Equations Based on the Operator Matrix of Fractional Bernstein Polynomials and Error Correction," International Journal of Innovative Computing and Applications, vol. 14, pp. 211-226, 2018.

[36] M. H. T. Alshbool, A. S. Bataineh, I. Hashim, O. R. Isik, "Solution of Fractional-Order Differential Equations Based on the Operational Matrices of New Fractional Bernstein Functions," Journal of King Saud University Science, vol. 29 no. 1,DOI: 10.1016/j.jksus.2015.11.004, 2017.

[37] T. Hussain, A. U. Awan, K. A. Abro, "A Passive Verses Active Exposure of Mathematical Smoking Model: a Role for Optimal and Dynamical Control," Nonlinear Engineering, vol. 11 no. 1, pp. 507-521, DOI: 10.1515/nleng-2022-0214, 2022.

[38] R. Ullah, M. Khan, G. Zaman, "Dynamical Features of a Mathematical Model on Smoking," Journal of Applied Environmental and Biological Sciences, vol. 6 no. 1, pp. 92-96, 2016.

[39] S. K. Sharma, A. N. Chatterjee, F. A. Basir, "Hopf Bifurcation and Optimal Control of HCV/HIV Co-infection Dynamics within Human: A Theoretical Study," Results in Control and Optimization, vol. 11,DOI: 10.1016/j.rico.2023.100234, 2023.

[40] A. N. Chatterjee, S. K. Sharma, F. Al Basir, "A Compartmental Approach to Modeling the Measles Disease: A Fractional Order Optimal Control Model," Fractal and Fractional, vol. 8 no. 8,DOI: 10.3390/fractalfract8080446, 2024.

[41] P. van den Driessche, J. Watmough, "Reproduction Numbers and Sub-threshold Endemic Equilibria for Compartmental Models of Disease Transmission," Mathematical Biosciences, vol. 180 no. 1-2, pp. 29-48, DOI: 10.1016/s0025-5564(02)00108-6, 2002.

[42] D. Matignon, "Stability Results for Fractional Differential Equations with Applications to Control Processing," Computational Engineering in Systems Applications, vol. 2, pp. 963-968, 1996.

[43] M. H. T. Alshbool, I. Hashim, "Multistage Bernstein Polynomials for the Solutions of the Fractional Order Stiff Systems," Journal of King Saud University Science, vol. 28 no. 4, pp. 280-285, DOI: 10.1016/j.jksus.2015.06.001, 2016.