\documentclass[10pt,twoside,fleqn]{article} \usepackage[margin=2cm]{geometry} \usepackage{graphicx} \usepackage{multirow} \usepackage [displaymath, mathlines]{lineno} \usepackage{titlesec} \usepackage{color} \usepackage{hyperref} \hypersetup{ colorlinks, linkcolor=blue, urlcolor=blue, citecolor=blue } \titlelabel{\thetitle.\quad} %\linenumbers \usepackage[font=small,format=plain,labelfont=bf,up]{caption} \usepackage{amsmath} \setlength{\mathindent}{0cm} \renewcommand\linenumberfont{\normalfont\bfseries\small} \renewcommand{\baselinestretch}{1} \setcounter{page}{1} \pagestyle{myheadings} \thispagestyle{empty} \markboth{\footnotesize \emph{\emph{International Journal of Scientific World}}}{\footnotesize \emph{\emph{International Journal of Scientific World}}} \date{} \begin{document} {\renewcommand{\arraystretch}{0.65} \begin{table}[ht] \begin{tabular}{ll} \multirow{8}{*}{\includegraphics[width=1.8cm]{logo}}&\\&{\scriptsize\emph{\textbf{International Journal of Scientific World}, 3 (x) (2015) xxx-xxx}}\\ &{\scriptsize\emph{www.sciencepubco.com/index.php/IJSW}}\\ &{\scriptsize\emph{\copyright Science Publishing Corporation}}\\ &{\scriptsize\emph{doi: }}\\ &{\scriptsize\emph {Research Paper}}\\ &{\scriptsize\emph {\textbf}}} \end{tabular} \end{table}} %\renewcommand{\arraystretch}{1} \centerline{} \centerline{} \centerline{} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \centerline {\huge{\bf Optimal control strategies in square root }} \centerline{} \centerline{\huge{\bf dynamics of smoking model}} \centerline{} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% My definition \newcommand{\mvec}[1]{\mbox{\bfseries\itshape #1}} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \centerline{\bf {Anwar Zeb $^{1*}$, Fiza Bibi $^{1}$, Gul Zaman $^{2}$}} \centerline{} {\small \centerline{\emph{$^{1}$ Department of Mathematics, COMSATS Institute of Information, Technology, Abbottabad, Pakistan }} \centerline{\emph{$^{2}$ Department of Mathematics, University of Malakand, Chakdara, Dir(Lower), Khyber Pakhtunkhwa, Pakistan}} \centerline{\emph{*Corresponding author E-mail: anwar55.ciit@yahoo.com}}} \centerline{} \centerline{} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \newtheorem{Theorem}{Theorem}[section] \newtheorem{Definition}[Theorem]{Definition} \newtheorem{Corollary}[Theorem]{Corollary} \newtheorem{Lemma}[Theorem]{Lemma} \newtheorem{Proof}[Theorem]{Proof} \newtheorem{Example}[Theorem]{Example} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \noindent \hspace{-3 pt}{\scriptsize \textbf{ Copyright \copyright 2015 Anwar Zeb et. al. 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.}} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \smallskip \noindent \hrulefill %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \noindent \textbf{Abstract}\\ \centerline{} In a recent paper [Anwar Zeb, Gul Zaman, Shaher Momani, Square-root Dynamics of a Giving Up Smoking Model, Appl. Math. Model., 37 (2013) 5326-5334], the authors presented a new model of giving up smoking model. In this paper, we introduce three control variables in the form of education campaign, anti-smoking gum, and anti-nicotive drugs/medecine for the eradication of smoking in a community. Using the optimal control theory, the optimal levels of the three controls are characterized, and then the existence and uniqueness for the optimal control pair are established. In order to do this, we minimize the number of potential and occasional smokers and maximize the number of quit smokers. We use Pontryagin's maximum principle to characterize the optimal levels of the three controls. The resulting optimality system is solved numerically by Matlab. \centerline{} \noindent {\footnotesize \emph{\textbf{Keywords}}: \emph{Mathematical model; Square root dynamics; Non-standard; finite difference scheme; Numerical analysis; Optimal control.}} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \noindent\hrulefill %============================= \section{Introduction} %============================= Modelling is a science which needs creative ability linked to a deep knowledge of the whole variety of methods offered by applied mathematics. Indeed, the design of a model has to be precisely related to the methods to be used to deal with the mathematical problems generated by the application of the model. Smoking is known to be the biggest cause of both preventable and premature not only in the US but also worldwide. Smoking-related diseases are cause of over 440,000 deaths in the US and about 105,000 UK annually \cite{Cas}. The life expectancy of smoker is cut short by 10-12 years and more than half of all smokers die from smoking-related diseases. Comparative smoking facts show that the risk of heart attack is 70\% high among smokers than among non-smokers. The incidence of lung cancer is ten times greater in smokers than non-smokers and one out of ten people will die from this disease \cite{Cas, Sharomi,Ham}. Some 80\% of smokers will at one time be diagnosed with heart disease, emphysema or chronic bronchitis of diseases attributable to tobacco habit, among 29\% are from lung cancer and 24\% are caused by heart disease \cite{Sharomi}. As cigarette smoke contain over 4,000 chemical compounds and toxins, which causes the above harmful infection to human health. Several authors did a lot of work in order to understand the dynamics of smoking. Castillo-Garsow et al. \cite{Cas} proposed a simple mathematical model for giving up smoking. They consider a system with total constant population which is divided into three classes: potential smokers (P), smokers (S), and quit smokers (Q). Sharomi and Gumel \cite{Sharomi} developed the above model by introducing mild and chain classes. They presented the development and public health impact of smoking related illnesses. Zaman \cite{Zam3} extended the giving-up smoking further by taking into account the occasional smoker comportment and presented dynamical interaction in an integer order. Then Zeb et. al \cite{Zeb} introduced the square root dynamics in the giving up smoking model. The model presented by Zeb et. al \cite{Zeb} is \begin{equation}\label{eq1.1} \begin{array}{rcl} &&\frac{dP}{dt} =\lambda-\beta\sqrt{PL}-(d+\mu)P,\cr\cr &&\frac{dL}{dt} =\beta\sqrt{PL}-(\gamma+d+\mu)L,\cr\cr &&\frac{dS}{dt} =\gamma L-(\delta+d+\mu)S,\cr\cr &&\frac{dQ}{dt}=\delta S-(\mu+d)Q. \end{array} \end{equation} The total population is $N(t)=P(t)+L(t)+S(t)+Q(t)$ along with the conservative law% \begin{equation}\label{eq1.2} \frac{dN}{dt} =\lambda-(\mu+d)N, \end{equation} Here $\mu$ is natural death rate, $\gamma$ is recover rate from infection, $\beta$ is transmission coefficient, $\delta$ is quit rate of smoking, $d$ represent death rate for potential smokers, occasional smoker, smoker and quit smoker related to smoking disease. Beside this optimal control theory is another area of mathematics that is used extensively in controlling the spread of infectious diseases. It is a power full mathematical tool that can be used to make decisions involving complex biological situation. On the behalf of this to control the spread of smoking in the community, we will consider possible control variables to decrease the attitude towards smoking. The optimality is taken to be minimize the number of light and chain smokers and maximize the number of quit smokers in a community. First, we will show the existence of an optimal control for the control problem and then we will derive the optimality system. To determine the optimal strategy for my optimal problem we will use the Pontryagin Maximum Principle. By using this principle, we will derive the optimality system consisting of the state and adjoint equations and will solve numerically the system by using an iterative method. The paper is organized as follows. In Section 2, we analyze the existence and stability of equilibria and use Lyaponuv function theory to present the global stability of disease-free equilibrium and use Dulac Criteria for the global stability of endemic equilibrium. A control system for the optimality and its existence, and the optimal control are derived in Section 3. Parameters estimation and numerical results are discussed in Section 4: Finally, we give conclusion. \section{Analysis of optimal control strategy} In this section, we apply the optimal control strategy to control the spread of smoking in the community. It is a power full mathematical tool that can be used to make decisions involving complex biological situation. In order to control the spread of smoking in the community, we will consider three possible control variables to decrease the attitude towards smoking. Our control variables represent education campaign $u_{1}(t),$ anti-smoking gum $u_{2}(t)$ and anti-nicotive drug/medecine $u_{3}(t).$ The control variables satisfy conditions $u_{1}\in[0,0.9]$ and $u_{2},u_{3}\in[0,1]$with $u_{1}(t)\leq u_{2}(t)\leq u_{3}(t)$. We introduced the control variables in the system (\ref{eq1.1}) is given by \begin{equation}\label{eq2.1} \begin{array}{rcl} &&\frac{dP}{dt} =\lambda-\beta\sqrt{PL}-(d+\mu)P-(1-u_1)P,\cr\cr &&\frac{dL}{dt} =\beta\sqrt{PL}-(\gamma+d+\mu+u_2)L+pu_3S,\cr\cr &&\frac{dS}{dt} =\gamma L-(\delta+d+\mu)S-(p+q)u_3S,\cr\cr &&\frac{dQ}{dt} =\delta S-(\mu+d)Q+(1-u_1)P+u_2L+qu_3S, \end{array} \end{equation} where $p$ and $q$ show the probabilities such that $p,q\in[0,1]$ and $p+q\leq1$. In this paper, we assume that the control functions $u_{1}(t), u_{2}(t),$ and $u_{3}(t)$ are bounded and Lebesague integrable function. We consider the cost (objective) function as follows \begin{equation}\label{eq2.2} J[u(t)]= \int_{0}^{t_{f}}[A_{1}S(t)-A_{2}Q(t)+\frac{1}{2}[r_{1}u_{1}^{2}(t)+r_{2}u_{2}^{2}(t)+r_{3}u_{3}^{2}(t)] ]dt+A_{3}P(t_f)+A_{4}L(t_f). \end{equation} We will investigate the existence and uniqueness of optimal control for the proposed model. Our aim is to find control $u_1^*, u_2^* $ and $u_3^*$, such that \begin{equation}\label{eq2.3} J[u_1^*,u_2^*,u_3^*]=\underset{\Omega}{min}{[J(u_1,u_2,u_3)],} \end{equation} where $\Omega=\{(u_1,u_2,u_3)\in L^{1}(0,t_f)|a_i\leq u_i\leq b_i, i=1,2,3\}$, $a_i,b_i, i=1,2,3,$ are fixed non-negative constants and $(u_1,u_2,u_3)\in L^{1}(0,t_f)$. Now we will prove the existence of the optimal control for the system (\ref{eq2.1}) and then derive the optimality system. To do this we construct the Hamiltonian function $H$ with respect to control variables is given by \begin{equation}\label{eq2.4} H={A_{1}S(t)-A_{2}Q(t)+\frac{1}{2}[r_{1}u_{1}^{2}(t)+r_{2}u_{2}^{2}(t)+r_{3}u_{3}^{2}(t)] +A_{3}P(t_f)+A_{4}L(t_f)+\sum_{i=1}^{4}\lambda_if_i,} \end{equation} where $f_i$ for $i=1, 2, 3, 4$ is the right side of the differential equations of the system (\ref{eq1.1}). Now we are able to state the results on the existence of the optimal control pair for the system (\ref{eq1.1}). \begin{Theorem} There exist control variables ${u}^{*}=(u_1^{*},u_2^{*},u_3^{*}\in\Omega)$ for the control problem (\ref{eq2.1})such that $$\underset{(u_1, u_2, u_3)\in\Omega}{min}J(u_1, u_2, u_3)=J(u_1^{*}, u_2^{*}, u_3^{*}).$$ \end{Theorem} \begin{Proof} For the proof of this theorem see the results \cite{Zam3}. \end{Proof} Now use the Pontryagin Maximum Principle to obtain the necessary conditions for optimal controls is given by \begin{equation}\label{eq2.5} \begin{array}{rcl} &&\frac{dX}{dt}=\frac{\partial H(t, X,u,\lambda)}{\partial \lambda},\cr\cr &&\frac{\partial H(t,X,u,\lambda)}{\partial u}=0,\cr\cr &&\lambda^{'}(t)=-\frac{H(t, X, u,\lambda)}{\partial X}. \end{array} \end{equation} By using these necessary conditions, we will find the optimal control system, the adjoint system and the characterization of the control variables. \begin{Theorem} Let $P^{*},L^{*},S^{*},Q^{*}$ be optimal state solutions with associated optimal control variables $(u_1^{*},u_2^{*},u_3^{*})$ for the optimal control problem (\ref{eq2.1}). Then there exist adjoint variables $\lambda_{i}$, for $i=1,2,3,4,$ satisfying \begin{equation}\label{eq2.6} \begin{array}{rcl} &&\lambda^{'}_1(t)=A_{2}(1-u_{1})\lambda_4+A_3\lambda_1(\frac{\beta}{2}\sqrt{\frac{L}{P}}+d+\mu+(1-u_1))-A_4\lambda_2\frac{\beta}{2}\sqrt{\frac{L}{P}},\cr\cr &&\lambda_{2}^{'}(t)=-A_1\lambda_3\gamma+A_2\lambda_4u_2+A_3\lambda_1\frac{\beta}{2}\sqrt{\frac{P}{L}}-A_4\lambda_2(\frac{\beta}{2}\sqrt{\frac{P}{L}} -(\gamma+d+\mu+u_2)),\cr\cr &&\lambda_3^{'}(t)=-\frac{H}{S}=A_1\lambda_3(\delta+d+\mu+(p+q)u_3)+A_2\lambda_4(\delta+qu_3S)-A_4\lambda_2pu_3,\cr\cr &&\lambda_{4}^{'}(t)=-\frac{H}{Q}=-A_2\lambda_4(\mu+d), \end{array} \end{equation} with transversality conditions $\lambda_i(t_f)=0, i=1,2,3,4.$ We also obtain the optimal controls $(u_1^{*},u_2^{*},u_3^{*})$ as \begin{equation}\label{eq2.7} \begin{array}{rcl} &&u_1^{*}=min{\{max{\{0,\frac{P(\lambda_4-\lambda_1)}{r_1}\}},u_{1max}}\},\cr\cr &&u_2^{*}=min{\{max{\{0,\frac{L(\lambda_2-\lambda_4)}{r_2}\}},u_{2max}}\},\cr\cr &&u_3^{*}=min{\{max{\{0,\frac{pS(\lambda_3-\lambda_2)+qS(\lambda_3-\lambda_4)}{r_3}\}},u_{3max}}\}. \end{array} \end{equation}\end{Theorem} \begin{Proof} To determine the adjoint equations and the transversality conditions, we use the Hamiltonian (\ref{eq2.4}). From setting $P(t)=P^{*}(t), L(t)=L^{*}(t), S(t)=S^{*}(t), Q(t)=Q^{*}(t)$ and differentiating the Hamiltonian (\ref{eq2.4}) with respect to $P(t), L(T),S(t),Q(t)$, respectively, we obtain equation (\ref{eq2.6}). By solving equations $\frac{H}{u_1}=0$, $\frac{H}{u_2}=0$ and $\frac{H}{u_3}=0$ on the interior of the control set and using the optimality conditions and the property of the control space $U$, we can derive (\ref{eq2.7}). \end{Proof} Therefore, taking the state system to gather with the adjoint system, the optimal control, and the transversality conditions, we have the following optimality system: \begin{equation}\label{eq2.8} \begin{array}{rcl} &&\frac{dP}{dt}=\lambda-\beta\sqrt{PL}-(d+\mu)P-(1-u_{1}^{*})P,\cr\cr &&\frac{dL}{dt}=\beta\sqrt{PL}-(\gamma+d+\mu+u_{2}^{*})L+pu_3^{*}S,\cr\cr &&\frac{dS}{dt}=\gamma L-(\delta+d+\mu)S-(p+q)u_3^{*}S,\cr\cr &&\frac{dQ}{dt}=\delta S-(\mu+d)Q+(1-u_{1}^{*})P+u_{2}^{*}L+qu_3^{*}S, \end{array} \end{equation} \begin{equation}\label{eq2.9} \begin{array}{rcl} &&\lambda^{'}_1=A_{2}(1-u_{1}^{*})\lambda_4+A_3\lambda_1(\frac{\beta}{2}\sqrt{\frac{L}{P}}+d+\mu+(1-u_1^{*})) -A_4\lambda_2\frac{\beta}{2}\sqrt{\frac{L}{P}},\cr\cr &&\lambda_{2}^{'}=-A_1\lambda_3\gamma+A_2\lambda_4u_2^{*}+A_3\lambda_1\frac{\beta}{2}\sqrt{\frac{P}{L}}-A_4\lambda_2(\frac{\beta}{2}\sqrt{\frac{P}{L}} -(\gamma+d+\mu+u_2^{*})),\cr\cr &&\lambda_3^{'}=A_1\lambda_3(\delta+d+\mu+(p+q)u_3^{*})+A_2\lambda_4(\delta+qu_3^{*}S)-A_4\lambda_2pu_3^{*},\cr\cr &&\lambda_{4}^{'}=-A_2\lambda_4(\mu+d), \end{array} \end{equation} $\lambda_{i}(t_{f}=0) , i=1,2,3,4,$ and $P(0)=P_{0},L(0)=L_{0}, S(0)=S_{0}, Q(0)=Q_{0}.$ \section{Numerical illustration} \subsection{The improved GSS1 method} The resolution of the optimality system is created improving the Gauss-Seidel-like implicit finite-difference method developed by Gumel et al. \cite{1}, presented in \cite{2} and denoted GSS1 method. It consists on discretizing the interval $[t_0; t_{end}]$ the points $t_k=kh + t_0~~~k = 0,1,..., n,$ where $h$ is the time step. Next, we define the state and adjoint variables $P(t),L(t),S(t),Q(t),$ $\lambda_{i},$ and the control $u_{j}(t),$ in terms of nodal points $P^k,L^k,S^k,Q^k,$ $\lambda_{i}^k,$ and $u_{j}^k,$ with $P^0,L^0,S^0,Q^0,$ $\lambda_{i}^0,$ and $u_{j}^0,$ as the state and adjoint variables and the controls at initial time $t_0$. $P^n,L^n,S^n,Q^n,$ $\lambda_{i}^n,$ and $u_{j}^n,$ as the state and adjoint variables and the controls at final time $t_{end}$ for $~i=1,2,3,4~~and~~j=1,2,3,$. As it is well known, the approximation of the time derivative by its first-order forward difference is given, for the first state variable $P(t)$, by$$\frac{dP(t)}{dt}=\lim_{h\rightarrow0}\frac{P(t+h)-P(t)}{h}$$ We use GSS1 to adapt it to our case as following: we visit the variables one by one by blocking all other value to the most recently calculated \begin{equation}\label {eq1} \frac{P^{k+1}-P^{k}}{h}=\lambda-\beta\sqrt{P^{k+1}L^k}-(d+\mu)P^{k+1}-(1-u_1^{k})P^{k+1}, \end{equation} \begin{equation}\label {eq2} \frac{L^{k+1}-L^{k}}{h}=\beta\sqrt{P^{k+1}L^{k+1}}-(\gamma+d+\mu+u_2)L^{k+1}+pu_3^{k}S^{k}, \end{equation} \begin{equation}\label {eq3} \frac{S^{k+1}-S^{k}}{h}=\gamma L^{k+1}-(\delta+d+\mu)S^{k+1}-(p+q)u_3^{k}S^{k+1}, \end{equation} \begin{equation}\label {eq4} \frac{Q^{k+1}-Q^{k}}{h}=\delta S^{k+1}-(\mu+d)Q^{k+1}+(1-u_1^{k})P^{k+1}+u_2^{k}L^{k+1}+qu_3^{k}S^{k+1}, \end{equation} By applying an analogous technology, we approximate the time derivative of the adjoint variables by their first-order backward-difference and we use the appropriated scheme as follows: \begin{equation*}\label{eq3.2} \begin{array}{rcl} &&\frac{\lambda_1^{n-k}-\lambda_1^{n-k-1}}{h}=A_{2}(1-u_{1}^{k})\lambda_4^{n-k}+A_3\lambda_1^{n-k-1}(\frac{\beta}{2}\sqrt{\frac{L^{k+1}}{P^{k+1}}}+d+\mu+(1-u_1^{k})) -A_4\lambda_2^{n-k}\frac{\beta}{2}\sqrt{\frac{L^{k+1}}{P^{k+1}}},\cr\cr &&\frac{\lambda_2^{n-k}-\lambda_2^{n-k-1}}{h}=-A_1\lambda_3^{n-k}\gamma+A_2\lambda_4^{n-k}u_2^{k}+A_3\lambda_1^{n-k}\frac{\beta}{2}\sqrt{\frac{P^{k+1}}{L^{k+1}}}-A_4\lambda_2^{n-k-1}(\frac{\beta}{2}\sqrt{\frac{P^{k+1}}{L^{k+1}}} -(\gamma+d+\mu+u_2^{k})),\cr\cr &&\frac{\lambda_3^{n-k}-\lambda_3^{n-k-1}}{h}=A_1\lambda_3^{n-k-1}(\delta+d+\mu+(p+q)u_3^{k})+A_2\lambda_4^{n-k}(\delta+qu_3^{k}S)-A_4\lambda_2^{n-k}pu_3^{k},\cr\cr &&\frac{\lambda_4^{n-k}-\lambda_4^{n-k-1}}{h}=-A_2\lambda_4^{n-k-1}(\mu+d). \end{array} \end{equation*} Hence, we can establish an algorithm to solve the optimality system and then to compute the optimal control. \subsection{Algorithm} Step 1: $$P(0)\leftarrow P_0,~L(0)\leftarrow L_0, ~ S(0)\leftarrow S_0,~Q(0)\leftarrow Q_0,$$ $$\lambda_1(t_{end})\leftarrow 0,~\lambda_2(t_{end})\leftarrow 0,~\lambda_3(t_{end})\leftarrow 0,~\lambda_4(t_{end})\leftarrow 0 $$ Step 2: For $k = 1,2,...,n$ do: Taking the transformation of variables% \[ v^{k+1}=\sqrt{P^{k+1}}, \] in equation (\ref {eq1}), we obtain quadratic equation for $v^{k+1},$% \[ (1+h(\mu+d+1-u_1^{k})){v^{2(k+1)}}+(\beta h\sqrt{L^{k}})v^{k+1}-(P^{k}+h \lambda)=0. \] Since our goal is to calculate, $P^{k+1}$ from knowledge of $(\lambda ,\mu,\beta,P^{k},L^{k})$ only the $P^{k+1}$ variable is required under the transformation equation. Solution of the above quadratic equation is given by% \begin{equation*} v^{k+1}=\left[ \frac{1}{2(1+h(\mu+d+1-u_1^{k}))}\right] \left[ -(\beta h \sqrt{L^{k}})+\sqrt{(\beta h\sqrt{L^{k}})^{2}+4(1+h(\mu+d+1-u_1^{k}))(P^{k}% +h\lambda)}\right] \end{equation*} Similarly, the remaining equations of system (\ref{eq1}-\ref{eq4}) can be solved {for variable} at $(k+1)th$ time step:% \begin{align*} L^{k+1} & =\left[ \frac{1}{(1+h(\gamma+\mu+d))}\right] \left( (\beta h\sqrt{L^{k}})v^{k+1}+L^{k}+ph(u_3S)^{k}\right) \\ S^{k+1} & =\left[ \frac{1}{(1+h(\delta+\mu+d+(p+q)u_3^{k}))}\right] \left( \gamma h L^{k+1}+S^{k}\right) \\ Q^{k+1} & =\left[ \frac{1}{(1+h(\mu+d))}\right] \left(h\{(\delta+qu_3^{k})S^{k+1}+(1-u_1^{k})P^{k+1}+u_2^{k}L^{k+1}\}+Q^{k}\right)\\ \end{align*} \begin{equation*}\label{eq3.2} \begin{array}{rcl} &&\lambda_1^{n-k-1}=\frac{\lambda_1^{n-k}+h(A_4\lambda_2^{n-k}\frac{\beta}{2}\sqrt{\frac{L^{k+1}}{P^{k+1}}}-A_{2}(1-u_{1}^{k})\lambda_4^{n-k})}{1+h(A_3(\frac{\beta}{2}\sqrt{\frac{L^{k+1}}{P^{k+1}}}+d+\mu+(1-u_1^{k})))} ,\cr\cr &&\lambda_2^{n-k-1}=\frac{\lambda_2^{n-k}+h(A_1\lambda_3^{n-k}\gamma-A_2\lambda_4^{n-k}u_2^{k}-A_3\lambda_1^{n-k}\frac{\beta}{2}\sqrt{\frac{P^{k+1}}{L^{k+1}}})}{1-h(A_4(\frac{\beta}{2}\sqrt{\frac{P^{k+1}}{L^{k+1}}}-(\gamma+d+\mu+u_2^{k})))} ,\cr\cr &&\lambda_3^{n-k-1}=\frac{\lambda_3^{n-k}+h(A_4\lambda_2^{n-k}pu_3^{k}-A_2\lambda_4^{n-k}(\delta+qu_3^{k}S))}{1+h(A_1(\delta+d+\mu+(p+q)u_3^{k}))},\cr\cr &&\lambda_4^{n-k-1}=\frac{\lambda_4^{n-k}}{1-h(A_2(\mu+d))}, \end{array} \end{equation*} $$\theta_1^{k+1}=\frac{P^{k+1}(\lambda_4^{n-k-1}-\lambda_1^{n-k-1})}{r_1},$$ $$\theta_2^{k+1}=\frac{L^{k+1}(\lambda_2^{n-k-1}-\lambda_4^{n-k-1})}{r_2},$$ $$\theta_3^{k+1}=\frac{S^{k+1}\{p(\lambda_3^{n-k-1}-\lambda_2^{n-k-1})+q(\lambda_3^{n-k-1}-\lambda_4^{n-k-1})\}}{r_3}$$ $$u_1^{k+1}=min{\{max{\{0,\theta_1^{k+1}\}},u_{1max}}\},$$ $$u_2^{k+1}=min{\{max{\{0,\theta_2^{k+1}\}},u_{2max}}\},$$ $$u_3^{k+1}=min{\{max{\{0,\theta_3^{k+1}\}},u_{3max}}\}$$ end for Step 3: For $k = 0,1,...,n$, write: $P^*(t_k) = P_k,~~L^*(t_k) = L_k, ~~S^*(t_k) = S_k, ~~Q^*(t_k) = Q_k,~~u_i^*(t_k) = u_k,~i=1,2,3.$ end for \begin{figure}[hbtp] \centering \includegraphics[height=53mm]{cr1.jpg} \includegraphics[height=53mm]{cr2.jpg} \includegraphics[height=53mm]{cr3.jpg} \includegraphics[height=53mm]{cr4.jpg} \caption{The plot represents the Potential smokers, Occasional smokers, Chain smokers and Quit smokers for both control and without control in time $\alpha=0.3,0.5,0.7,1$.} \label{fig1234} \end{figure} \newpage \section{Conclusion} In a recent paper [Anwar Zeb, Gul Zaman, Shaher Momani, Square-root Dynamics of a Giving Up Smoking Model, Appl. Math. Model., 37 (2013) 5326-5334], the authors presented a new model of giving up smoking model. In this paper, we introduce three control variables in the form of education campaign, anti-smoking gum, and anti-nicotive drugs/medecine for the eradication of smoking in a community. Using the optimal control theory, the optimal levels of the three controls are characterized, and then the existence and uniqueness for the optimal control pair are established. In order to do this, we minimize the number of potential and occasional smokers and maximize the number of quit smokers. We use Pontryagin's maximum principle to characterize the optimal levels of the three controls. The resulting optimality system is solved numerically by Matlab. %========================== \section*{Acknowledgements} %========================== {We wish to thank the anonymous referees for their thorough reading and constructive comments.} \begin{thebibliography}{99} {\small \bibitem {Cas}C. Castillo-Garsow, G. Jordan-Salivia and A. Rodriguez Herrera ``Mathematical Models for Dynamics of Tobacco Use, Recovery and Relapse" \textit{Technical Report Series BU-1505-M, Cornell Uneversity}, \textbf{(2000)}. \bibitem {Sharomi}O. Sharomi and A.B. Gumel ``Curtailing smoking dynamics: A mathematical modeling approach" \textit{Applied Mathematics and Computation,} 195 \textbf{(2008)} 475-499. \bibitem {Ham}O.K. Ham ``Stages and Processes of Smoking Cessation Among Adolescents" \textit{ West J. Nurs. Res.,} 29 \textbf{(2007)} 301-315. \bibitem {Zam3}G. Zaman ``Qualitative behavior of giving up smoking models" \textit{Bulletin of the Malaysian Mathematical Sciences Society,} 34 \textbf{(2011)} 403-415. \bibitem{Zeb} A. Zeb, G. Zaman and S. Momani``Square-root dynamics of a giving up smoking model" \textit{Applied Mathematical Modelling,} 37 \textbf{(2013)} 5326?5334 \bibitem {Rad}A.G. Radwan, K. Moaddy and S. Momani ``Stability and non-standard finite difference method of the generalized Chua's circuit" \textit{Computer and Mathematics with Applications,} \textbf{62} (3) (2011) 961-970. \bibitem {Kirschner}D. Kirschner, S. Lenhart, S. Serbin, \textit{Optimal control of chemotherapy of HIV,} J. Math. Biol., \textbf{35} (1997) 775-792. \bibitem {Zaman1}G. Zaman, Y.H. Kang, I.H. Jung,\textit{Optimal treatment of an SIR epidemic model with time delay,} Bio., \textbf{98} (1) (2009) 43-50. \bibitem {Kamien}M.I. Kamien, N.L. Schwartz, \textit{Dynamics Optimization,} The Calculus of Variations and Optimal Control in Economics and Management, 1991. \bibitem {Zaman2}G. Zaman, H. Jung, \textit{Optimal vaccination and treatment in the SIR epidemic model,} Proc. KSIAM, \textbf{3} (2007) 31-33. \bibitem {Riewe}F. Riewe, \textit{Nonconservative Lagrangian and Hamiltonian mechanics,} Phys. Rev. E., \textbf{53} (1996) 1890-1899. \bibitem {1} A.B. Gumel, P.N. Shivakumar, B.M. Sahai,\textit{A mathematical model for the dynamics of HIV-1 during the typical course of infection}, Nonlinear Anal., Theory Methods Appl., Ser. A, Theory Methods, \textbf{47} (3) (2001) 1773?1783. \bibitem {2} J. Karrakchou, M. Rachik, S. Gourari, \textit{Optimal control and infectiology: Application to an HIV/AIDS model,} Appl. Math. Comput., \textbf{177} (2006) 807?818. } \end{thebibliography} \end{document}