\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}
