On Non-Symmetric Fractal-Fractional Modeling for Ice Smoking: Mathematical Analysis of Solutions

Abstract

Drugs have always been one of the most important concerns of families and government officials at all times, and they have caused irreparable damage to the health of young people. Given the importance of this great challenge, this article discusses a non-symmetric fractal-fractional order ice-smoking mathematical model for the existence results, numerical results, and stability analysis. For the existence of the solution of the given ice-smoking model, successive iterative sequences are defined. The uniqueness of the solution Hyers–Ulam (HU) stability is established with the help of the existing definitions and theorems in functional analysis. By the utilization of two-step Lagrange polynomials, we provide numerical solutions and provide a comparative numerical analysis for different values of the fractional order and fractal order. The numerical simulations show the applicability of the scheme and future prediction and the effects of fractal-fractional orders simultaneously.

1. Introduction

Arab traders discovered and recognized the poppy plant (Opium) for the first time in Southwestern China and India in the 7th century. It was used as a medicine for various diseases. The cultivation of this plant was rare, but it became prevalent in the 15th century when European traders boosted their trades to these regions, see [1]. The European merchants were exchanging gold and silver with the Chinese and easterners for spices and silk. In [2], Cady says that in the 17th century, the Chinese people started smoking a mixture of tobacco and opium. This provided an opportunity for the Europeans to become rich and recoup their money from the Chinese by supplying opium from colonial India in the late 18th century. According to [1], the British started to supply 15 million opium addicts in China. The cultivation of the poppy crop spread in India and even to Myanmar and Laos, bringing profit to the Europeans. Millions of Asians became addicted to opium, and up to 1970s, all of the opium produced in this region was consumed by opium-addicted Asians.

In this era, the narcotics syndicate arrived, and they started the refinement of opium into heroin. The Sicilian Mafia introduced heroin into Europe and propagated it there. Legally produced opium from Turkey was taken to France and Marseille for refinement into heroin and transferred to the United States and Western Europe. It is a confirmed reality that drugs have ruined the lives of the millions people throughout the world. Germany, Japan, Pakistan, and the United States are the countries where the proportion of drug addicts is alarming. Morphine, churse, heroin, ice (methamphetamine), cannabis, and cocaine are the very common drugs. Among these drugs, ice claims the most addicts worldwide. According to a paper from Dr. Nowshad Khan and Shah Fahad [3], for the first time in world history, ice was used as drug in Japan in 1919; it was prepared in laboratories in 1960, although, as Philip Jenkins in [4] states, it did not become the main concern of the media and politicians until 1989 and 1990. Philip further explains that ice (smokable crystal methamphetamine) was a threat to American society. It was used by the army persons in World War II for prolonged duties. Ice stimulates one’s hormones and speeds up one’s activities by up to one thousand times. In the beginning, it provides extreme pleasure, but with the passage of time, it destroys cells in the body. Nowshad and Fahad state that doctors have reported that ice increases one’s heart rate and blood pressure, and it gradually makes humans lazy and weak.

To study the effects of such dangerous drugs in more detail, many different articles have been published in which researchers conduct various types of scientific research based on various mathematical data and models. Some of them can be mentioned, including [5,6,7]. In 2013, Zeb et al. [8] designed a square-root-structured model of giving up smoking and analyzed it via the finite difference method. In the same year, Huo et al. [9] investigated the effect of relapse on the cessation of smoking in the context of a mathematical model. Recently, in 2018, Zeb et al. published two papers on the dynamics of cigarette smoking [10,11]. Similar to them, one can mention the other works done by other researchers such as [12,13,14,15,16]. In the meantime, we can even see the high efficiency of mathematical models in studying the properties and dynamics of various processes and diseases, including the mathematical models of thermostat control [17], anthrax in animals [18], mosaic disease [19], Q fever [20], memristor-based circuits [21], Mump virus [22], the hybrid system of p-Laplacial operators [23], and other biologic and engineering models [24,25,26,27,28,29,30,31,32,33].

Our main purpose in this paper is the analysis of a mathematical system of the ice-smoking model with the help of the fractal-fractional operators. This specifies the motivation of our research; the application of the fractal-fractional operators for this model specifies the novelty of our research because we extend the previous results to a generalized system of fractal-fractional IVPs, and their simulations of the solutions of the mentioned system give more accurate results.

In this manuscript, the desired model of the fractal-fractional system is described in Section 2, and some preliminaries are stated in Section 3. Via the successive iterative sequences and limit points, we investigate the existence criterion in Section 4, and further, we verify its uniqueness in Section 5. By defining HU-stable solutions, we analyze this qualitative property for the given fractal-fractional ice-smoking model in Section 6. Furthermore, the two-step Lagrange polynomials help us to derive a numerical solution, and then we simulate the comparative results for different values of parameters of these fractal-fractional operators in Section 7. Lastly, we complete the research via the conclusion in Section 8.

2. Description of the Ice-Smoking Model

Recently, in 2021, Zeb and Alzahrani [34] studied the ice-smoking model based on the finite scheme along with the linear differential equation approach, which takes the form:

$$\left\{\begin{array}{c}\dot{S}\left(t\right)=\lambda -\beta SC-\mu S,\hfill \\ \dot{C}\left(t\right)=\beta SC-(\gamma +\mu +{\rho }_1)C+aCR,\hfill \\ \dot{R}\left(t\right)=\gamma C-(\mu +{\rho }_2)R-aCR,\hfill \\ \dot{Q}\left(t\right)={\rho }_{1}C+{\rho }_{2}R-\mu Q.\hfill \end{array}\right.$$

We intend to discuss the same model by the non-symmetric fractal-fractional order structure with existence, uniqueness, stability, and numerical simulations. We proceed as follows:

$$\left\{\begin{array}{c}{}_0^{FFM}{D}_t^{{\nu }_1,{\sigma }^{*}}S\left(t\right)=\lambda -\beta SC-\mu S,\hfill \\ {}_0^{FFM}{D}_t^{\nu ,{\sigma }^{*}}C\left(t\right)=\beta SC-(\gamma +\mu +{\rho }_1)C+aCR,\hfill \\ {}_0^{FFM}{D}_t^{{\nu }_1,{\sigma }^{*}}R\left(t\right)=\gamma C-(\mu +{\rho }_2)R-aCR,\hfill \\ {}_0^{FFM}{D}_t^{{\nu }_1,{\rho }^{*}}Q\left(t\right)={\rho }_{1}C+{\rho }_{2}R-\mu Q,\hfill \end{array}\right.$$

(1)

with initial conditions

$$S\left(0\right)=S^0,C\left(0\right)=C^0,R\left(0\right)=R^0,Q\left(0\right)=Q^0,$$

with $S^0,C^0,R^0,Q^0\ge 0$. Here, ${}_0^{FFM}{D}_t^{{\nu }_1,{\sigma }^{*}}$ stands for the fractal-fractional derivative in the Atangana–Baleanu sense for ${\nu }_1,{\sigma }^{*}\in (0,1]$. In the model (1), the whole population is grouped into 4 classes. They are $S\left(t\right)$, which illustrates potential smokers and potential non-smokers. The mentioned compartment or class is increased under the rate $\lambda$, denoting the recruitment rate.

The compartment or category C(t) specifies the chain smokers. It is increased when potential smokers begin to smoke under an incidence rate or under the contact rate among the potential smokers $\beta S\left(t\right)C\left(t\right)$ and $aC\left(t\right)R\left(t\right)$ of the relapse smokers who have returned to smoking. Note that some other individuals will leave this category under the rates $\gamma C\left(t\right),{\rho }_{1}C\left(t\right)$ and $\mu C\left(t\right)$. Additionally, the category $C\left(t\right)$ is increased under the rate $\gamma C\left(t\right)$ in which $\gamma$ is the rate at which regular smokers transfer to the reversion category. Moreover, $C\left(t\right)$ is decreased under the rates $\mu R\left(t\right)$, ${\rho }_{2}R\left(t\right)$. Accordingly, $Q\left(t\right)$ represents the permanent quitters.

In this model, $\lambda$ is a parameter denoting the birth or migration rate to the host population. The rate $\mu$ shows the natural death in all four categories. The rate ${\rho }_1$ shows the number of regular smokers who quit, and the rate $rh{o}_2$ shows the number of reversion persons who quit. $\beta$ is the incidence rate in relation to the susceptible persons to the class of the regular smokers. Additionally, a is the relapse rate, and the rate $\alpha$ shows the number of individuals who relapsed and became regular smokers.

3. Preliminaries

In this article, we presume the Banach’s space $\left\{\varphi \right(t)\in C([0,1]\to \mathbb{R}\left)\right\}$ with the norm $\parallel \varphi \parallel ={max}_{t\in [0,1]}\left|\varphi \left(t\right)\right|$. Here, we highlight the basic notions of the fractal-fractional calculus. This derivative is based on the Mittag–Leffler kernel, which is a non-singular kernel. This derivative is recently recommended by several authors for the dynamical problems and has many advantages too. Like other fractional derivatives, it possesses the potential to have the classical results as well new solutions for the fractal and fractional orders. Regarding the fractal-fractional operators and their basic notions, see [35,36].

Definition 1

([35,36]). Consider $\varphi \in C\left(\right(a,b),\mathbb{R})$ which is fractal differentiable on $(a,b)$ of order $0<{\varrho }^{*}\le 1$. The fractal-fractional derivation operator for ϕ in the Atangana–Baleanu settings of order $0<{\kappa }_1\le 1$, with the generalized kernel of the Mittag–Leffler type, is introduced as

$${}_0^{FFM}{D}_t^{{\kappa }_1,{\varrho }^{*}}\varphi \left(t\right)=\frac{AB\left({\kappa }_1\right)}{1-{\kappa }_1}\frac{d}{d{t}^{{\varrho }^{*}}}{\int }_0^t\varphi \left(s\right)E_{{\kappa }_1}\left[-\frac{{\kappa }_1}{1-{\kappa }_1}{(t-s)}^{{\kappa }_1}\right]ds,$$

where $AB\left({\kappa }_1\right)=1-{\kappa }_1+{\displaystyle \frac{{\kappa }_1}{\Gamma {\kappa }_1}}$, and

$$\frac{d\varphi \left(s\right)}{d{s}^{{\varrho }^{*}}}=\underset{t\to s}{lim}\frac{\varphi \left(t\right)-\varphi \left(s\right)}{t^{{\varrho }^{*}}-s^{{\varrho }^{*}}}.$$

Definition 2

([35,36]). Let ϕ be the same function considered above. Then, the fractal-fractional integration operator in the Atangana–Baleanu settings for ϕ of order $0<{\kappa }_1\le 1$ with the kernel of Mittag–Leffler type is given by

$${}_0^{FFM}{I}_t^{{\kappa }_1,{\varrho }^{*}}\varphi \left(t\right)=\frac{{\kappa }_1{\varrho }^{*}}{AB\left({\kappa }_1\right)\Gamma {\kappa }_1}{\int }_0^t{s}^{{\varrho }^{*}-1}\varphi \left(s\right){(t-s)}^{{\kappa }_1-1}ds+\frac{{\varrho }^{*}(1-{\kappa }_1)t^{{\varrho }^{*}-1}}{AB\left({\kappa }_1\right)}\varphi \left(t\right),$$

where $AB\left({\kappa }_1\right)=1-{\kappa }_1+{\displaystyle \frac{{\kappa }_1}{\Gamma {\kappa }_1}}.$

By making use of a successive iterative technique, we shall accomplish the proof for the existence criterion to the fractal-fractional model (1). For this, we operate the integral given in Definition 2 from [36] on the proposed model (1), and we have

$$\begin{array}{cc}\hfill S\left(t\right)-S\left(0\right)& =\frac{{\nu }_1{\sigma }^{*}}{AB\left({\nu }_1\right)\Gamma {\nu }_1}{\int }_0^t{s}^{{\sigma }^{*}-1}{(t-s)}^{{\nu }_1-1}[\lambda -\beta SC-\mu S]ds\hfill \\ \hfill & +\frac{{\sigma }^{*}(1-{\nu }_1)t^{{\sigma }^{*}-1}}{AB\left({\nu }_1\right)}[\lambda -\beta SC-\mu S],\hfill \\ \hfill C\left(t\right)-C\left(0\right)& =\frac{{\nu }_1{\sigma }^{*}}{AB\left({\nu }_1\right)\Gamma {\nu }_1}{\int }_0^t{s}^{{\sigma }^{*}-1}{(t-s)}^{{\nu }_1-1}[\beta SC-(\gamma +\mu +{\rho }_1)C+aCR]ds\hfill \\ \hfill & +\frac{{\sigma }^{*}(1-{\nu }_1)t^{{\sigma }^{*}-1}}{AB\left({\nu }_1\right)}[\beta SC-(\gamma +\mu +{\rho }_1)C+aCR],\hfill \\ \hfill R\left(t\right)-R\left(0\right)& =\frac{{\nu }_1{\sigma }^{*}}{AB\left({\nu }_1\right)\Gamma {\nu }_1}{\int }_0^t{s}^{{\sigma }^{*}-1}{(t-s)}^{{\nu }_1-1}[\gamma C-(\mu +{\rho }_2)R-aCR]ds\hfill \\ \hfill & +\frac{{\sigma }^{*}(1-{\nu }_1)t^{{\sigma }^{*}-1}}{AB\left({\nu }_1\right)}[\gamma C-(\mu +{\rho }_2)R-aCR],\hfill \\ \hfill Q\left(t\right)-Q\left(0\right)& =\frac{{\nu }_1{\sigma }^{*}}{AB\left({\nu }_1\right)\Gamma {\nu }_1}{\int }_0^t{s}^{{\sigma }^{*}-1}{(t-s)}^{{\nu }_1-1}[{\rho }_{1}C+{\rho }_{2}R-\mu Q]ds\hfill \\ \hfill & +\frac{{\sigma }^{*}(1-{\nu }_1)t^{{\sigma }^{*}-1}}{AB\left({\nu }_1\right)}[{\rho }_{1}C+{\rho }_{2}R-\mu Q].\hfill \end{array}$$

(2)

Regard the functions ${\mathcal{V}}_i$ for $i=1,2,\dots ,4$ or $i\in {\mathbb{N}}_1^4$, given below:

$$\left\{\begin{array}{c}{\mathcal{V}}_1(t,S)=\lambda -\beta SC-\mu S\hfill \\ {\mathcal{V}}_2(t,C)=\beta SC-(\gamma +\mu +{\rho }_1)C+aCR,\hfill \\ {\mathcal{V}}_3(t,R)=\gamma C-(\mu +{\rho }_2)R-aCR,\hfill \\ {\mathcal{V}}_4(t,Q)={\rho }_{1}C+{\rho }_{2}R-\mu Q.\hfill \end{array}\right.$$

(3)

4. Existence Criteria

In the present part, for establishing the desired theorem on the existence property, we provide an assumption:

$\left(H^{*}\right)$ All $S\left(t\right)$, $S^{*}\left(t\right)$, $C\left(t\right)$, $C^{*}\left(t\right)$, $R\left(t\right)$, $R^{*}\left(t\right)$, $Q\left(t\right)$, $Q^{*}\left(t\right)$ $\in L[0,1]$ are continuous so that $\parallel S\parallel \le a_1,\parallel C\parallel \le a_2,\parallel R\parallel \le a_3,\parallel Q\parallel \le a_4$ for some positive constants $a_1,a_2,a_3,a_4>0$. Furthermore, we define the following constants: ${\varphi }_1=\beta a_2+\mu$, ${\varphi }_2=\beta a_1+(\gamma +\mu +{\rho }_1)+a{a}_3$, ${\varphi }_3=\mu +{\rho }_2+a{a}_2$, ${\varphi }_4=\mu$.

Theorem 1.

The Lipschitz condition is valid for the kernels ${\mathcal{V}}_i$ for $i\in {\mathbb{N}}_1^4$ if the assumption ($H^{*}$) fulfills and ${\varphi }_i<1$ for $i\in {\mathbb{N}}_1^4$.

Proof. 

We first check ${\mathcal{V}}_1(t,S)$ for the Lipschitz property. For this, we are helped by $\left(H^{*}\right)$ and (3) and obtain

$$\begin{array}{ccc}\hfill \parallel {\mathcal{V}}_1\left(S\right)-{\mathcal{V}}_1\left(S^{*}\right)\parallel & =& \parallel (\lambda -\beta SC-\mu S)-(\lambda -\beta S^{*}C-\mu S^{*})\parallel \hfill \\ & \le & (\beta \parallel S-S^{*}\parallel \parallel C\parallel +\mu \parallel S-S^{*}\parallel \hfill \\ & \le & (\beta a_2+\mu )\parallel S-S^{*}\parallel \hfill \\ & =& {\varphi }_1\parallel S-S^{*}\parallel ,\hfill \end{array}$$

(4)

where ${\varphi }_1=\beta a_2+\mu$. Hence, ${\mathcal{V}}_1$ satisfies the Lipchitz condition with the Lipchitz-constant ${\varphi }_1$. Similarly, for ${\mathcal{V}}_2(t,C)$, we have

$$\begin{array}{ccc}\hfill \parallel {\mathcal{V}}_2\left(C\right)-{\mathcal{V}}_2\left(C^{*}\right)\parallel & =& \parallel (\beta SC-(\gamma +\mu +{\rho }_1)C+aCR)-(\beta S{C}^{*}-(\gamma +\mu +{\rho }_1)C^{*}+a{C}^{*}R)\parallel \hfill \\ & \le & \left[\beta a_1+(\gamma +\mu +{\rho }_1)+a{a}_3\right]\parallel C-C^{*}\parallel \hfill \\ & =& {\varphi }_2\parallel C-C^{*}\parallel ,\hfill \end{array}$$

(5)

where ${\varphi }_2=\beta a_1+(\gamma +\mu +{\rho }_1)+a{a}_3$. Hence, ${\mathcal{V}}_2$ fulfills the Lipschitz property with constant ${\varphi }_2$. For ${\mathcal{V}}_3(t,R)$, we also have

$$\begin{array}{ccc}\hfill \parallel {\mathcal{V}}_3\left(R\right)-{\mathcal{V}}_3\left(R^{*}\right)\parallel & =& \parallel (\gamma C-(\mu +{\rho }_2)R-aCR)-(\gamma C-(\mu +{\rho }_2)R^{*}-aC{R}^{*})\parallel \hfill \\ & \le & \left[\parallel (\mu +{\rho }_2)(R^{*}-R)+aC(R^{*}-R)\parallel \right]\hfill \\ & \le & (\mu +{\rho }_2+a\parallel C\parallel )\parallel R-R^{*}\parallel \hfill \\ & =& (\mu +{\rho }_2+a{a}_2)\parallel R-R^{*}\parallel \hfill \\ & =& {\varphi }_3\parallel R-R^{*}\parallel ,\hfill \end{array}$$

(6)

where ${\varphi }_3=\mu +{\rho }_2+a{a}_2$. This implies that ${\mathcal{V}}_3$ satisfies the Lipchitz condition via constant ${\varphi }_3$. Now, for ${\mathcal{V}}_4(t,Q)$, we have

$$\begin{array}{ccc}\hfill \parallel {\mathcal{V}}_4\left(Q\right)-{\mathcal{V}}_4\left(Q^{*}\right)\parallel & =& \parallel ({\rho }_{1}C+{\rho }_{2}R-\mu Q)-({\rho }_{1}C+{\rho }_{2}R-\mu Q^{*})\parallel \hfill \\ & \le & \parallel \mu (Q^{*}-Q)\parallel \hfill \\ & =& \mu \parallel Q-Q^{*}\parallel \hfill \\ & =& {\varphi }_4\parallel Q-Q^{*}\parallel ,\hfill \end{array}$$

(7)

where ${\varphi }_4=\mu$. Thus, ${\mathcal{V}}_4$ also fulfills the Lipschitz property with constant ${\varphi }_4$. Thus, from (4)–(7), we have that ${\mathcal{V}}_i$ for $i=1,2,\dots ,4$, satisfy the Lipschitz property, and the result is accomplished. □

Let us assume:

$$\left\{\begin{array}{c}{\displaystyle S\left(t\right)-S\left(0\right)=\frac{{\nu }_1{\sigma }^{*}}{AB\left({\nu }_1\right)\Gamma {\nu }_1}{\int }_0^t{(t-s)}^{{\nu }_1-1}{s}^{{\sigma }^{*}-1}{\mathcal{V}}_1(s,S\left(s\right))ds+\frac{{\sigma }^{*}(1-{\nu }_1)}{AB\left({\nu }_1\right)}{t}^{{\sigma }^{*}-1}{\mathcal{V}}_1(t,S\left(t\right)),}\hfill \\ {\displaystyle C\left(t\right)-C\left(0\right)=\frac{{\nu }_1{\sigma }^{*}}{AB\left({\nu }_1\right)\Gamma {\nu }_1}{\int }_0^t{(t-s)}^{{\nu }_1-1}{s}^{{\sigma }^{*}-1}{\mathcal{V}}_2(s,C\left(s\right))ds+\frac{{\sigma }^{*}(1-{\nu }_1)}{AB\left({\nu }_1\right)}{t}^{{\sigma }^0*-1}{\mathcal{V}}_2(t,C\left(t\right)),}\hfill \\ {\displaystyle R\left(t\right)-R\left(0\right)=\frac{{\nu }_1{\sigma }^{*}}{AB\left({\nu }_1\right)\Gamma {\nu }_1}{\int }_0^t{(t-s)}^{{\nu }_1-1}{s}^{{\sigma }^{*}-1}{\mathcal{V}}_3(s,R\left(s\right))ds+\frac{{\sigma }^{*}(1-{\nu }_1)}{AB\left({\nu }_1\right)}{t}^{{\sigma }^{*}-1}{\mathcal{V}}_3(t,R\left(t\right)),}\hfill \\ {\displaystyle Q\left(t\right)-R\left(0\right)=\frac{{\nu }_1{\sigma }^{*}}{AB\left({\nu }_1\right)\Gamma {\nu }_1}{\int }_0^t{(t-s)}^{{\nu }_1-1}{s}^{{\sigma }^{*}-1}{\mathcal{V}}_4(s,Q\left(s\right))ds+\frac{{\sigma }^{*}(1-{\nu }_1)}{AB\left({\nu }_1\right)}{t}^{{\sigma }^{*}-1}{\mathcal{V}}_4(t,Q\left(t\right)).}\hfill \end{array}\right.$$

Now, we define the following recursive formulas for the model (1):

$$\begin{array}{ccc}\hfill S_n\left(t\right)-S\left(0\right)& =& \frac{{\nu }_1{\sigma }^{*}}{AB\left({\nu }_1\right)\Gamma {\nu }_1}{\int }_0^t{(t-s)}^{{\nu }_1-1}{s}^{{\sigma }^{*}-1}{\mathcal{V}}_1(s,S_{n-1}\left(s\right))ds\hfill \\ & +& \frac{{\sigma }^{*}(1-{\nu }_1)}{AB\left({\nu }_1\right)}{t}^{{\sigma }^{*}-1}{\mathcal{V}}_1(t,S_{n-1}\left(t\right)),\hfill \end{array}$$

$$\begin{array}{ccc}\hfill C_n\left(t\right)-C\left(0\right)& =& \frac{{\nu }_1{\sigma }^{*}}{AB\left({\nu }_1\right)\Gamma {\nu }_1}{\int }_0^t{(t-s)}^{{\nu }_1-1}{s}^{{\sigma }^{*}-1}{\mathcal{V}}_2(s,C_{n-1}\left(s\right))ds\hfill \\ & +& \frac{{\sigma }^{*}(1-{\nu }_1)}{AB\left({\nu }_1\right)}{t}^{{\sigma }^{*}-1}{\mathcal{V}}_2(t,C_{n-1}\left(t\right)),\hfill \end{array}$$

$$\begin{array}{ccc}\hfill R_n\left(t\right)-R\left(0\right)& =& \frac{{\nu }_1{\sigma }^{*}}{AB\left({\nu }_1\right)\Gamma {\nu }_1}{\int }_0^t{(t-s)}^{{\nu }_1-1}{s}^{{\sigma }^{*}-1}{\mathcal{V}}_3(s,R_{n-1}\left(s\right))ds\hfill \\ & +& \frac{{\sigma }^{*}(1-{\nu }_1)}{AB\left({\nu }_1\right)}{t}^{{\sigma }^{*}-1}{\mathcal{V}}_3(t,R_{n-1}\left(t\right)).\hfill \end{array}$$

$$\begin{array}{ccc}\hfill Q_n\left(t\right)-Q\left(0\right)& =& \frac{{\nu }_1{\sigma }^{*}}{AB\left({\nu }_1\right)\Gamma {\nu }_1}{\int }_0^t{(t-s)}^{{\nu }_1-1}{s}^{{\sigma }^{*}-1}{\mathcal{V}}_4(s,Q_{n-1})\left(s\right))ds\hfill \\ & +& \frac{{\sigma }^{*}(1-{\nu }_1)}{AB\left({\nu }_1\right)}{t}^{{\sigma }^{*}-1}{\mathcal{V}}_4(t,Q_{n-1})\left(t\right)),\hfill \end{array}$$

Theorem 2.

Under assumption $\left(H^{*}\right)$, the fractal-fractional ice-smoking model (1) has a solution if we have

$$\Delta =max[{\varphi }_1,{\varphi }_2,{\varphi }_3,{\varphi }_4]<1.$$

(8)

Proof. 

We define four functions as follows:

$$\left\{\begin{array}{c}\mathcal{U}{1}_n\left(t\right)=S_{n+1}\left(t\right)-S\left(t\right),\hfill \\ \mathcal{U}{2}_n\left(t\right)=C_{n+1}\left(t\right)-C\left(t\right),\hfill \\ \mathcal{U}{3}_n\left(t\right)=R_{n+1}\left(t\right)-R\left(t\right),\hfill \\ \mathcal{U}{4}_n\left(t\right)=Q_{n+1}\left(t\right)-Q\left(t\right).\hfill \end{array}\right.$$

(9)

Then, we find that

$$\begin{array}{ccc}\hfill \parallel \mathcal{U}{1}_n\left(t\right)\parallel & =& \frac{{\nu }_1{\sigma }^{*}}{AB\left({\nu }_1\right)\Gamma {\nu }_1}{\int }_0^t{(t-s)}^{{\nu }_1-1}{s}^{{\sigma }^{*}-1}\parallel [{\mathcal{V}}_1(s,S_n\left(s\right))-{\mathcal{V}}_1(s,S\left(s\right))]\parallel ds\hfill \\ & +& \frac{{\sigma }^{*}(1-{\nu }_1)}{AB\left({\nu }_1\right)}{t}^{{\sigma }^{*}-1}\parallel [{\mathcal{V}}_1(t,S_n\left(t\right))-{\mathcal{V}}_1(t,S\left(t\right))]\parallel \hfill \\ & \le & \frac{{\nu }_1{\sigma }^{*}}{AB\left({\nu }_1\right)\Gamma {\nu }_1}{\int }_0^t{(t-s)}^{{\nu }_1-1}{s}^{{\sigma }^{*}-1}{\varphi }_1\parallel S_n-S\parallel ds\hfill \\ & +& \frac{{\sigma }^{*}(1-{\sigma }_1)}{AB\left({\nu }_1\right)}{t}^{{\sigma }^{*}-1}{\varphi }_1\parallel S_n-S\parallel \hfill \\ & \le & [\frac{{\nu }_1{\sigma }^{*}\Gamma \left({\sigma }^{*}\right)}{AB\left({\nu }_1\right)\Gamma ({\nu }_1+{\sigma }^{*})}+\frac{{\sigma }^{*}(1-{\nu }_1)}{AB\left({\nu }_1\right)}]{\varphi }_1\parallel S_n-S\parallel \hfill \\ & \le & [\frac{{\nu }_1{\sigma }^{*}\Gamma \left({\sigma }^{*}\right)}{AB\left({\nu }_1\right)\Gamma ({\nu }_1+{\sigma }^{*})}+\frac{{\sigma }^{*}(1-{\nu }_1)}{AB\left({\nu }_1\right)}]{\varphi }_1\parallel S_n-S\parallel \hfill \\ & \le & {[\frac{{\nu }_1\Gamma ({\sigma }^{*}+1)}{AB\left({\nu }_1\right)\Gamma ({\nu }_1+{\sigma }^{*})}+\frac{{\sigma }^{*}(1-{\nu }_1)}{AB\left({\nu }_1\right)}]}^n{\varphi }_1^n\parallel S_1-S\parallel .\hfill \end{array}$$

(10)

In which, for ${\varphi }_1<1$ and as $n\to \infty$, we have $S_n\to S$. So, $\mathcal{U}{1}_n\to 0$ as $n\to \infty$. Similarly,

$$\begin{array}{ccc}\hfill \parallel \mathcal{U}{2}_n\left(t\right)\parallel & \le & {[\frac{{\nu }_1\Gamma ({\sigma }^{*}+1)}{AB\left({\nu }_1\right)\Gamma ({\nu }_1+{\sigma }^{*})}+\frac{{\sigma }^{*}(1-{\nu }_1)}{AB\left({\nu }_1\right)}]}^n{\varphi }_2^n\parallel C_1-C\parallel ,\hfill \end{array}$$

(11)

$$\begin{array}{ccc}\hfill \parallel \mathcal{U}{3}_n\left(t\right)\parallel & \le & {[\frac{{\nu }_1\Gamma ({\sigma }^{*}+1)}{AB\left({\nu }_1\right)\Gamma ({\nu }_1+{\sigma }^{*})}+\frac{{\sigma }^{*}(1-{\nu }_1)}{AB\left({\nu }_1\right)}]}^n{\varphi }_3^n\parallel R_1-R\parallel ,\hfill \end{array}$$

(12)

$$\begin{array}{ccc}\hfill \parallel \mathcal{U}{4}_n\left(t\right)\parallel & \le & {[\frac{{\nu }_1\Gamma ({\sigma }^{*}+1)}{AB\left({\nu }_1\right)\Gamma ({\nu }_1+{\sigma }^{*})}+\frac{{\sigma }^{*}(1-{\nu }_1)}{AB\left({\nu }_1\right)}]}^n{\varphi }_4^n\parallel Q_1-Q\parallel .\hfill \end{array}$$

(13)

By (10)–(13), when $n\to \infty$, then $\mathcal{U}{i}_n\left(t\right)\to 0,i\in {\mathbb{N}}_2^4$, for ${\varphi }_i<1,(i=2,\dots ,4)$. Ultimately, the ice-smoking system (1) has a solution. □

5. Unique Solution

For our suggested ice-smoking model (1), we follow the analysis of the uniqueness property.

Theorem 3.

The fractal-fractional ice-smoking model (1) possesses one solution exactly if $\left(H^{*}\right)$ is satisfied and the following holds:

$$\begin{array}{c}\hfill \left[\frac{{\nu }_1\Gamma ({\sigma }^{*}+1)}{AB\left({\nu }_1\right)\Gamma ({\nu }_1+{\sigma }^{*})}+\frac{{\sigma }^{*}(1-{\nu }_1)}{AB\left({\nu }_1\right)}\right]{\varphi }_i\le 1,i\in {\mathbb{N}}_1^4.\end{array}$$

(14)

Proof. 

We assume that the conclusion of the theorem is not valid. In other words, another solution exists for the supposed ice-smoking (I M) model (1) in the fractal-fractional settings. Hence, $S^{*}\left(t\right),C^{*}\left(t\right),R^{*}\left(t\right),Q^{*}\left(t\right)$ is another solution with $S^{*}\left(0\right)=S^0,C^{*}\left(0\right)=C^0,R^{*}\left(0\right)=R^0,Q^{*}\left(0\right)=Q^0$ s.t.

$$\begin{array}{ccc}\hfill S^{*}\left(t\right)-S^{*}\left(0\right)& =& \frac{{\nu }_1{\sigma }^{*}}{AB\left({\nu }_1\right)\Gamma {\nu }_1}{\int }_0^t{(t-s)}^{{\nu }_1-1}{s}^{{\sigma }^{*}-1}{\mathcal{V}}_1(s,S^{*}\left(s\right))ds,\hfill \\ & +& \frac{{\sigma }^{*}(1-{\nu }_1)}{AB\left({\nu }_1\right)}{t}^{{\sigma }^{*}-1}{\mathcal{V}}_1(t,S^{*}\left(t\right)),\hfill \end{array}$$

(15)

and similarly,

$$\begin{array}{ccc}\hfill C^{*}\left(t\right)-C^{*}\left(0\right)& =& \frac{{\nu }_1{\sigma }^{*}}{AB\left({\nu }_1\right)\Gamma {\nu }_1}{\int }_0^t{(t-s)}^{{\nu }_1-1}{s}^{{\sigma }^{*}-1}{\mathcal{V}}_2(s,C^{*}\left(s\right))ds\hfill \\ & +& \frac{{\sigma }^{*}(1-{\nu }_1)}{AB\left({\nu }_1\right)}{t}^{{\sigma }^{*}-1}{\mathcal{V}}_2(t,C^{*}\left(t\right)),\hfill \end{array}$$

(16)

$$\begin{array}{ccc}\hfill R^{*}\left(t\right)-R^{*}\left(0\right)& =& \frac{{\nu }_1{\sigma }^{*}}{AB\left({\nu }_1\right)\Gamma {\nu }_1}{\int }_0^t{(t-s)}^{{\nu }_1-1}{s}^{{\sigma }^{*}-1}{\mathcal{V}}_3(s,R^{*}\left(s\right))ds,\hfill \\ & +& \frac{{\sigma }^{*}(1-{\nu }_1)}{AB\left({\nu }_1\right)}{t}^{{\sigma }^{*}-1}{\mathcal{V}}_3(t,R^{*}\left(t\right)),\hfill \end{array}$$

(17)

$$\begin{array}{ccc}\hfill Q^{*}\left(t\right)-Q^{*}\left(0\right)& =& \frac{{\nu }_1{\sigma }^{*}}{AB\left({\nu }_1\right)\Gamma {\nu }_1}{\int }_0^t{(t-s)}^{{\nu }_1-1}{s}^{{\sigma }^{*}-1}{\mathcal{V}}_4(s,Q^{*}\left(s\right))ds,\hfill \\ & +& \frac{{\sigma }^{*}(1-{\nu }_1)}{AB\left({\nu }_1\right)}{t}^{{\sigma }^{*}-1}{\mathcal{V}}_4(t,Q^{*}\left(t\right)).\hfill \end{array}$$

(18)

Now, we write

$$\begin{array}{cc}\hfill \parallel S-S^{*}\parallel & =\frac{{\nu }_1{\sigma }^{*}}{AB\left({\nu }_1\right)\Gamma {\nu }_1}{\int }_0^t{(t-s)}^{{\nu }_1-1}{s}^{{\sigma }^{*}-1}\parallel {\mathcal{V}}_1-{\mathcal{V}}_1\left(S^{*}\right)\parallel \hfill \\ \hfill & +\frac{{\sigma }^{*}(1-{\nu }_1)}{AB\left({\nu }_1\right)}{t}^{{\sigma }^{*}-1}\parallel {\mathcal{V}}_1\left(S\right)-{\mathcal{V}}_1\left(S^{*}\right)\parallel \hfill \\ \hfill & \le \frac{{\nu }_1{\sigma }^{*}}{AB\left({\nu }_1\right)\Gamma {\nu }_1}{\int }_0^t{(t-s)}^{{\nu }_1-1}{s}^{{\sigma }^{*}-1}{\varphi }_1\parallel S-S^{*}\parallel \hfill \\ \hfill & +\frac{{\sigma }^{*}(1-{\nu }_1)}{AB\left({\nu }_1\right)}{t}^{{\sigma }^{*}-1}{\varphi }_1\parallel S-S^{*}\parallel \hfill \\ \hfill & \le \left[\frac{{\kappa }_1{\sigma }^{*}\Gamma \left({\sigma }^{*}\right)}{AB\left({\nu }_1\right)\Gamma ({\nu }_1+{\sigma }^{*})}+\frac{{\sigma }^{*}(1-{\nu }_1)}{AB\left({\nu }_1\right)}\right]{\varphi }_1\parallel S-S^{*}\parallel ,\hfill \end{array}$$

so

$$\left[1-\left[\frac{{\nu }_1{\sigma }^{*}\Gamma \left({\sigma }^{*}\right)}{AB\left({\nu }_1\right)\Gamma ({\nu }_1+{\sigma }^{*})}+\frac{{\sigma }^{*}(1-{\nu }_1)}{AB\left({\nu }_1\right)}\right]{\varphi }_1\right]\parallel S-S^{*}\parallel \le 0.$$

(19)

The above inequality (19) is true if $\parallel S-S^{*}\parallel =0$; accordingly, $S=S^{*}$. Similarly, from

$$\parallel C-C^{*}\parallel \le \left[\frac{{\nu }_1{\sigma }^{*}\Gamma \left({\sigma }^{*}\right)}{AB\left({\nu }_1\right)\Gamma ({\nu }_1+{\sigma }^{*})}+\frac{{\sigma }^{*}(1-{\nu }_1)}{AB\left({\nu }_1\right)}\right]{\varphi }_2\parallel C-C^{*}\parallel ,$$

we arrive at

$$\left[1-\left[\frac{{\nu }_1{\sigma }^{*}\Gamma \left({\sigma }^{*}\right)}{AB\left({\nu }_1\right)\Gamma ({\nu }_1+{\sigma }^{*})}+\frac{{\sigma }^{*}(1-{\nu }_1)}{AB\left({\nu }_1\right)}\right]{\varphi }_2\right]\parallel C-C^{*}\parallel \le 0.$$

This implies $\parallel C-C^{*}\parallel =0$ and $C=C^{*}$. Additionally,

$$\left[1-\left[\frac{{\nu }_1{\sigma }^{*}\Gamma \left({\sigma }^{*}\right)}{AB\left({\nu }_1\right)\Gamma ({\nu }_1+{\sigma }^{*})}+\frac{{\sigma }^{*}(1-{\nu }_1)}{AB\left({\nu }_1\right)}\right]{\varphi }_3\right]\parallel R-R^{*}\parallel \le 0.$$

This inequality is true, if $\parallel R-R^{*}\parallel =0$; accordingly, $R=R^{*}$. In similar manner, the inequality

$$\left[1-\left[\frac{{\nu }_1{\sigma }^{*}\Gamma \left({\sigma }^{*}\right)}{AB\left({\nu }_1\right)\Gamma ({\nu }_1+{\sigma }^{*})}+\frac{{\sigma }^{*}(1-{\nu }_1)}{AB\left({\nu }_1\right)}\right]{\varphi }_4\right]\parallel Q-Q^{*}\parallel \le 0.$$

is valid if $\parallel Q-Q^{*}\parallel =0$, which gives $Q=Q^{*}$. Therefore, the ice-smoking model (1) contains a unique solution. □

6. Hyers–Ulam Stability

The notion of Hyers–Ulam-Stability (HU-stability) is investigated for the solutions of the suggested ice-smoking model (1).

Definition 3.

The fractal-fractional ice-smoking system (1) is termed as HU-stable if $\exists {\eta }_i>0,i\in {\mathbb{N}}_1^4$ provided that $\forall {\varsigma }_i>0,i\in {\mathbb{N}}_1^4$ and for each $(S^{*},C^{*},R^{*},Q^{*})$ satisfying

$$\left\{\begin{array}{c}{|}_0^{FFM}{D}_t^{{\nu }_1,{\sigma }^{*}}{S}^{*}\left(t\right)-{\mathcal{V}}_1(t,S^{*})|\le {\xi }_1,\hfill \\ {|}_0^{FFM}{D}_t^{{\nu }_1,{\sigma }^{*}}{C}^{*}\left(t\right)-{\mathcal{V}}_1(t,C^{*})|\le {\xi }_2,\hfill \\ {|}_0^{FFM}{D}_t^{{\nu }_1,{\sigma }^{*}}{R}^{*}\left(t\right)-{\mathcal{V}}_1(t,R^{*})|\le {\xi }_3,\hfill \\ {|}_0^{FFM}{D}_t^{{\nu }_1,{\sigma }^{*}}{Q}^{*}\left(t\right)-{\mathcal{V}}_1(t,Q^{*})|\le {\xi }_4,\hfill \end{array}\right.$$

(20)

$(S,C,R,Q)$ exists, satisfying the ice-smoking system (1); further, we have

$$\left\{\begin{array}{c}\parallel S-S^{*}\parallel \le {\eta }_1{\xi }_1,\hfill \\ \parallel C-C^{*}\parallel \le {\eta }_2{\xi }_2,\hfill \\ \parallel R-R^{*}\parallel \le {\eta }_3{\xi }_3,\hfill \\ \parallel Q-Q^{*}\parallel \le {\eta }_4{\xi }_4.\hfill \end{array}\right.$$

where ${\mathcal{V}}_i,i\in {\mathbb{N}}_1^4$ are introduced in (3).

Remark 1.

Consider that the function $S^{*}$ is a solution of the first inequality (20) iff a continuous map $h_1$ exists (depending on $S^{*}$) so that (a) $|h_1\left(t\right)|<{\xi }_1$, and

$$\left(b\right){}_0^{FFM}{D}_t^{{\nu }_1,{\sigma }^{*}}{S}^{*}\left(t\right)={\mathcal{V}}_1(t,S^{*})+h_1\left(t\right).$$

Similarly, one can indicate such a definition for each of the solutions to the inequalities (20) by finding $h_i$ for $i\in {\mathbb{N}}_2^4$.

Theorem 4.

Let the hypothesis $\left(H^{*}\right)$ be true. Then, the fractal-fractional ice-smoking model (1) is HU-stable if

$$\left[\frac{{\nu }_1\Gamma ({\sigma }^{*}+1)}{AB\left({\nu }_1\right)\Gamma ({\nu }_1+{\sigma }^{*})}+\frac{{\sigma }^{*}(1-{\nu }_1)}{AB\left({\nu }_1\right)}\right]{\varphi }_i\le 1,i\in {\mathbb{N}}_1^4.$$

Proof. 

Let ${\xi }_1>0$ and the function $S^{*}$ be arbitrary so that

$${|}_0^{FFM}{D}_t^{{\nu }_1,{\sigma }^{*}}{S}^{*}\left(t\right)-{\mathcal{V}}_1(t,S^{*})|\le {\xi }_1.$$

In view of Remark 1, we have a function such as $h_1$ with $|h_1\left(t\right)|<{\xi }_1$, which satisfies

$${}_0^{FFM}{D}_t^{{\nu }_1,{\sigma }^{*}}{S}^{*}\left(t\right)={\mathcal{V}}_1(t,S^{*})+h_1\left(t\right).$$

Accordingly, we obtain

$$\begin{array}{cc}\hfill S^{*}\left(t\right)& =S^0+\frac{{\nu }_1{\sigma }^{*}}{AB\left({\nu }_1\right)\Gamma \left({\nu }_1\right)}{\int }_0^t{(t-s)}^{{\nu }_1-1}{s}^{{\sigma }^{*}-1}{\mathcal{V}}_1(s,S^{*}\left(s\right))ds+\frac{{\sigma }^{*}(1-{\nu }_1)}{AB\left({\nu }_1\right)}{t}^{{\sigma }^{*}-1}{\mathcal{V}}_1(t,S^{*}\left(t\right))\hfill \\ \hfill & +\frac{{\nu }_1{\sigma }^{*}}{AB\left({\nu }_1\right)\Gamma \left({\nu }_1\right)}{\int }_0^t{(t-s)}^{{\nu }_1-1}{s}^{{\sigma }^{*}-1}{h}_1\left(s\right)ds+\frac{{\sigma }^{*}(1-{\nu }_1)}{AB\left({\nu }_1\right)}{t}^{{\sigma }^{*}-1}{h}_1\left(t\right).\hfill \end{array}$$

Consider S as the unique solution of the fractal-fractional ice-smoking model (1). Then, it becomes

$$S\left(t\right)=S^0+\frac{{\nu }_1{\sigma }^{*}}{AB\left({\nu }_1\right)\Gamma \left({\nu }_1\right)}{\int }_0^t{(t-s)}^{{\nu }_1-1}{s}^{{\sigma }^{*}-1}{\mathcal{V}}_1(s,S\left(s\right))ds+\frac{{\sigma }^{*}(1-{\nu }_1)}{AB\left({\nu }_1\right)}{t}^{{\sigma }^{*}-1}{\mathcal{V}}_1(t,S\left(t\right)).$$

Hence,

$$\begin{array}{cc}\hfill |S^{*}\left(t\right)-S\left(t\right)|& \le \frac{{\nu }_1{\sigma }^{*}}{AB\left({\nu }_1\right)\Gamma \left({\nu }_1\right)}{\int }_0^t{(t-s)}^{{\nu }_1-1}{s}^{{\sigma }^{*}-1}|{\mathcal{V}}_1(s,S^{*}\left(s\right))-{\mathcal{W}}_1(s,S\left(s\right))|ds\hfill \\ \hfill & +\frac{{\sigma }^{*}(1-{\nu }_1)}{AB\left({\nu }_1\right)}{t}^{{\sigma }^{*}-1}|{\mathcal{V}}_1(t,S^{*}\left(t\right))-{\mathcal{V}}_1(t,S\left(t\right))|\hfill \\ \hfill & +\frac{{\nu }_1{\sigma }^{*}}{AB\left({\nu }_1\right)\Gamma \left({\nu }_1\right)}{\int }_0^t{(t-s)}^{{\nu }_1-1}{s}^{{\sigma }^{*}-1}\left|h_1\left(s\right)\right|ds\hfill \\ \hfill & +\frac{{\sigma }^{*}(1-{\nu }_1)}{AB\left({\nu }_1\right)}{t}^{{\sigma }^{*}-1}\left|h_1\left(t\right)\right|\hfill \\ \hfill & \le \left[\frac{{\nu }_1{\sigma }^{*}\Gamma \left({\sigma }^{*}\right)}{AB\left({\nu }_1\right)\Gamma ({\nu }_1+{\sigma }^{*})}+\frac{{\sigma }^{*}(1-{\nu }_1)}{AB\left({\nu }_1\right)}\right]{\varphi }_1|S^{*}\left(t\right)-S\left(t\right)|\hfill \\ \hfill & +\left[\frac{{\kappa }_1{\varrho }^{*}\Gamma \left({\varrho }^{*}\right)}{AB\left({\kappa }_1\right)\Gamma ({\kappa }_1+{\varrho }^{*})}+\frac{{\sigma }^{*}(1-{\nu }_1)}{AB\left({\nu }_1\right)}\right]{\xi }_1.\hfill \end{array}$$

In consequence,

$$\begin{array}{c}\hfill \parallel S^{*}-S\parallel \le {\displaystyle \frac{\left[\frac{{\nu }_1\Gamma ({\sigma }^{*}+1)}{AB\left({\nu }_1\right)\Gamma ({\nu }_1+{\sigma }^{*})}+{\displaystyle \frac{{\sigma }^{*}(1-{\nu }_1)}{AB\left({\nu }_1\right)}}\right]{\xi }_1}{1-\left[{\displaystyle \frac{{\nu }_1\Gamma ({\sigma }^{*}+1)}{AB\left({\nu }_1\right)\Gamma ({\nu }_1+{\sigma }^{*})}}+{\displaystyle \frac{{\sigma }^{*}(1-{\sigma }_1)}{AB\left({\nu }_1\right)}}\right]{\varphi }_1}}.\end{array}$$

(21)

If, we take

$${\eta }_1:=\frac{\left[{\displaystyle \frac{{\nu }_1\Gamma ({\sigma }^{*}+1)}{AB\left({\nu }_1\right)\Gamma ({\nu }_1+{\sigma }^{*})}}+{\displaystyle \frac{{\sigma }^{*}(1-{\nu }_1)}{AB\left({\nu }_1\right)}}\right]}{1-\left[{\displaystyle \frac{{\nu }_1\Gamma ({\sigma }^{*}+1)}{AB\left({\nu }_1\right)\Gamma ({\nu }_1+{\sigma }^{*})}}+{\displaystyle \frac{{\sigma }^{*}(1-{\nu }_1)}{AB\left({\nu }_1\right)}}\right]{\varphi }_1},$$

then $\parallel S^{*}-S\parallel \le {\eta }_1{\xi }_1$. Similarly, we have

$$\begin{array}{c}\hfill \parallel C^{*}-C\parallel \le {\eta }_2{\xi }_2,\parallel R^{*}-R\parallel \le {\eta }_3{\xi }_3,\parallel Q^{*}-Q\parallel \le {\eta }_4{\xi }_4.\end{array}$$

Thus, the fractal-fractional ice-smoking model (1) is HU-stable, which ends the argument. □

7. Numerical Algorithm

In this section, we describe the numerical scheme in relation to the fractal-fractional ice-smoking model (1). For this, we have received help from the technique regarding the two-step Lagrange polynomials. For the numerical scheme, consider the linear general differential equation ${}^{FFM}{}_0{D}_t^{{\nu }_1,{\sigma }^{*}}\vartheta \left(t\right)=\mathcal{V}(t,\vartheta \left(t\right))$, where $\vartheta \left(0\right)={\vartheta }_0$ is the initial value. The latter equation can be rewritten with respect to the Atangana–Baleanu fractla-fractional derivative as ${}^{AB}{}_0{D}_t^{{\nu }_1,{\sigma }^{*}}\vartheta \left(t\right)={\sigma }^{*}{t}^{{\sigma }^{*}-1}\mathcal{V}(t,\vartheta \left(t\right))=\mathcal{Y}(t,\vartheta \left(t\right)).$ With the help of the fractal-fractional integral operator that has a kernel of the generalized Mittag–Leffler type, we obtain

$$\vartheta \left(t\right)=\vartheta \left(0\right)+\frac{1-{\nu }_1}{AB\left({\nu }_1\right)}\mathcal{Y}(t,\vartheta \left(t\right))+\frac{{\nu }_1}{AB\left({\nu }_1\right)\Gamma {\nu }_1}{\int }_0^t{\varpi }^{{\sigma }^{*}-1}{(t-\varpi )}^{{\nu }_1-1}\mathcal{Y}(\varpi ,\vartheta \left(\varpi \right))d\varpi .$$

By replacing t with $t_{n+1}$, we have

$${\vartheta }^{n+1}=\vartheta \left(0\right)+\frac{1-{\nu }_1}{AB\left({\nu }_1\right)}\mathcal{Y}(t_n,\vartheta \left(t_n\right))+\frac{{\nu }_1}{AB\left({\nu }_1\right)\Gamma {\nu }_1}{\int }_0^{t_{n+1}}{\varpi }^{{\sigma }^{*}-1}{(t_{n+1}-\varpi )}^{{\nu }_1-1}\mathcal{Y}(\vartheta ,\vartheta \left(\varpi \right))d\varpi .$$

(22)

According to two-step Lagrange polynomials, we have

$$\begin{array}{cc}\hfill \mathcal{V}(x,\vartheta (t\left)\right)& =\frac{(x-t_{\ell -1})\mathcal{V}(t_{\ell },\vartheta \left(t_{\ell }\right))}{t_{\ell }-t_{\ell -1}}-\frac{(x-t_{\ell })\mathcal{V}(t_{\ell -1},\vartheta \left(t_{\ell -1}\right))}{t_{\ell }-t_{\ell -1}}\hfill \\ \hfill & =\frac{\mathcal{V}(t_{\ell },\vartheta \left(t_{\ell }\right)(x-t_{\ell -1})}{t_{\ell }-t_{\ell -1}}-\frac{\mathcal{V}(t_{\ell -1},\vartheta \left(t_{\ell -1}\right))(x-t_{\ell })}{t_{\ell }-t_{\ell -1}}\hfill \\ \hfill & =\frac{\mathcal{V}(t_{\ell },{\vartheta }_{\ell })(x-t_{\ell -1})}{h}-\frac{\mathcal{V}(t_{\ell -1},{\vartheta }_{\ell -1})(x-t_{\ell })}{h}.\hfill \end{array}$$

In this case, if we use the aforesaid Lagrange polynomial to (22), we obtain

$$\begin{array}{cc}\hfill {\vartheta }^{n+1}& =\vartheta \left(0\right)+\frac{1-{\nu }_1}{AB\left({\nu }_1\right)}\mathcal{V}(t_n,\vartheta \left(t_n\right))\hfill \\ \hfill & +\frac{{\nu }_1}{AB\left({\nu }_1\right)\Gamma {\nu }_1}\sum _{\ell =1}^n[\frac{\mathcal{V}(t_{\ell },\vartheta \left(t_{\ell }\right))}{h}{\int }_{t_{\ell }}^{t_{\ell +1}}(\varpi -t_{\ell -1}){(t_{n+1}-\varpi )}^{{\nu }_1-1}d\varpi \hfill \\ \hfill & -\frac{\mathcal{V}(t_{\ell -1},\vartheta \left(t_{\ell -1}\right))}{h}{\int }_{t_{\ell }}^{t_{n+1}}(\varpi -t_{\ell }){(t_{n+1}-\varpi )}^{{\nu }_1-1}d\varpi ].\hfill \end{array}$$

Further, we solve the above integral equation and obtain

$$\begin{array}{cc}\hfill {\vartheta }^{n+1}& ={\vartheta }_0+\frac{1-{\nu }_1}{AB\left({\nu }_1\right)}\mathcal{V}(t_n,\vartheta \left(t_n\right))\hfill \\ \hfill & +\frac{{\nu }_1{h}^{{\nu }_1}}{AB\left({\nu }_1\right)\Gamma ({\nu }_1+2)}\sum _{\ell =1}^n\left[\mathcal{V}(t_{\ell },\vartheta \left(t_{\ell }\right))\right({(n+1-\ell )}^{{\nu }_1}(n-\ell +2+{\nu }_1)\hfill \\ \hfill & -{(n-\ell )}^{{\nu }_1}(n-\ell +2+2{\nu }_1){\displaystyle )}\hfill \\ \hfill & -\mathcal{V}(t_{\ell -1},{\vartheta }_{\ell -1})\left({(n+1-\ell )}^{{\nu }_1+1}-(n-\ell +1+{\nu }_1){(n-\ell )}^{{\nu }_1}\right)].\hfill \end{array}$$

Inserting the value of $\mathcal{Y}(t,\vartheta (t\left)\right)$, it becomes

$$\begin{array}{cc}\hfill {\vartheta }^{n+1}& ={\vartheta }_0+{\sigma }^{*}{t}_n^{{\sigma }^{*}-1}\frac{1-{\nu }_1}{AB\left({\nu }_1\right)}\mathcal{V}(t_n,\vartheta \left(t_n\right))\hfill \\ \hfill & +\frac{{\sigma }^{*}{h}^{{\nu }_1}}{AB\left({\nu }_1\right)\Gamma ({\nu }_1+2)}\sum _{\ell =1}^n\left[t_{\ell }^{{\sigma }^{*}-1}\mathcal{V}(t_{\ell },\vartheta \left(t_{\ell }\right))\right({(n+1-\ell )}^{{\nu }_1}(n-\ell +2+{\nu }_1)\hfill \\ \hfill & -{(n-\ell )}^{{\nu }_1}(n-\ell +2+2{\nu }_1){\displaystyle )}\hfill \\ \hfill & -t_{\ell -1}^{{\sigma }^{*}-1}\mathcal{V}(t_{\ell -1},{\vartheta }_{\ell -1})\left({(n+1-\ell )}^{{\nu }_1+1}-(n-\ell +1+{\nu }_1){(n-\ell )}^{{\nu }_1}\right)].\hfill \end{array}$$

Thus, by assuming

$$\begin{array}{cc}\hfill {\Phi }_1(n,\ell )& :={(n+1-\ell )}^{{\nu }_1}(n-\ell +2+{\nu }_1)-{(n-\ell )}^{{\nu }_1}(n-\ell +2+2{\nu }_1),\hfill \\ \hfill {\Phi }_2(n,\ell )& :={(n+1-\ell )}^{{\nu }_1+1}-(n-\ell +1+{\nu }_1){(n-\ell )}^{{\nu }_1},\hfill \end{array}$$

the numerical scheme for the integral system (2) is obtained as

$$\begin{array}{cc}\hfill S\left(t_{n+1}\right)& =S\left(0\right)+{\sigma }^{*}{t}_n^{{\sigma }^{*}-1}\frac{1-{\nu }_1}{AB\left({\nu }_1\right)}{\mathcal{V}}_1(t_n,S\left(t_n\right))+\frac{{\sigma }^{*}{h}^{{\nu }_1}}{AB\left({\nu }_1\right)\Gamma ({\nu }_1+2)}\hfill \\ \hfill & \times \sum _{\ell =1}^n\left[t_{\ell }^{{\sigma }^{*}-1}{\mathcal{V}}_1(t_{\ell },S\left(t_{\ell }\right)){\Phi }_1(n,\ell )-t_{\ell -1}^{{\sigma }^{*}-1}{\mathcal{V}}_1(t_{\ell -1},S\left(t_{\ell -1}\right)){\Phi }_2(n,\ell )\right],\hfill \end{array}$$

$$\begin{array}{cc}\hfill C\left(t_{n+1}\right)& =C\left(0\right)+{\sigma }^{*}{t}_n^{{\sigma }^{*}-1}\frac{1-{\nu }_1}{AB\left({\nu }_1\right)}{\mathcal{V}}_2(t_n,C\left(t_n\right))+\frac{{\sigma }^{*}{h}^{{\nu }_1}}{AB\left({\nu }_1\right)\Gamma ({\nu }_1+2)}\hfill \\ \hfill & \times \sum _{\ell =1}^n\left[t_{\ell }^{{\sigma }^{*}-1}{\mathcal{V}}_2(t_{\ell },C\left(t_{\ell }\right)){\Phi }_1(n,\ell )-t_{\ell -1}^{{\sigma }^{*}-1}{\mathcal{V}}_2(t_{\ell -1},C\left(t_{\ell -1}\right)){\Phi }_2(n,\ell )\right],\hfill \end{array}$$

$$\begin{array}{cc}\hfill R\left(t_{n+1}\right)& =R\left(0\right)+{\sigma }^{*}{t}_n^{{\sigma }^{*}-1}\frac{1-{\nu }_1}{AB\left({\nu }_1\right)}{\mathcal{V}}_3(t_n,R\left(t_n\right))+\frac{{\sigma }^{*}{h}^{{\nu }_1}}{AB\left({\nu }_1\right)\Gamma ({\nu }_1+2)}\hfill \\ \hfill & \times \sum _{\ell =1}^n\left[t_{\ell }^{{\sigma }^{*}-1}{\mathcal{V}}_3(t_{\ell },R\left(t_{\ell }\right)){\Phi }_1(n,\ell )-t_{\ell -1}^{{\sigma }^{*}-1}{\mathcal{V}}_3(t_{\ell -1},R\left(t_{\ell -1}\right)){\Phi }_2(n,\ell )\right],\hfill \end{array}$$

$$\begin{array}{cc}\hfill Q\left(t_{n+1}\right)& =Q\left(0\right)+{\sigma }^{*}{t}_n^{{\sigma }^{*}-1}\frac{1-{\nu }_1}{AB\left({\nu }_1\right)}{\mathcal{V}}_4(t_n,Q\left(t_n\right))+\frac{{\sigma }^{*}{h}^{{\nu }_1}}{AB\left({\nu }_1\right)\Gamma ({\nu }_1+2)}\hfill \\ \hfill & \times \sum _{\ell =1}^n\left[t_{\ell }^{{\sigma }^{*}-1}{\mathcal{V}}_4(t_{\ell },R\left(t_{\ell }\right)){\Phi }_1(n,\ell )-t_{\ell -1}^{{\sigma }^{*}-1}{\mathcal{V}}_4(t_{\ell -1},Q\left(t_{\ell -1}\right)){\Phi }_2(n,\ell )\right].\hfill \end{array}$$

Computational Results

In this part, the numerical results about the fractal-fractional model (1) are provided with the help of our numerical scheme. The basic data were taken from the work given in [34]. We have: $\mu =0.0001,\beta =0.0005,\alpha =0.0002,{\rho }_1=0.004,{\rho }_2=0.002,\gamma =0.007,\lambda =0.9$, and the fractal fractional orders ${\nu }_1,{\sigma }^{*}=1.0,0.95,0.90,0.85$.

In these figures, we can see the effect of fractal-fractional orders and dimensions on the numerical solutions of the given model. In Figure 1, by increasing the fractal-fractional orders and tending to the integer order, the slope of the graphs decrease, and after 120 days, all of them become stable. In Figure 2, the opposite of this happens. By increasing the fractal-fractional orders and tending to the integer order, the slope of the graphs is positive and it increases, and after 120 days, all of them become stable.

In Figure 3, we see the same behavior in the diagrams as in Figure 1. Finally, Figure 4 depicts an increase in the number of the quitter smokers by increasing the time, and all of the diagrams experience a slight slope by tending ${\nu }_1,{\sigma }^{*}$ to the integer-order.

8. Conclusions

In this article, we have studied a fractal-fractional order ice-smoking mathematical model (1) for the qualitative analysis and computational aspect. The qualitative work is based on the fixed-point approach, while the numerical simulations are obtained with the help of the Lagrange’s interpolation polynomial. For the numerical simulations, we analyzed the fractal-fractional order ice-smoking mathematical model (1) for the orders ${\nu }_1={\sigma }^{*}=1.0,0.95,0.90,0.85$. For further new development/continuation of the study on the subject area, we can use other fractional operators that have exponential decay type kernels for some new ideas.