\documentclass[reqno]{amsart}
\usepackage{hyperref}
\usepackage{graphicx}

\AtBeginDocument{{\noindent\small
\emph{Electronic Journal of Differential Equations},
Vol. 2015 (2015), No. 277, pp. 1--15.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu}
\thanks{\copyright 2015 Texas State University.}
\vspace{8mm}}

\begin{document}
\title[\hfilneg EJDE-2015/277\hfil 
 Spatially-structured environmental-economic model]
{Dynamics and optimal control for a spatially-structured \\
environmental-economic model}

\author[D. La Torre, D. Liuzzi, T. Malik, O. Sharomi, R. Zaki 
\hfil EJDE-2015/277\hfilneg]
{Davide La Torre, Danilo Liuzzi, Tufail Malik, \\ 
Oluwaseun Sharomi, Rachad Zaki}

\address{Davide La Torre \newline
Department of Economics, Management and Quantitative Methods, 
University of Milan, Milan, Italy. \newline
Department of Applied Mathematics and Sciences, 
Khalifa University of Science, Technology and Research, 
Abu Dhabi, United Arab Emirates}
\email{davide.latorre@unimi.it, davide.latorre@kustar.ac.ae}


\address{Danilo Liuzzi \newline
Department of Economics, Management and Quantitative Methods,
 University of Milan, Milan, Italy. \newline
Department of Applied Mathematics and Sciences, 
Khalifa University of Science, Technology and Research,
 Abu Dhabi, United Arab Emirates}
\email{danilo.liuzzi@unimi.it, danilo.liuzzi@kustar.ac.ae}

\address{Tufail Malik \newline
Department of Applied Mathematics and Sciences,
Khalifa University of Science, Technology and Research\\
Abu Dhabi, United Arab Emirates}
\email{tufail.malik@kustar.ac.ae}

\address{Oluwaseun Sharomi (corresponding author) \newline
Department of Applied Mathematics and Sciences,
Khalifa University of Science, Technology and Research\\
Abu Dhabi, United Arab Emirates}
\email{oluwaseun.sharomi@kustar.ac.ae}

\address{Rachad Zaki \newline
Department of Applied Mathematics and Sciences,
Khalifa University of Science, Technology and Research\\
Abu Dhabi, United Arab Emirates}
\email{rachad.zaki@kustar.ac.ae}

\thanks{Submitted June 1, 2015. Published November 4, 2015.}
\subjclass[2010]{35K57, 49K20, 91B76}
\keywords{Solow model; pollution; economic growth; spatial diffusion}

\begin{abstract}
 A deterministic model of economic growth and pollution accumulation,
 in the form of a system of partial differential equations,
 is designed and analyzed. The model assumes pollution as a by-product of
 production. The stock of pollution has a negative impact on production.
 The accumulation of pollution is dampened by a share of the investments,
 in the form of an environmental tax. We consider a linear region where
 both capital and pollution can diffuse. This economic-environmental
 model is described by a pair of partial  differential equations whose
 dynamics and steady state characteristics with respect to time and
 space are studied. Then we look at this ambient environment from the
 point of view of a social planner  who can act on the consumption and
 taxation, also functions of time and space, considering the dynamics of
 capital and pollution as constraints.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{definition}[theorem]{Definition}
\allowdisplaybreaks

\section{Introduction}

In this paper we connect two recent strands of economic literature. 
The first strand considers the joint evolution of economic growth and pollution. 
In the recent decades a growing interest both  in the preservation of the 
natural environment and in the long run sustainability of the  economic  
growth has compelled researchers to devise  models that could investigate 
the econo-environment interactions, make predictions, and contrive recommendations
for the formulation of optimal policies. Early works in this regard can be 
traced back to the  70s and early 80s, see for example 
\cite{Becker1972,Brock1973,Forster1973,Gruver1976, Keeler1971,Luptacik1982, Maler1974}. 
Environmental  pollution enters   neoclassical growth models both as a joint 
product and as a source of disutility. There are many different approaches 
in the literature to mathematically model the  interconnections between the 
economy and the environment. Nevertheless they can be grouped in two clusters. 
It is possible to consider pollution as an input of production, assuming that 
the more the pollution is allowed, the less costly are the techniques of 
production;  alternatively, a damage function can represent  the negative 
effect the pollution can have on production.  
See \cite{Bovemberg1995,Mohtadi1996,Rubio2000,Smulders1996} for examples of 
the former approach and 
\cite{Capasso2013,Sebast2015, Anita201536,Capasso20103858,LaTorre2015} 
for an example of the latter. 
Our model formulation is in line with those 
presented in \cite{Capasso2013,Anita201536}; with respect to them, 
in our formulation, the level of taxation is affecting the level of pollution 
(not the level of physical capital) and we have an extra term in the
 pollution equation that justifies the positive effect of taxation on 
abatement activities. This makes our model more complicated, due to the 
presence of quasimonotone functions. Furthermore, the present model does 
not have an integral term and the objective function only includes a combination
of consumption of physical capital and the level of  pollution.

 The second strand of the literature is represented by the recent attempts 
into extending some models of economic growth into the spatial dimension, 
taking into account both the temporal and the spatial dynamics of capital. 
The idea of this spatio-temporal approach has been mutated from the so-called 
New Economic Geography, whose founding father is the nobel awarded economist 
Paul Krugman \cite{Krugman1991,Krugman1993}.

 In a nutshell, we present a model that combines the evolution of the economy 
and its effects on the environment in the natural spatio-temporal ambience. 
The vector that connects the economy and the environment is pollution. 
The production of capital is negatively affected by the stock of pollution, 
which in turn is reduced by a tax proportional to the level of  production. 
The  environmental tax can be an answer to the pressing objective of ensuring 
a reasonable and sustainable level of pollution. This is the rationale for the 
inclusion of such `green tax' in our model, in order to allow the policy makers 
to employ a portion of the investment to the reduction of the quantity of 
pollution per unit of production (see \cite{Xepapadeas1997} for an overview on 
applied principles of environmental policies). For models involving two connected  
modules, applying the same ideas to a slightly different context,  
see  \cite{LaTorre20101017,Marsiglio20121197,LaTorre2011}.

 This article is organized as follows: The model is described in Section 2. 
The dynamic analysis of the model is presented in Section 3.
 An associated optimal control problem is formulated and analyzed in Section 4.
 Section 5 presents the numerical simulations, and the results are discussed 
in Section 6.

\section{The Model}

In the first part of the paper we stick to the hypothesis that  the decisions 
about the consumption share of production  ($c$) and the choice of the 
environmental tax share ($\tau$) are not available to the policy makers 
($c$ and $\tau$ are given parameters). We drop this limitation in a following 
section, where it will be possible for the policy maker to choose the 
consumption and taxation following the path of an optimal control approach. 
Two partial differential equations are the building blocks of our model, given by:
\begin{equation}\label{model}
\begin{gathered}
\frac{\partial k}{\partial t}(x,t)=d_k\Delta k 
+\frac{g(k)(1-\tau-c)}{1+\theta_p p}-\delta_k k, \\
\frac{\partial p}{\partial t}(x,t)=d_p\Delta p 
+\frac{\sigma g(k)}{1+\theta_\tau \tau}-\delta_p p,
\end{gathered}
\end{equation}
where $k$ and $p$ designate capital and pollution, respectively. 
The reaction-diffusion system \eqref{model} describes the mutual 
interaction between the economy and the environment. The first equation 
of System \eqref{model} takes into account  the evolution of the capital $k$. 
The dynamics of the capital at  position $x$ depend on the production function 
acting at $x$ and the contribution of the diffusive term, $\Delta k $. 
The capital can diffuse in space, which means that the producers can decide 
to move their plants to positions where they anticipate better returns. 
The higher the diffusion coefficient $d_k$, the easier it is to move the 
capital from place to place ($d_k$ describes the level of free circulation 
of the capital that the policy maker allows between positions $x$ and $x+dx$). 
 In the neoclassical economic theory better returns  realize where less capital 
has been accumulated, thanks to the law of decreasing returns, an immediate 
consequence of the convexity of the production function. This is still true 
for a convex-concave production function, provided that  the level of 
accumulated capital is above the poverty trap threshold. 
Not all the outcome of the process of production is invested in the accumulation 
of new capital. A fraction of this outcome, $c$, is dedicated to consumption. 
 As a pollution counteracting measure, another share of production, $\tau$,  
is devoted to an environmental tax whose amount would be invested to reduce 
the impact of pollution on the environment.  This is a version of the Solow Model 
of Growth Theory \cite{Solow1956,Solow1974} and of its recent extension 
to the spatial dimension \cite{Camacho2004},  where a constant fraction 
of production, $1-c-\tau$, is invested in the accumulation of the capital. 
It is important to underline how the parameter $A$, the total factor productivity, 
reads in our model. Usually this parameter conveys information about the
 effect of technology on the production function, and it can be either 
exogenously or endogenously determined. In our model, 
$A=\frac{1}{1+\theta_pp}$; that is, we allow the pollution to negatively 
impact the production through a term that can be interpreted as a 
factor inhibiting the performance of the economy, or a damage function. 
This formulation says that if pollution is zero, then there is no 
externality on production. Otherwise the production decreases proportionally 
with the increase in pollution.  One can imagine, for example, that high 
levels of pollution, by destroying environmental amenities, make disconsolate 
and/or less productive human beings (see \cite{Zivin2012} for a study of
 the effect of the environmental pollution on the productivity of labor).
 The depreciation of the capital is taken into account by the term $\delta_kk$,
 a rather standard assumption.

The second equation in  \eqref{model} describes the accumulation and diffusion
of pollution $p$. Pollution is a by-product of the production, as is clear 
from the source term  $\frac{\sigma g(k)}{1+\theta_{\tau}\tau}$.  
The level of pollution at position $x$ is also given by the amount of 
pollution reaching $x$ through the process of diffusion. 
The current model considers a composite pollutant, whose diffusion 
properties are summed up by the diffusion coefficient $d_p$. 
One can think about a combination of greenhouse gases. Moreover  we observe 
that the flow of pollution per unit of production,  
$\frac{\sigma}{1+\theta_{\tau}\tau}$, depends on the level of  
the environmental tax share. The model assumes that the resources 
collected through the environmental taxation are employed to develop 
cleaner industrial processes or abatement activities that facilitate 
the reduction of pollution level, given the same amount of  production. 
The term $\delta_pp$ describes the self-cleaning capacity of the environment. 
In a more realistic approach the exponential decay of pollution must be 
accompanied by a nonlinear feedback term that recounts the possible 
irreversibilities and hysteresis connected to the environmental degradation. 
Indeed, it is not always possible to restore the initial conditions of the 
environment by stopping the economic activity that deteriorated the 
environment in the first place.

 The general form of the production function $g(k)$ is 
\begin{equation}\label{prodfunc}
g(k)=\frac{\alpha_1k^n}{1+\alpha_2 k^n},\quad
n\in\mathbb{Z},\;\alpha_1,\alpha_2\geq 0.
\end{equation}
We consider two different shapes of the production function:
\begin{itemize}
  \item $S-$shaped production function with $n=2$ and $\alpha_2\neq0$;
  \item Concave production function with $n=1$ and $\alpha_2\neq0$.
\end{itemize}
The standard neoclassical (concave) production function, with $n=1$ 
and $\alpha_2\neq0$,  has  convenient  properties from the point 
of view of the neoclassical economic theory: Positivity and decreasing 
return to capital. Yet there are circumstances that are not properly 
modeled if some departures from pure concavity are not allowed for. 
Hence, following the famous idea of Skiba \cite{Skiba1978}, we allow 
for the existence of non-concavity. With  $n=2$ and $\alpha_2\neq0$ 
the function $g(k)$  is a convex-concave ($S-$shaped) production function, 
meaning that for values of $k$ up to a certain threshold the function 
exhibits convexity, and then concavity (it has the so-called $S-$shaped form).  
In terms of the economic literature, the function exhibits increasing and then 
decreasing return to capital. The $S-$shaped curve  is not a pure neoclassical 
production function because it does not respect the law of diminishing return 
for all the values of $k$ (the second derivative is not always negative), 
but it gives rise to richer dynamics.

\section{Dynamic Analysis}

Let $D=\Omega\times[0,T]$, where $\Omega$ is a bounded domain. Consider the problem
\begin{equation} \label{model1}
\begin{gathered}
\frac{\partial k}{\partial t}(x,t)=d_k\Delta k +f_1(k,p) \quad \text{in }D,
\\
\frac{\partial p}{\partial t}(x,t)=d_p\Delta p +f_2(k,p) \quad \text{in }D,
\end{gathered}
\end{equation}
where
\[
f_1(k,p)=\frac{g(k)(1-\tau-c)}{1+\theta_p p}-\delta_k k,
\quad
f_2(k,p)=\frac{\sigma g(k)}{1+\theta_\tau \tau}-\delta_p p,
\]
with $g$ given by \eqref{prodfunc}.

\begin{table}[htb]
\caption{Description and values of the model parameters} \label{table1}
\begin{center}
\begin{tabular}{l|l|l}
\hline
Parameters & Description & Typical value (range)\\
\hline
$\frac{\alpha_1}{\alpha_2}$ & Maximum production level  & 10\\
$1/\alpha_2$ & Half-saturation constant & $\alpha_2=0.1$\\
$\delta_k$ & Capital depreciation rate & 0.02\\
$\theta_p$ & Trade-off parameter & 0.001\\
$\theta_\tau$ & Trade-off parameter & 2\\
$\sigma$ & Trade-off parameter & 2\\
$\delta_p$ & Pollution abatement rate & 0.4\\
$\tau$ & Green tax rate & Varies\\
\hline
\end{tabular}
\end{center}
\end{table}

 The boundary conditions are the homogeneous Neumann conditions
\begin{equation} \label{bc}
\begin{gathered}
\frac{\partial k}{\partial n}=0  \quad \text{for }x\in\partial\Omega,\;t\in[0,T],
\\
\frac{\partial p}{\partial n}=0 \quad  \text{for }x\in\partial\Omega,\;t\in[0,T],
\end{gathered}
\end{equation}
where $n$ designates a unit outward normal vector to $\partial\Omega$. 
The initial conditions are given by
\begin{equation} \label{ic}
\begin{gathered}
k(x,0)=k_0(x) \quad \text{for }x\in\Omega, \\
p(x,0)= p_0(x) \quad\text{for }x\in\Omega,
\end{gathered}
\end{equation}
where functions $k_0$ and $p_0$ are smooth. The goal is to prove the existence
of a solution of this problem (we will also prove the uniqueness in this case) 
and to study the steady-state solutions. We start by introducing several 
definitions and results that will be useful for the study of 
Problem \eqref{model1}--\eqref{ic}. The process we will follow is based 
on finding what is referred to as lower and upper solutions in order to 
prove the existence and, under suitable hypotheses, uniqueness of the 
solution of Problem \eqref{model1}--\eqref{ic}.
This technique is introduced in details in \cite{pao92} for the systems
of parabolic as well as elliptic partial differential equations.
The current model contains mixed quasimonotone functions, which will be
elaborated in the following section.

\subsection*{Preliminary results}

To prove the existence of a solution of System \eqref{model1}, subject to 
the boundary and initial conditions \eqref{bc} and \eqref{ic}, we need to 
introduce some preliminary results. These results can be found in more 
details in \cite{pao92}. Consider the system
\begin{equation}\label{system}
\begin{gathered}
\frac{\partial u_i}{\partial t}=\mathcal{L}_i u_i+ g_i(x,t,u_1,u_2) \quad \text{in }D,
\\
\mathcal{B}_i u_i=h_i(x,t) \quad \text{for }x\in\partial\Omega,\;t\in[0,T],
\\
u_i(x,0)=u_{i,0}(x) \quad \text{for }x\in\Omega,
\end{gathered}
\end{equation}
for $i=1,2$, where $u_{i,0}\in C_{L^\infty}^0(\Omega)$ and the operators 
$\mathcal{L}_i$ are uniformly elliptic with H\"{o}lder continuous coefficients 
and having the form
\[
\mathcal{L}_i\equiv\sum_{j,l=1}^n a_{j,l}^{(i)}(x,t)
\frac{\partial^2}{\partial x_j\partial x_l}
+\sum_{j=1}^n b_j^{(i)}(x,t)\frac{\partial}{\partial x_j},
\]
and $\mathcal{B}_i$ are defined as
\[
\mathcal{B}_i\equiv\beta_i(x,t)\frac{\partial}{\partial n}+\gamma_i(x,t),
\]
where $\beta_i$ and $\gamma_i$ are continuous for $i=1,2$, $\beta_i\geq0$,
$\gamma_i\geq 0$, $\beta_i+\gamma_i>0$, and $n$ being the unit outward
normal vector to $\partial\Omega$.

\begin{definition} \label{def1} \rm
$(g_1,g_2)=(g_1(r,s), g_2(r,s))$ is called mixed quasimonotone when
\[%\label{mixed}
\frac{\partial g_1}{\partial s}\leq 0 \quad\text{and}\quad
\frac{\partial g_2}{\partial r}\geq 0,
\]
or vice versa.
\end{definition}

\begin{definition} \label{def2} \rm
Suppose that $(g_1,g_2)$ is mixed quasimonotone. We call 
$\overline{u}=(\overline{u_1},\overline{u_2})$ and
 $\underline{u}=(\underline{u_1},\underline{u_2})$ in 
$\left(C_{L^\infty}^{2,1}(\Omega\times[0,T])\right)^2$ coupled ordered upper 
and lower solutions of \eqref{system}
 if they satisfy the following relations for $i=1,2$:
\begin{gather*} %\label{uplow}
\overline{u}\geq \underline{u},
\\
\mathcal{B}_i\overline{u}_i \geq h_i(x,t) \geq \mathcal{B}_i\underline{u}_i,
\\
\overline{u}_i(x,0)\geq u_{i,0}(x)\geq  \underline{u}_i(x,0),
\\
\frac{\partial \overline{u}_1}{\partial t}-\mathcal{L}_1 \overline{u}_1- g_1(x,t,\overline{u}_1,\underline{u}_2)
\geq 0
\geq \frac{\partial \underline{u}_1}{\partial t}-\mathcal{L}_1 \underline{u}_1
 - g_1(x,t,\underline{u}_1,\overline{u}_2),
\\
\frac{\partial \overline{u}_2}{\partial t}-\mathcal{L}_2 \overline{u}_2
 - g_2(x,t,\overline{u}_1,\overline{u}_2)
\geq 0
\geq \frac{\partial \underline{u}_2}{\partial t}-\mathcal{L}_2 \underline{u}_2
 - g_2(x,t,\underline{u}_1,\underline{u}_2).
\end{gather*}
\end{definition}


\begin{theorem} \label{thm3.1}
Let $(\overline{u}_1,\overline{u}_2)$ and $(\underline{u}_1,\underline{u}_2)$ 
be coupled upper and lower solutions of \eqref{system} and let $(g_1,g_2)$ 
be mixed quasimonotone in
\[
\langle \underline{u},\overline{u}\rangle
:=\big\{(u_1,u_2)\in\big(C_{L^\infty}^{2,1}(\Omega\times[0,T])\big)^2:
(\underline{u}_1,\underline{u}_2)\leq(u_1,u_2)\leq(\overline{u}_1,\overline{u}_2)\big\}.
\]
Then, Problem \eqref{system} has a unique solution $u$ in
$\langle \underline{u},\overline{u}\rangle$.
\end{theorem}

\subsection*{Existence and uniqueness of a solution}

\begin{theorem} \label{thm3.2}
Problem \eqref{model1} subject to the conditions \eqref{bc} and \eqref{ic}
 has a solution and this solution is unique.
\end{theorem}

\begin{proof}
We first observe that $f=(f_1(k,p),f_2(k,p))$ is mixed quasimonotone. Indeed
\begin{gather*}
\frac{\partial f_1}{\partial p}
= g(k)(1-\tau-c)\frac{(-\theta_p)}{(1+\theta_p p)^2}\leq 0,\\
\frac{\partial f_2}{\partial k}
= \frac{\sigma}{1+\theta_\tau \tau}g'(k)
=\frac{\sigma}{1+\theta_\tau \tau}\frac{n\alpha_1 k^{n-1}}{(1+\alpha_2k^n)^2}\geq 0.
\end{gather*}
We now prove the existence of lower and upper solutions of \eqref{model1}, 
denoted by $(\underline{k},\underline{p})$ and $(\overline{k},\overline{p})$ 
respectively, satisfying the following conditions:
\begin{gather*} %\label{uplow_kp}
\overline{k}\geq \underline{k}, \quad 
\overline{p}\geq \underline{p},
\\
\frac{\partial \overline{k}}{\partial n} 
\geq 0 \geq \frac{\partial \underline{k}}{\partial n},
\\
\frac{\partial \overline{p}}{\partial n}
\geq 0 \geq \frac{\partial \underline{p}}{\partial n},
\\
\overline{k}(x,0)
\geq k_0(x)
\geq \underline{k}(x,0),
\\
\overline{p}(x,0)
\geq p_0(x)
\geq \underline{p}(x,0),
\\
\frac{\partial \overline{k}}{\partial t}-d_k\Delta\overline{k}
- f_1(\overline{k},\underline{p})
\geq 0
\geq \frac{\partial \underline{k}}{\partial t}-d_k\Delta\underline{k}
 - f_1(\underline{k},\overline{p}),
\\
\frac{\partial \overline{p}}{\partial t}-d_p\Delta \overline{p}
 - f_2(\overline{k},\overline{p})
\geq 0
\geq \frac{\partial \underline{p}}{\partial t}-d_p\Delta \underline{p}
 - f_2(\underline{k},\underline{p}).
\end{gather*}
Let $\underline{k}=0$ and note that it satisfies 
$ \frac{\partial \underline{k}}{\partial t}
\leq d_k\Delta \underline{k}+ f_1(\underline{k},\overline{p})$. 
Using $\underline{k}=0$, $\underline{p}$ satisfies
\begin{gather*}
\frac{\partial \underline{p}}{\partial t}
\leq d_p\Delta \underline{p}- \delta_{p}\underline{p},
\\
\frac{\partial \underline{p}}{\partial n}
\leq 0,
\\
\underline{p}(x,0)
\leq p_0.
\end{gather*}
It is reasonable to choose $\underline{p}=\min\{0,\inf_{\Omega}p_0\}$, 
which leads to $\underline{p}=0$ (since $p_0(x,t)\geq 0$ for all valid $(x,t)$,
 by definition). Considering $\overline{k}$ independent of $x$ and noticing 
that $\underline{p}=0$ the upper solution $\overline{k}$ satisfies
\begin{gather*}
\frac{\partial \overline{k}}{\partial t}
\geq g(\overline{k})(1-\tau-c)-\delta_{k}\overline{k},
\\
\frac{\partial \overline{k}}{\partial n}
\geq 0,
\\
\overline{k}(x,0)
\geq k_0.
\end{gather*}
Denoting $a:=\alpha_1(1-\tau-c)$,
\[
\frac{\partial \overline{k}}{\partial t}
\geq h(\overline{k})-\delta_{k}\overline{k},
\]
where $h(x)=\frac{ax^n}{1+\alpha_2x^n}$. Since $h$ is increasing
for $x>0$ and $\lim_{x\to\infty}h(x)=\frac{a}{\alpha_2}$,
it is sufficient to look for a $\overline{k}$ that satisfies
\[
\frac{\partial \overline{k}}{\partial t}
\geq \frac{a}{\alpha_2}-\delta_{k}\overline{k}.
\]
Therefore, we consider the initial value problem
\begin{equation}\label{ivp}
\begin{gathered}
y'(t)=\frac{a}{\alpha_2}-\delta_k y,
\\
y(0)=\sup_\Omega k_0.
\end{gathered}
\end{equation}
The solution of \eqref{ivp} is
\[
y(t)= \frac{a}{\delta_k\alpha_2}+\Big(\sup_\Omega k_0
 -\frac{a}{\delta_k\alpha_2}\Big)e^{-\delta_k t}.
\]
Hence we can choose
\[
\overline{k}
=\max\big\{\frac{\alpha_1(1-\tau-c)}{\delta_k\alpha_2}, \sup_\Omega k_0\big\}.
\]
Finally, suppose $\overline{p}$ is also space-independent.
Then $\overline{p}$ must satisfy
\[
\frac{\partial \overline{p}}{\partial t}
\geq \frac{\sigma g(\overline{k})}{1+\theta_\tau \tau}-\delta_{p}\overline{p},
\]
along with the corresponding boundary and initial conditions.
To identify $\overline{p}$, we similarly look for a solution of
\begin{equation}\label{ivp_p}
\begin{gathered}
z'(t)=\mu-\delta_{p}z,
\\
z(0)=\sup_\Omega{p_0},
\end{gathered}
\end{equation}
where
\[
\mu=\frac{\sigma g(\overline{k})}{1+\theta_\tau \tau}.
\]
The solution of \eqref{ivp_p} is
\[
z(t)=\frac{\mu}{\delta_{p}}+\Big(\sup_\Omega{p_0}-\frac{\mu}{\delta_{p}}\Big)
 e^{-\delta t}.
\]
One can therefore choose
\[
\overline{p}
= \max\Big\{\frac{\sigma g(\overline{k})}{\delta_{p}(1+\theta_\tau \tau)},
\sup_\Omega{p_0}\Big\}.
\]
We can now apply Theorem \ref{thm3.1} to deduce the existence of a unique
solution $u$ to \eqref{model1}-\eqref{ic}, with
\[
u\in \Big\langle
(0,0),
\Big(\max\big\{\frac{\alpha_1(1-\tau-c)}{\delta_k\alpha_2}, \sup_\Omega k_0\big\} ,
\max\big\{\frac{\sigma g(\overline{k})}{\delta_{p}(1+\theta_\tau \tau)},
  \sup_\Omega{p_0}\big\}
\Big) \Big\rangle.
\]
\end{proof}

A solution profile of $k(x,t)$ and $p(x,t)$ is provided in Figures \ref{fig1} 
and \ref{fig2}.
\begin{figure}[htb]
\begin{center}
 \includegraphics[width=0.8\textwidth]{fig1} % cbshapedpde1
    \end{center}
\caption{Plot of $k(x,t)$ and $p(x,t)$ using the concave production function}
\label{fig1}
\end{figure}

\begin{figure}[htb]
\begin{center} 
\includegraphics[width=0.8\textwidth]{fig2} % sshapedpde1
\end{center}
 \caption{Plot of $k(x,t)$ and $p(x,t)$ using the $S-$shaped production function}
\label{fig2} 
\end{figure}

The steady state solution of \eqref{model1}-\eqref{ic}
 is a pair of smooth functions $(k^\star,p^\star)$ satisfying
\begin{gather*} %\label{steady}
d_k\Delta k^\star +f_1(k^\star,p^\star)=0 \quad \text{in }\Omega,
\\
d_p\Delta p^\star +f_2(k^\star,p^\star)=0 \quad \text{in }\Omega,
\\
\frac{\partial k^\star}{\partial n}=\frac{\partial p^\star}{\partial n}=0
 \quad \text{on }\partial \Omega.
\end{gather*}
This problem, although involving elliptic partial differential equations, 
can be addressed in a similar way as \eqref{model1}-\eqref{ic}
 by looking for upper and lower solutions in order to deduce the existence 
of solutions of the problem; that is, equilibrium solutions of 
\eqref{model1}-\eqref{ic}. The tools needed for treating this problem 
are similar to the ones used in the previous section and can be 
found in \cite{pao92}.
The steady state solution profiles of $k^*$ and $p^*$ are provided in 
Figures \ref{fig3} and \ref{fig4}.

\begin{figure}[htb]
\begin{center}
\includegraphics[width=0.8\textwidth]{fig3} % cbshapedsteady
\end{center}
    \caption{Plot of $k^*$ and $p^*$ using the concave production function}
\label{fig3}
\end{figure}

\begin{figure}[htb]
\begin{center}
\includegraphics[width=0.8\textwidth]{fig4} % sshapedsteady
\end{center}
    \caption{Plot of $k^*$ and $p^*$ using the $S-$shaped production function}
\label{fig4}
\end{figure}


\section{Optimal control}

So far the consumption share  $c$ and the taxation share $\tau$ have been 
considered exogenous, meaning that their values have been treated as the 
result of exogenous choices. Now we want to allow a social planner to choose 
these parameters in such a way to optimize an objective function. 
The objective function takes into account consumption as a source of utility 
and pollution as a source of disutility. The planner's problem reads
\[
\max_{\{ c(x,t),\tau(x,t)\}}\int_{0}^{T}
\int_a^b[\theta c(x,t)g(k(x,t))-\gamma p(x,t)]\,dx\,dt,
\]
subject to
\begin{equation}\label{pde1}
\begin{gathered}
\frac{\partial}{\partial t}k(x,t)
=d_k\frac{\partial^2 k(x,t)}{\partial x^2}+\frac{g(k(x,t))
[1-\tau(x,t)-c(x,t)]}{1+\theta_pp(x,t)}-\delta_k k(x,t),
\\
\frac{\partial}{\partial t}p(x,t)
=d_p\frac{\partial^2 p(x,t)}{\partial x^2}
+\frac{\sigma g(k(x,t))}{1+\theta_\tau \tau(x,t)}-\delta_pp(x,t).
\end{gathered}
\end{equation}
The planner can choose whatever values for $c$ and $\tau$ (s)he 
considers optimal, provided that these two control variables are picked in 
a reasonable control set. $c$ and $\tau$ are shares, so they need 
to belong to the set $[0,1]$. Moreover it is plausible to think that there 
is a bottom limit for the investment in new capital that cannot be crossed. 
In other words, the control set has been built in such a way that consumption 
and tax shares do not exhaust the investment; that is, 
$ c+\tau \leq\theta_{c \tau}$. We chose $\theta_{c \tau}=0.8$, but other values 
do not affect our qualitative results. The Hamiltonian function $\mathcal{H}$ is
\begin{equation}\label{hamilt}
\begin{aligned}
\mathcal{H}
&=\theta c(x,t)g(k(x,t))-\gamma p(x,t)\\
&\quad +\lambda_k(x,t)\Big[d_k\frac{\partial^2 k(x,t)}{\partial x^2}+\frac{g(k(x,t))[1-\tau(x,t)-c(x,t)]}{1+\theta_pp(x,t)}-\delta_kk(x,t)\Big]\\
&\quad +\lambda_p(x,t)\Big[d_p\frac{\partial^2 p(x,t)}{\partial x^2}+\frac{\sigma g(k(x,t))}{1+\theta_\tau \tau(x,t)}-\delta_pp(x,t)\Big].
\end{aligned}
\end{equation}
The first order conditions to the previous problem are
\begin{equation}\label{pollution}
\begin{gathered}
\begin{aligned}
\frac{\partial }{\partial t}\lambda_{k}(x,t)
&=-d_k\frac{\partial^2 \lambda_{k}(x,t)}{\partial x^2}-\theta c(x,t)g_k(k(x,t))\\
&\quad-\lambda_k(x,t)g_k(k(x,t))\frac{[1-\tau(x,t)-c(x,t)]}{1+\theta_pp(x,t)}\\
&\quad -\lambda_p(x,t)g_k(k(x,t))\frac{\sigma}{1+\theta_{\tau}\tau(x,t)}
 +\delta_k \lambda_k(x,t),
\end{aligned}\\
\begin{aligned}
\frac{\partial }{\partial t}\lambda_{p}(x,t)
&=-d_p\frac{\partial^2 \lambda_{p}(x,t)}{\partial x^2}+ \gamma\\
&\quad +\lambda_k(x,t)g(k(x,t))\frac{\theta_p[1-\tau(x,t)-c(x,t)]}{[1+\theta_pp(x,t)]^2}\\
&\quad +\delta_p \lambda_{p}(x,t),
\end{aligned}
\end{gathered}
\end{equation}
subject to homogeneous Neumann boundary conditions and the final conditions
\[
\lambda_k(x,T)=\lambda_p(x,T)=0.
\]

\section{Numerical Simulations}
\subsection*{Description of the algorithm}
We implemented the forward-backward sweep method for System \eqref{pde1}, 
\eqref{hamilt} and \eqref{pollution} as follows:

(1) Choose an initial guess: $(c^{(0)},\tau^{(0)})=(c^{(0)}(t),\tau^{(0)}(t))$.
\smallskip

(2) Iterate for $j\ge0$: Using the spectral method, we solved
\begin{align*}
\frac{\partial k^{(j+1)}(x,t)}{\partial t}
&=d_k\frac{\partial^2 k^{(j+1)}(x,t)}{\partial x^2}
 +\frac{g(k^{(j+1)}(x,t))[1-\tau^{(j)}(x,t)-c^{(j)}(x,t)]}{1+\theta_pp^{(j+1)}(x,t)}\\
&\quad -\delta_k k^{(j+1)}(x,t),\\
\frac{\partial p^{(j+1)}(x,t)}{\partial t}
&=d_p\frac{\partial^2 p^{(j+1)}(x,t)}{\partial x^2}
 +\frac{\sigma g(k^{(j+1)}(x,t))}{1+\theta_\tau \tau^{(j)}(x,t)}
 -\delta_pp^{(j+1)}(x,t),
\end{align*}
subject to
\begin{gather*}
k^{(j+1)}(x,0)=k_0(x)  \quad\text{in }\Omega,\\
p^{(j+1)}(x,0)=p_0(x)  \quad \text{in }\Omega,\\
\frac{\partial k^{(j+1)}(x,t)}{\partial n}=0  \quad \text{for }x\in\partial \Omega,\\
\frac{\partial p^{(j+1)}(x,t)}{\partial n}=0 \quad  \text{for }x\in\partial \Omega,
\end{gather*}
from $t=0$ to $t=T$. We reversed the equations \eqref{pollution} in time, via the change of variable $\bar{t}=T-t$, turning the problem into a forward problem with zero initial conditions. Then, we solved
\begin{align*}
\frac{\partial \lambda_{k}^{(j+1)}(x,\bar{t})}{\partial \bar{t}}
&=d_k\frac{\partial^2 \lambda_{k}^{(j+1)}(x,\bar{t})}{\partial x^2}
 +\theta c^{(j)}(x,\bar{t})g_k(k^{(j+1)}(x,\bar{t}))\\
&\quad +\lambda_k^{(j+1)}(x,\bar{t})g_k(k^{(j+1)}(x,\bar{t}))
 \frac{[1-\tau^{(j)}(x,\bar{t})-c^{(j)}(x,\bar{t})]}{1+\theta_pp^{(j+1)}(x,\bar{t})}\\
&\quad +\lambda_p^{(j+1)}(x,\bar{t})g_k(k^{(j+1)}(x,\bar{t}))
 \frac{\sigma}{1+\theta_{\tau}\tau^{(j)}(x,\bar{t})}-\delta_k\lambda_k^{(j+1)}
 (x,\bar{t}),
\end{align*}
\begin{align*}
\frac{\partial \lambda_{p}^{(j+1)}(x,\bar{t})}{\partial \bar{t}}
&=d_p\frac{\partial^2 \lambda_{p}^{(j+1)}(x,\bar{t})}{\partial x^2}-\gamma\\
&\quad -\lambda_k^{(j+1)}(x,\bar{t})g(k^{(j+1)}(x,\bar{t}))\frac{\theta_p[1-\tau^{(j)}(x,\bar{t})-c^{(j)}(x,\bar{t})]}{[1+\theta_pp^{(j+1)}(x,\bar{t})]^2}\\
&\quad -\delta_p \lambda_{p}^{(j+1)}(x,\bar{t}).
\end{align*}
subject to
\begin{gather*}
\lambda_{k}^{(j+1)}(x,0)=0 \quad \text{in }\Omega,\\
\lambda_{k}^{(j+1)}(x,0)=0 \quad \text{in }\Omega,\\
\frac{\partial \lambda_{k}^{(j+1)}(x,\bar{t})}{\partial n}=0  
\quad \text{for }x\in\partial \Omega,\\
\frac{\partial \lambda_{k}^{(j+1)}(x,\bar{t})}{\partial n}=0 
\quad \text{for }x\in\partial \Omega,
\end{gather*}
from $\bar{t}=0$ to $\bar{t}=T$. We used the `fmincon' function of 
MATLAB (dedicated to finding the minimum of a constrained nonlinear
 multivariable function) defined below
 \[
 \min_{x}f(x)\text{ such that }
 \begin{cases}
 c(x)\leq 0,\\
 ceq(x) = 0,\\
 A\cdot x\le b,\\
 Aeq\cdot x = beq, \\
 lb \le x \le ub,
 \end{cases}
 \]
to determine the values of $c(x,t)$ and $\tau(x,t)$ that maximize
$\mathcal{H}$. We achieved this by finding the values of $c(x,t)$ and
 $\tau(x,t)$ that minimized $-\mathcal{H}$. Here, we define
\begin{gather*}
A=  \begin{pmatrix}
    1 & 0 \\
    0 & 1 \\
    1 & 1 \\
    -1 & -1
  \end{pmatrix}, \quad
b=  \begin{pmatrix}
    1 \\
    1 \\
    0.8 \\
    0
  \end{pmatrix}, \\
lb=  \begin{pmatrix}
    0 \\
    0
  \end{pmatrix}, \quad
ub=  \begin{pmatrix}
    1 \\
    1
  \end{pmatrix}, \quad
x=  \begin{pmatrix}
    c^{(j+1)}(x,t) \\
    \tau^{(j+1)}(x,t)
  \end{pmatrix}.
\end{gather*}

(3)  We checked convergence by computing the difference between the 
respective values of $c(x,t)$, $\tau(x,t)$, $k(x,t)$, $p(x,t)$, 
$\lambda_k(x,t)$ and $\lambda_p(x,t)$ in two consecutive iterations. 
If the maximum of the $L^2$-norm of the difference was negligibly small, 
we output the current function as a solution, otherwise we continued iterating.

\begin{figure}[htb]
\begin{center}
 \includegraphics[width=0.8\textwidth]{fig5} % spaceplot-cb
\end{center}
\caption{Plot of $k(x,t), p(x,t), c(x,t)$ and $\tau(x,t)$ using the
 concave production function} \label{fig5}
\end{figure}

\begin{figure}[htb]
\begin{center}
 \includegraphics[width=0.8\textwidth]{fig6} %spaceplot-ss
\end{center}
\caption{Plot of $k(x,t), p(x,t), c(x,t)$ and $\tau(x,t)$ using the
$S$-shaped production function} \label{fig6}
\end{figure}



\section{Conclusions}

In this article we have studied the dynamics in space and time of  
a coupled environment-growth model. In the first part the analysis focused
 on the existence and uniqueness of solutions. This  first part can be 
considered a  short-run analysis of the economy with fixed values for 
the control variables. Then we  allowed for these variables to be chosen by 
a social planner. The social planner wants to optimize the objective function 
given the partial differential equations describing the coupled environmental 
economic system. We used a generalized version of Pontryagin's maximum principle 
as in \cite{Xepapadeas1997}, but our analysis went through the complete 
dynamical study of the solutions. In the spatial environmental economic literature, 
see \cite{Krugman1991,Krugman1993,Fujita1999,Fujita2002,Camacho2004,Brito2004,
Brito2005,Camacho2008,Boucekkine2009,capasso2010,Boucekkine2013}, 
only the study of the solutions around a space-homogeneous steady state was performed.
 By applying the Sweep Algorithm to our framework, we have been able to perform 
numerical simulations of all the spatio-temporal path.

As we can see from Figures \ref{fig5} and \ref{fig6}, both in the case of an 
concave production function and in the case of a $S$-shaped function, the diffusion 
creates spatial homogeneity: The spatial heterogeneous  initial time profiles 
of capital and pollution are smoothed out. As for the temporal dynamics, we see 
that both capital and pollution grow along a sustainable path. The social planner 
finds it optimum to let consumption share  grow over time, while the green taxation 
share decreases simultaneously. In other words, it is optimum to first dedicate 
the major part of the available investment resources to abatement activities, 
and then progressively let consumption increase while the taxation share is
 being reduced. This model proposes a precise suggestion for the consumption 
and taxation policy to be followed when the time horizon is set to $T$: 
Consumption is supposed to increase slowly, giving the time to the abatement 
activities financed by the  taxation share to do its job, namely driving the
 growth path of capital and pollution toward a sustainable outcome 
(this is in-line with the definition of green taxation that is mainly devoted 
to pollution abatement) and with the recommendations of the Organization for 
Economic Co-operation and Development (OECD) on green growth and taxation. 
As pointed out in \cite{OECD}: ``Environmentally related taxes are increasingly 
being used in OECD economies and can provide significant incentives for 
innovation, as firms and consumers seek new, cleaner solutions in response 
to the price put on pollution. These incentives also make it commercially 
attractive to invest in R\&D activities to develop technologies and consumer 
products with a lighter environmental footprint''.  According to the economic
 conclusion of our model, there exist optimal paths for consumption and taxation 
that allow to reach a sustainable level of pollution, together with a hopefully 
satisfactory level of capital. This is an interesting result from the economic 
perspective and extends similar results in the literature (see \cite{Brock2010}).

\subsection*{Acknowledgements} 
This research was supported by Khalifa University Internal Research Fund 
(Grant No. 210032).


\begin{thebibliography}{00}

\bibitem{Capasso2013} S.~Anita, V.~Capasso, H.~Kunze, D.~La Torre;
 \emph{Optimal control and   long-run dynamics for a spatial economic 
growth model with physical capital   accumulation and pollution diffusion}, 
Applied Mathematics Letters  \textbf{8}   (2013), no.~26, 908--912.

\bibitem{Sebast2015}  S.~Anita, V.~Capasso, H.~Kunze, D.~La Torre;
\emph{Dynamics and control of an integro-differential system of
  geographical economics}, Annals of the Academy of Romanian Scientists: Series
  on Mathematics and its Applications \textbf{7} (2015), no.~1, 8--26.

\bibitem{Anita201536}  S.~Anita, V.~Capasso, H.~Kunze, D.~La Torre;
 \emph{Dynamics and optimal control in a spatially structured economic
  growth model with pollution diffusion and environmental taxation}, Applied
  Mathematics Letters \textbf{42} (2015), 36--40.

\bibitem{Becker1972} R.~A. Becker;
 \emph{Intergenerational equity: The capital-environment
  trade-off}, Journal of Environmental Economics and Management \textbf{9}
  (1982), no.~2, 165--185.

\bibitem{Boucekkine2013} R.~Boucekkine, C.~Camacho, G.~Fabbri;
\emph{Spatial dynamics and   convergence: The spatial \{AK\} model}, 
Journal of Economic Theory   \textbf{148} (2013), no.~6, 2719--2736.

\bibitem{Boucekkine2009} R.~Boucekkine, C.~Camacho, B.~Zou;
 \emph{Bridging the gap between growth
  theory and the new economic geography: the spatial ramsey model},
  Macroeconomics Dynamics \textbf{13} (2009), no.~1, 20--45.

\bibitem{Bovemberg1995} A.~L. Bovenberg and S.~Smulders;
\emph{Environmental quality and   pollution-augmenting technological 
change in a two-sector endogenous growth   model}, Journal of 
Public Economics \textbf{57} (1995), no.~3, 369--391.

\bibitem{Brito2004} P.~Brito;
 \emph{The dynamics of growth and distribution in a spatially
  heterogeneous world}, UECE-ISEG, Technical University of Lisbon (2004).

\bibitem{Brock1973} W.~A. Brock;
 \emph{A polluted golden age}, Economics of Natural and
  Environmental Resources (Gordon Breach, New York) (1973).

\bibitem{Brock2010} W.~A. Brock, M.~S. Taylor;
\emph{The green solow model}, Journal of Economic
  Growth \textbf{15} (2010), no.~2, 127--153 (English).

\bibitem{Brito2005} W. A. Brock, M. S. Taylor;
 \emph{Economic growth and the environment: a review
  of theory and empirics}, (In: Aghion, P., Durlauf, S., Eds., Handbook of
  Economic Growth, 1749-1821) (2005).

\bibitem{Camacho2008} C.~Camacho, B.~Zou, M.~Briani;
 \emph{On the dynamics of capital
  accumulation across space}, European Journal of Operational Research
  \textbf{186} (2008), no.~2, 451--465.

\bibitem{Capasso20103858} V.~Capasso, R.~Engbers,  D.~La Torre;
\emph{On a spatial solow model with   technological diffusion and 
nonconcave production function}, Nonlinear   Analysis: 
Real World Applications \textbf{11} (2010), no.~5, 3858--3876.

\bibitem{capasso2010} V.~Capasso, R.~Engbers,  D.~La Torre;
\emph{On a spatial solow model with technological diffusion and
  nonconcave production function}, Nonlinear Analysis:
 Real World Applications
  \textbf{11} (2010), no.~5, 3858--3876.

\bibitem{Forster1973} B.~A. Forster;
 \emph{Optimal capital accumulation in a polluted environment},
  Southern Economic Journal \textbf{39} (1973), no.~4, pp. 544--547.

\bibitem{Fujita1999} M.~Fujita, P.~Krugman, A.~Venables;
 \emph{The spatial economy. cities},
  regions and international trade (MIT Press) (1999).

\bibitem{Fujita2002} M. Fujita  J.F . Thisse;
 \emph{Economics of agglomeration}, (Cambridge
  University Press) (2002).

\bibitem{Gruver1976} G.~W. Gruver;
 \emph{Optimal investment in pollution control capital in a
  neoclassical growth context}, Journal of Environmental Economics and
  Management \textbf{3} (1976), no.~3, 165--177.

\bibitem{Keeler1971} E.~Keeler, M.~Spence, R.~Zeckhauser;
 \emph{The optimal control of   pollution}, Journal of Economic
 Theory \textbf{4} (1971), 19--34.

\bibitem{Krugman1991} P.~Krugman;
\emph{Increasing returns and economic geography}, Journal of
  Political Economy \textbf{99} (1991), no.~3, 483--499.

\bibitem{Krugman1993} P.~Krugman;
\emph{On the number and location of cities}, European Economic Review
  \textbf{37} (1993), no.~2, 293--298.

\bibitem{Luptacik1982} M.~Luptacik, U.~Schubert;
 \emph{Optimal economic growth and the   environment}, 
Economic Theory of Natural Resources (Physica-Verlag, Wurzburg, Wien), 1982.

\bibitem{LaTorre2015} D.~La Torre, D.~Liuzzi, S.~Marsiglio;
 \emph{Pollution diffusion and  abatement activities across space 
and over time}, Mathematical Social Sciences 78 (2015), 48--63.

\bibitem{LaTorre20101017} D.~La Torre, S.~Marsiglio;
 \emph{Endogenous technological progress in a   multi-sector growth model}, 
Economic Modelling \textbf{27} (2010), no.~5,  1017--1028.

\bibitem{Maler1974} K.~G. Maler;
\emph{Environmental economics: A theoretical inquiry}, Johns
  Hopkins University Press, Baltimore, 1974.

\bibitem{Marsiglio20121197} S.~Marsiglio, D.~La Torre;
 \emph{Population dynamics and utilitarian   criteria in the 
lucas–uzawa model}, Economic Modelling \textbf{29} (2012),
  no.~4, 1197--1204.

\bibitem{LaTorre2011} S.~Marsiglio, D.~La Torre,  F.~Privileggi;
\emph{Fractals and   self-similarity in economics: the case of a stochastic 
two sector growth   model}, Image Analysis and Stereology \textbf{30} (2010), 
no.~3, 143--151.

\bibitem{Mohtadi1996} H.~Mohtadi;
\emph{Environment, growth, and optimal policy design}, Journal of
  Public Economics \textbf{63} (1996), no.~1, 119--140.

\bibitem{OECD} OECD;
\emph{Green growth and taxation},
http://www.oecd.org/greengrowth/greengrowthandtaxation.htm. Date
  accessed: 2015-10-19.

\bibitem{pao92} C.~V. Pao;
\emph{Nonlinear parabolic and elliptic equations}, Plenum Press, New
  York, 1992.

\bibitem{Rubio2000} S.~J. Rubio, J.~Aznar;
 \emph{Sustainable growth and environmental policies},
  Fondazione Enrico Mattei Discussion Paper (2000).

\bibitem{Skiba1978} A.~A. Skiba;
 \emph{Optimal growth with a convex-concave production function},
  Econometrica \textbf{46} (1978), no.~3, 527--539.

\bibitem{Smulders1996}; Sjak Smulders, Raymond Gradus;
 \emph{Pollution abatement and long-term   growth}, 
European Journal of Political Economy \textbf{12} (1996), no.~3, 505--532.

\bibitem{Solow1956} R.~M. Solow;
 \emph{A contribution to the theory of economic growth}, The
  Quarterly Journal of Economics \textbf{70} (1956), no.~1, 65--94.

\bibitem{Solow1974} R. M. Solow;
 \emph{Intergenerational equity and exhaustible resources}, Review
  of Economic Studies \textbf{41} (1974), 29--45.

\bibitem{Xepapadeas1997} A.~Xepapadeas;
 \emph{Advanced principles in environmental policy}, Edward Elgar
  Publishers, Cheltenam, 1997.

\bibitem{Zivin2012} J.~G. Zivin and M.~Neidell;
\emph{The impact of pollution on worker   productivity}, 
American Economic Review \textbf{102} (2012), no.~7, 3652--3673.

\bibitem{Camacho2004} B.~Zou, C.~Camacho;
 \emph{{The spatial Solow model}}, Economics Bulletin
  \textbf{18} (2004), no.~2, 1--11.

\end{thebibliography}


\end{document}