\section{Introduction}

In the paper, we consider an optimal control problem in a dynamical system which motion is described by a linear differential equation with the Caputo fractional derivative of an order $\alpha \in (0, 1)$.
The goal of control is to minimize a given Bolza-type cost functional, which consists of two terms.
One of them evaluates the state vector of the system realized at a fixed terminal time $\vartheta$, and the other is an integral evaluation of a control on the whole time interval $[t_0, \vartheta]$.
We are interested in finding optimal, or, at least, $\varepsilon$-optimal, open-loop controls, as well as in constructing optimal feedback (closed-loop) controls, which are formalized within the framework of positional strategies \cite{Krasovskii_Subbotin_1988,Krasovskii_Krasovskii_1995} (see also \cite{Gomoyunov_2019_Trudy_Eng,Gomoyunov_2019_SIAM}).

In order to solve the problem, we propose an approach based on its reduction to some auxiliary optimal control problem in a dynamical system described by a first-order ordinary differential equation and further application of the methods and results of the optimal control theory widely developed for such systems.
The reduction relies on a suitable notion of a finite-dimensional informational image of an infinite-dimensional position \cite{Gomoyunov_2019_SIAM} of the original fractional-order system, which can be treated as a special prediction of a motion of this system at the time $\vartheta$.

This approach is closely related to a functional interpretation of control processes \cite{Krasovskii_1963} (see also \cite{Lukoyanov_Gomoyunov_2019_DGAA} and the references therein) and was previously developed for linear functional-differential systems of retarded \cite{Lukoyanov_Reshetova_1998_Eng,Gomoyunov_Plaksin_2015_IFAC} and neutral \cite{Gomoyunov_Lukoyanov_2018_Trudy_Eng} types and also for linear systems with control delays \cite{Gomoyunov_Lukoyanov_2012_Eng}.
However, in contrast to these studies, the auxiliary dynamical system obtained in the present paper may not satisfy the assumptions that are usually made in the optimal control theory.
More precisely, the right-hand side of the corresponding differential equation, in general, has a singularity at the time $\vartheta$ and, therefore, is unbounded.
This circumstance is explained by some special properties of the fundamental matrix solution of linear fractional-order differential equations (see, e.g., \cite{Idczak_Kamocki_2011,Bourdin_2018,Gomoyunov_2019_FCAA_2}).
In order to overcome this difficulty, we propose to introduce a small parameter $\eta \in (0, \vartheta - t_0)$ and shift the terminal time in the auxiliary optimal control problem from $\vartheta$ to $\vartheta_\eta = \vartheta - \eta$.
After such a modification, the auxiliary problem, on the one hand, can still be used to solve, at least approximately, the original problem, and, on the other hand, it already meets the typical assumptions from the optimal control theory.

The rest of the paper is organized as follows.
In Sect.~\ref{section_OCP}, we formulate the optimal control problem and study the question of finding optimal and $\varepsilon$-optimal open-loop controls.
We introduce the notions of a position of the system and its informational image, derive the auxiliary optimal control problem, and provide a connection between the original and auxiliary problems.
After that, we consider the auxiliary optimal control problem with the shifted terminal time $\vartheta_\eta$ and establish its connection with the original problem.
Sect.~\ref{section_OFC} is devoted to constructing optimal feedback controls.
We introduce the notion of a positional control strategy and show how to obtain an optimal positional control strategy on the basis of optimal positional strategies in the auxiliary problems with the shifted terminal times $\vartheta_\eta$.
Both sections contain illustrative examples.
The proofs of the statements are given in Appendix.

\section{Optimal Control Problem}
\label{section_OCP}

\subsection{Preliminaries}

Let $n \in \mathbb{N}$, $\alpha \in (0, 1)$, and $t_0$, $\vartheta \in \mathbb{R}$ such that $t_0 < \vartheta$ be fixed throughout the paper.
Let $\mathbb{R}^n$ and $\mathbb{R}^{n \times n}$ be the spaces of $n$-dimensional vectors and $(n \times n)$-matrices, and let $\operatorname{Id}_n \in \mathbb{R}^{n \times n}$ stand for the identity matrix.
By $\|\cdot\|$ and $\langle \cdot, \cdot \rangle$, we denote the Euclidean norm and the inner product in $\mathbb{R}^n$.
The corresponding norm in $\mathbb{R}^{n \times n}$ is also denoted by $\|\cdot\|$.

Let $t_\ast \in [t_0, \vartheta]$.
Let us consider a function $x: [t_0, t_\ast] \rightarrow \mathbb{R}^n$, for which we also use the notation $x(\cdot)$.
The (left-sided) Riemann--Liouville fractional integral $(I^\alpha x)(t)$ and Caputo fractional derivative $(^C D^\alpha x)(t)$ of the order $\alpha$ of $x(\cdot)$ at $t \in [t_0, t_\ast]$ are defined respectively by
\begin{equation*}
    \begin{array}{c}
        \displaystyle
        (I^\alpha x) (t)
        = \frac{1}{\Gamma(\alpha)} \int_{t_0}^{t} \frac{x(\tau)}{(t - \tau)^{1 - \alpha}} \, \mathrm{d} \tau, \\[0.5em]
        \displaystyle
        (^C D^\alpha x)(t)
        = \frac{\mathrm{d}}{\mathrm{d} t} \big(I^{1 - \alpha} (x(\cdot) - x(t_0)) \big) (t)
        = \frac{1}{\Gamma(1 - \alpha)} \frac{\mathrm{d}}{\mathrm{d} t} \int_{t_0}^{t} \frac{x(\tau) - x(t_0)}{(t - \tau)^\alpha} \, \mathrm{d} \tau,
    \end{array}
\end{equation*}
where $\Gamma$ is the gamma-function.
The basic properties of integrals and derivatives of fractional order can be found, e.g., in \cite{Samko_Kilbas_Marichev_1993,Kilbas_Srivastava_Trujillo_2006,Diethelm_2010}.
We say that $x(\cdot) \in AC^\alpha([t_0, t_\ast], \mathbb{R}^n)$, if there exists a (Lebesgue) measurable and essentially bounded function $\varphi: [t_0, t_\ast] \rightarrow \mathbb{R}^n$ such that $x(t) = x(t_0) + (I^\alpha \varphi) (t)$, $t \in [t_0, t_\ast]$.
Let us note that, in the case when $t_\ast = t_0$, the set $AC^\alpha([t_0, t_\ast], \mathbb{R}^n)$ can be identified with $\mathbb{R}^n$.

\subsection{Statement of the Problem}
\label{subsection_OCP_Statement_of_the_problem}

We consider a dynamical system which motion is described by the linear fractional differential equation
\begin{equation} \label{system}
    (^C D^\alpha x)(t) = A(t) x(t) + f(t, u(t)),
    \quad x(t) \in \mathbb{R}^n, \quad u(t) \in \mathbb{U}, \quad t \in [t_0, \vartheta].
\end{equation}
Here, $t$ is the time;
$x(t)$ and $u(t)$ are respectively the current values of the state and control vectors;
the set $\mathbb{U} \subset \mathbb{R}^r$ is compact, $r \in \mathbb{N}$;
$t_0$ and $\vartheta$ are the initial and terminal times.
We assume that the right-hand side of the differential equation in (\ref{system}) satisfies the following conditions:
\begin{description}
    \item[$(A.1)$]
        The function $A: [t_0, \vartheta] \rightarrow \mathbb{R}^{n \times n}$ is measurable and essentially bounded.

    \item[$(A.2)$]
        The function $f: [t_0, \vartheta] \times \mathbb{U} \rightarrow \mathbb{R}^n$ is continuous.
\end{description}

Let $R_x > 0$ be fixed.
We suppose that, at the initial time $t_0$, an initial value $x_0 \in B(R_x) = \{y \in \mathbb{R}^n: \|y\| \leq R_x\}$ of the state vector of system (\ref{system}) is given.
As the set $\mathcal{U}(t_0, \vartheta)$ of admissible controls $u(\cdot)$ on the time interval $[t_0, \vartheta)$, we consider the set of all measurable functions $u: [t_0, \vartheta) \rightarrow \mathbb{U}$.
By a motion of the system that corresponds to the initial value $x_0$ and a control $u(\cdot) \in \mathcal{U}(t_0, \vartheta)$, we mean a function $x(\cdot) \in AC^\alpha([t_0, \vartheta], \mathbb{R}^n)$ that satisfies the initial condition
\begin{equation} \label{initial_condition}
    x(t_0)
    = x_0
\end{equation}
and, together with $u(\cdot)$, satisfies the differential equation in (\ref{system}) for almost every $t \in [t_0, \vartheta]$.
Due to conditions $(A.1)$ and $(A.2)$, such a motion $x(\cdot)$ exists and is unique (see, e.g., \cite[Theorem~3.1]{Gomoyunov_2019_FCAA}), and we denote it by $x(\cdot \mid t_0, x_0, \vartheta, u(\cdot))$.

The goal of control is to minimize the cost functional
\begin{equation} \label{cost_functional}
    J(t_0, x_0, u(\cdot))
    = \sigma(x(\vartheta))
    + \int_{t_0}^{\vartheta} \chi(t, u(t)) \, \mathrm{d} \tau,
    \quad u(\cdot) \in \mathcal{U}(t_0, \vartheta),
\end{equation}
where $x(\cdot) = x(\cdot \mid t_0, x_0, \vartheta, u(\cdot))$.
We assume that the conditions below hold:
\begin{description}
    \item[$(A.3)$]
        The function $\sigma: \mathbb{R}^n \rightarrow \mathbb{R}$ is continuous.

    \item[$(A.4)$]
        The function $\chi: [t_0, \vartheta] \times \mathbb{U} \rightarrow \mathbb{R}$ is continuous.
\end{description}

The value of the optimal result in the control problem for system (\ref{system}) with initial condition (\ref{initial_condition}) and cost functional (\ref{cost_functional}) is defined by
\begin{equation} \label{Val}
    \rho(t_0, x_0)
    = \inf_{u(\cdot) \in \mathcal{U}(t_0, \vartheta)} J(t_0, x_0, u(\cdot)).
\end{equation}
A control $u^\circ(\cdot) \in \mathcal{U}(t_0, \vartheta)$ is called optimal in this problem, if the equality $J(t_0, x_0, u^\circ(\cdot)) = \rho(t_0, x_0)$ holds.

In order to find such an optimal control $u^\circ(\cdot)$, we propose to reduce the problem (\ref{system})--(\ref{cost_functional}) to some auxiliary optimal control problem in a dynamical system which motion is described by a first-order ordinary differential equation.
This reduction is based on some linear transformation, which is called the information image of a position of system (\ref{system}).

\subsection{Positions of the System}

According to \cite{Gomoyunov_2019_SIAM}, by a position of system (\ref{system}), we mean a pair $(t, w(\cdot))$ consisting of a time $t \in [t_0, \vartheta]$ and a function $w(\cdot) \in AC^\alpha([t_0, t], \mathbb{R}^n)$, $\|w(t_0)\| \leq R_x$, which is treated as a history of a motion of the system on the time interval $[t_0, t]$.
The set of all such positions is denoted by $G$.
Respectively, for every $x_0 \in B(R_x)$, the pair $(t_0, x_0) \in G$ is regarded as an initial position.

By analogy with Sect.~\ref{subsection_OCP_Statement_of_the_problem}, let us give a definition of motions of the system starting from an arbitrary position $(t_\ast, w_\ast(\cdot)) \in G$ and evolving on some time interval $[t_\ast, t^\ast]$, where $t^\ast \in [t_\ast, \vartheta]$.
Let us suppose that $t_\ast < t^\ast$.
Then, as the set $\mathcal{U}(t_\ast, t^\ast)$ of admissible controls $u(\cdot)$ on $[t_\ast, t^\ast)$, we consider the set of all measurable functions $u: [t_\ast, t^\ast) \rightarrow \mathbb{U}$.
A motion of the system generated from the position $(t_\ast, w_\ast(\cdot))$ by a control $u(\cdot) \in \mathcal{U}(t_\ast, t^\ast)$ is defined as a function $x(\cdot) \in AC^\alpha([t_0, t^\ast], \mathbb{R}^n)$ that satisfies the equality
\begin{equation} \label{w_ast}
    x(t)
    = w_\ast(t),
    \quad t \in [t_0, t_\ast],
\end{equation}
and, together with $u(\cdot)$, satisfies the differential equation in (\ref{system}) for almost every $t \in [t_\ast, t^\ast]$.
By the scheme from \cite[Proposition~2]{Gomoyunov_2019_DGAA}, one can prove that, owing to conditions $(A.1)$ and $(A.2)$, such a motion $x(\cdot) = x(\cdot \mid t_\ast, w_\ast(\cdot), t^\ast, u(\cdot))$ exists and is unique.
In the degenerate case when $t_\ast = t^\ast$, the motion $x(\cdot)$ is completely determined by (\ref{w_ast}), and there is no
need in considering controls $u(\cdot)$ and determining the set $\mathcal{U}(t_\ast, t^\ast)$.
However, it is convenient to formally say that this motion $x(\cdot) = x(\cdot \mid t_\ast, w_\ast(\cdot), t^\ast, u(\cdot))$ is generated from $(t_\ast, w_\ast(\cdot))$ by $u(\cdot) \in \mathcal{U}(t_\ast, t^\ast)$.
Further, for the motion $x(\cdot)$ and a time $t \in [t_0, t^\ast]$, we denote the corresponding position of the system by $(t, x_t(\cdot))$, where the function $x_t: [t_0, t] \rightarrow \mathbb{R}^n$ is defined by
\begin{equation} \label{x_t}
    x_t(\tau)
    = x(\tau),
    \quad \tau \in [t_0, t].
\end{equation}
Let us note that the inclusion $(t, x_t(\cdot)) \in G$ is valid.

Following \cite{Gomoyunov_2019_FCAA_2} (see also \cite{Bourdin_2018}), let us consider the fundamental solution matrix of the differential equation in (\ref{system}), which is a continuous function
\begin{equation} \label{Omega}
    \Omega \ni (t, \tau) \mapsto F(t, \tau) \in \mathbb{R}^{n \times n},
    \quad \Omega
    = \{ (t, \tau) \in [t_0, \vartheta] \times [t_0, \vartheta]: t \geq \tau\},
\end{equation}
such that, for every fixed $\tau \in [t_0, \vartheta]$, the function $[\tau, \vartheta] \ni t \mapsto F(t, \tau) \in \mathbb{R}^{n \times n}$ is a unique continuous solution to the integral equation
\begin{equation*}
    F(t, \tau)
    = \frac{\operatorname{Id}_n}{\Gamma(\alpha)}
    + \frac{(t - \tau)^{1 - \alpha}}{\Gamma(\alpha)} \int_{\tau}^{t} \frac{A(\xi) F(\xi, \tau)}{(t - \xi)^{1 - \alpha} (\xi - \tau)^{1 - \alpha}} \, \mathrm{d} \xi,
    \quad t \in [\tau, \vartheta].
\end{equation*}
Due to \cite[Theorem~5.2]{Gomoyunov_2019_FCAA_2}, for every motion $x(\cdot) = x(\cdot \mid t_\ast, w_\ast(\cdot), t^\ast, u(\cdot))$ of system (\ref{system}), where $(t_\ast, w_\ast(\cdot)) \in G$, $t^\ast \in [t_\ast, \vartheta]$, and $u(\cdot) \in \mathcal{U}(t_\ast, t^\ast)$, the representation formula below holds:
\begin{multline} \label{Cauchy_formula}
    x(t)
    = \Big( \operatorname{Id}_n + \int_{t_\ast}^{t} \frac{F(t, \tau) A(\tau)}{(t - \tau)^{1 - \alpha}} \, \mathrm{d} \tau \Big) w_\ast(t_\ast) \\
    + \frac{1}{\Gamma(1 - \alpha)}\int_{t_\ast}^{t} \frac{F(t, \tau)}{(t - \tau)^{1 - \alpha}}
    \Big( \alpha \int_{t_0}^{t_\ast} \frac{w_\ast(\xi) - w_\ast(t_0)}{(\tau - \xi)^{1 + \alpha}} \, \mathrm{d} \xi
    - \frac{w_\ast(t_\ast) - w_\ast(t_0)}{(\tau - t_\ast)^\alpha} \Big) \, \mathrm{d} \tau \\
    + \int_{t_\ast}^{t} \frac{F(t, \tau) f(\tau, u(\tau))}{(t - \tau)^{1 - \alpha}} \, \mathrm{d} \tau,
    \quad t \in [t_\ast, t^\ast].
\end{multline}

\subsection{Informational Image}

For a position $(t_\ast, w_\ast(\cdot)) \in G$, let us define the value $\mathcal{I}(t_\ast, w_\ast(\cdot)) \in \mathbb{R}^n$, called the informational image of $(t_\ast, w_\ast(\cdot))$, as follows.
Let us consider the linear homogenous fractional differential equation corresponding to (\ref{system})
\begin{equation} \label{system_homogeneous}
    (^C D^\alpha y)(t) = A(t) y(t), \quad y(t) \in \mathbb{R}^n, \quad t \in [t_\ast, \vartheta],
\end{equation}
with the initial condition
\begin{equation} \label{initial_condition_homogenous}
    y(t)
    = w_\ast(t), \quad t \in [t_0, t_\ast].
\end{equation}
By analogy with the above, there exists a unique solution to the Cauchy problem (\ref{system_homogeneous}) and (\ref{initial_condition_homogenous}), which is the function $y(\cdot) = y(\cdot \mid t_\ast, w_\ast(\cdot), \vartheta) \in AC^\alpha([t_0, \vartheta], \mathbb{R}^n)$ that satisfies the equality in (\ref{initial_condition_homogenous}) and the differential equation in (\ref{system_homogeneous}) for almost every $t \in [t_\ast, \vartheta]$.
Then, we put
\begin{equation} \label{informational_image}
    \mathcal{I}(t_\ast, w_\ast(\cdot))
    = y(\vartheta \mid t_\ast, w_\ast(\cdot), \vartheta).
\end{equation}
Let us note that, owing to representation formula (\ref{Cauchy_formula}), the informational image $\mathcal{I}(t_\ast, w_\ast(\cdot))$ can be defined explicitly:
\begin{multline} \label{informational_image_Cauchy}
    \mathcal{I}(t_\ast, w_\ast(\cdot))
    = \Big( \operatorname{Id}_n + \int_{t_\ast}^{\vartheta} \frac{F(\vartheta, \tau) A(\tau)}{(\vartheta - \tau)^{1 - \alpha}} \, \mathrm{d} \tau \Big) w_\ast(t_\ast) \\
    + \frac{1}{\Gamma(1 - \alpha)}\int_{t_\ast}^{\vartheta} \frac{F(\vartheta, \tau)}{(\vartheta - \tau)^{1 - \alpha}}
    \Big( \alpha \int_{t_0}^{t_\ast} \frac{w_\ast(\xi) - w_\ast(t_0)}{(\tau - \xi)^{1 + \alpha}} \, \mathrm{d} \xi
    - \frac{w_\ast(t_\ast) - w_\ast(t_0)}{(\tau - t_\ast)^\alpha} \Big) \, \mathrm{d} \tau.
\end{multline}

Taking the last term in (\ref{Cauchy_formula}) into account, let us introduce the function
\begin{equation} \label{f_ast}
    f_\ast(t, u)
    = \frac{F(\vartheta, t) f(t, u)}{(\vartheta - t)^{1 - \alpha}},
    \quad t \in [t_0, \vartheta), \quad u \in \mathbb{U}.
\end{equation}
Due to continuity of $F$ and $f$, the function $f_\ast: [t_0, \vartheta) \times \mathbb{U} \rightarrow \mathbb{R}^n$ is continuous, and the following estimate is valid:
\begin{equation} \label{f_ast_bound}
    \|f_\ast(t, u)\|
    \leq \frac{M_F M_f}{(\vartheta - t)^{1 - \alpha}},
    \quad t \in [t_0, \vartheta), \quad u \in \mathbb{U},
\end{equation}
where we denote
\begin{equation}\label{M_F_M_f}
    M_F = \max_{t \in [t_0, \vartheta]} \|F(\vartheta, t)\|,
    \quad M_f = \max_{(t, u) \in [t_0, \vartheta] \times \mathbb{U}} \|f(t, u)\|.
\end{equation}

The proposition below describes the dynamics of informational image (\ref{informational_image}) along motions of system (\ref{system}).
\begin{proposition} \label{Proposition_Informational_image}
    Let $x_0 \in B(R_x)$, $t^\ast \in [t_0, \vartheta]$, and $u(\cdot) \in \mathcal{U}(t_0, t^\ast)$.
    Let $x(\cdot) = x(\cdot \mid t_0, x_0, t^\ast, u(\cdot))$ be the corresponding motion of system $(\ref{system})$.
    Then, the equality below holds:
    \begin{equation} \label{Proposition_Informational_image_main}
        \mathcal{I}(t, x_t(\cdot))
        = \mathcal{I}(t_0, x_0)
        + \int_{t_0}^{t} f_\ast(\tau, u(\tau)) \, \mathrm{d} \tau,
        \quad t \in [t_0, t^\ast],
    \end{equation}
    where $x_t(\cdot)$ is defined by $x(\cdot)$ according to $(\ref{x_t})$.
\end{proposition}

The proof of the proposition is given in Appendix.

\subsection{Auxiliary Optimal Control Problem}

Based on Proposition~\ref{Proposition_Informational_image}, we consider the auxiliary dynamical system which motion is described by the differential equation
\begin{equation} \label{system_z}
    \dot{z}(t) = f_\ast(t, p(t)),
    \quad z(t) \in \mathbb{R}^n, \quad p(t) \in \mathbb{U}, \quad t \in [t_0, \vartheta),
\end{equation}
with the initial condition
\begin{equation} \label{initial_condition_z}
    z(t_0)
    = z_0
    = \mathcal{I}(t_0, x_0).
\end{equation}
Here, $\dot{z}(t) = \mathrm{d} z / \mathrm{d} t$;
$z(t)$ and $p(t)$ are respectively the current values of the state and control vectors in the auxiliary system;
the set $\mathbb{U}$ is the same as in original system (\ref{system});
the function $f_\ast$ is defined in (\ref{f_ast});
an initial value $z_0$ is determined by the informational image $\mathcal{I}(t_0, x_0)$ (see (\ref{informational_image})) of an initial position $(t_0, x_0) \in G$ of system (\ref{system}).

By a position of auxiliary system (\ref{system_z}), we mean a pair $(t, z) \in [t_0, \vartheta] \times \mathbb{R}^n$.
Let $(t_\ast, z_\ast) \in [t_0, \vartheta] \times \mathbb{R}^n$, $t^\ast \in [t_\ast, \vartheta]$, and $p(\cdot) \in \mathcal{U}(t_\ast, t^\ast)$.
A motion $z(\cdot)$ of this system generated from $(t_\ast, z_\ast)$ by $p(\cdot)$ is an absolutely continuous function $z: [t_\ast, t^\ast] \rightarrow \mathbb{R}^n$ that satisfies the equality $z(t_\ast) = z_\ast$ and, together with $p(\cdot)$, satisfies the differential equation in (\ref{system_z}) for almost every $t \in [t_\ast, t^\ast]$.
In view of the described above properties of the function $f_\ast$, such a motion exists and is unique, and we denote it by $z(\cdot \mid t_\ast, z_\ast, t^\ast, p(\cdot))$.

As a direct consequence of Proposition~\ref{Proposition_Informational_image}, we derive the following lemma, giving a connection between motions of the original and auxiliary systems.
\begin{lemma} \label{Lemma_Connection}
    Let $x_0 \in B(R_x)$, $t^\ast \in [t_0, \vartheta]$, and $u(\cdot) \in \mathcal{U}(t_0, t^\ast)$.
    Let $x(\cdot) = x(\cdot \mid t_0, x_0, t^\ast, u(\cdot))$ be the corresponding motion of system $(\ref{system})$ with initial condition $(\ref{initial_condition})$.
    Let $z(\cdot) = z(\cdot \mid t_0, z_0, t^\ast, p(\cdot))$ be the motion of auxiliary system $(\ref{system_z})$ with initial condition $(\ref{initial_condition_z})$ generated by the same control $p(t) = u(t)$, $t \in [t_0, t^\ast)$.
    Then, the equality below holds:
    \begin{equation*}
        \mathcal{I}(t, x_t(\cdot))
        = z(t),
        \quad t \in [t_0, t^\ast].
    \end{equation*}
\end{lemma}

Further, taking into account that (see (\ref{informational_image}))
\begin{equation} \label{informational_image_terminal}
    \mathcal{I}(\vartheta, w_\ast(\cdot))
    = w_\ast(\vartheta),
    \quad (\vartheta, w_\ast(\cdot)) \in G,
\end{equation}
we define the auxiliary cost functional to be minimized by $p(\cdot)$ as follows:
\begin{equation} \label{cost_functional_z}
    J_\ast(t_0, z_0, p(\cdot))
    = \sigma(z(\vartheta))
    + \int_{t_0}^{\vartheta} \chi(\tau, p(\tau)) \, \mathrm{d} \tau,
    \quad p(\cdot) \in \mathcal{U}(t_0, \vartheta),
\end{equation}
where $z(\cdot) = z(\cdot \mid t_0, z_0, \vartheta, p(\cdot))$, and the functions $\sigma$ and $\chi$ are taken from (\ref{cost_functional}).

The value of the optimal result in the auxiliary control problem for system (\ref{system_z}) with initial condition (\ref{initial_condition_z}) and cost functional (\ref{cost_functional_z}) is given by
\begin{equation*}
    \rho_\ast(t_0, z_0)
    = \inf_{p(\cdot) \in \mathcal{U}(t_0, \vartheta)} J_\ast(t_0, z_0, p(\cdot)),
\end{equation*}
and a control $p^\circ(\cdot) \in \mathcal{U}(t_0, \vartheta)$ is called optimal if $J_\ast(t_0, z_0, p^\circ(\cdot)) = \rho_\ast(t_0, z_0)$.

\begin{theorem} \label{Theorem_OCP}
    For any $x_0 \in B(R_x)$, the following statements are valid:
    \begin{description}
        \item[$i)$]
            A control $u^\circ(\cdot) \in \mathcal{U}(t_0, \vartheta)$ is optimal in the problem $(\ref{system})$--$(\ref{cost_functional})$ if and only if it is optimal in the auxiliary problem $(\ref{system_z})$, $(\ref{initial_condition_z})$, and $(\ref{cost_functional_z})$.
        \item[$ii)$]
            The optimal results in the original and auxiliary problems coincide, i.e., $\rho(t_0, x_0) = \rho_\ast(t_0, z_0)$.
    \end{description}
\end{theorem}

The proof of the theorem is given in Appendix.

Thus, the original optimal control problem in fractional-order system (\ref{system}) is reduced to the auxiliary optimal control problem in first-order system (\ref{system_z}).

\begin{remark} \label{Remark_Reduction}
    Let us assume that the function $\sigma$ from cost functional (\ref{cost_functional}) can be represented in the following form:
    \begin{equation} \label{Remark_Reduction_sigma}
        \sigma(x)
        = \mu( K (x - c)),
        \quad x \in \mathbb{R}^n.
    \end{equation}
    Here, $K$ is a $(n_\ast \times n)$-matrix, $n_\ast \in \mathbb{N}$, $n_\ast < n$;
    $c \in \mathbb{R}^n$;
    $\mu: \mathbb{R}^{n_\ast} \rightarrow \mathbb{R}$ is a continuous function.
    In particular, this is the case when $\sigma$ does not depend on some $n - n_\ast$ coordinates of $x$.
    Under this additional assumption, we can reduce the dimension of the state vector in the auxiliary system from $n$ to $n_\ast$ and simplify the auxiliary cost functional.
    Namely, we consider the system
    \begin{equation*}
        \dot{z}(t) = K f_\ast(t, p(t)),
        \quad z(t) \in \mathbb{R}^{n_\ast}, \quad p(t) \in \mathbb{U}, \quad t \in [t_0, \vartheta),
    \end{equation*}
    with the initial condition $z(t_0) = z_0 = K (\mathcal{I}(t_0, x_0) - c)$ and the cost functional
    \begin{equation*}
        J_\ast(t_0, z_0, p(\cdot))
        = \mu(z(\vartheta))
        + \int_{t_0}^{\vartheta} \chi(t, p(t)) \, \mathrm{d} t,
        \quad p(\cdot) \in \mathcal{U}(t_0, \vartheta).
    \end{equation*}
    One can show that, for this auxiliary optimal control problem, the result similar to Theorem~\ref{Theorem_OCP} takes place.
\end{remark}

However, it should be noted that, in general, the right-hand side of the differential equation in (\ref{system_z}) is unbounded, and, therefore, the auxiliary problem does not satisfy the assumptions that are usually made in the optimal control theory.
This complicates the application of the methods and results developed within this theory to solving the original optimal control problem.
In order to overcome this difficulty, in Sect.~\ref{subsection_OCP_Auxiliary_Optimal_Control_Problem_with_Parameter}, we propose to introduce a small parameter $\eta \in (0, \vartheta - t_0)$ and consider auxiliary system (\ref{system_z}) only up to the shifted terminal time $\vartheta_\eta = \vartheta - \eta$.
But before doing this, let us give an example illustrating Theorem~\ref{Theorem_OCP} and Remark~\ref{Remark_Reduction}.

\begin{example}
    Following \cite[Sect.~4.2]{Bergounioux_Bourdin_2019}, let us consider the optimal control problem described by the system
    \begin{equation} \label{Example_1_system}
        \begin{array}{c}
            \left\{
                \begin{array}{l}
                    (^C D^\alpha x_1)(t) = x_2(t) + \cos(u(t)), \\[0.5em]
                    (^C D^\alpha x_2)(t) = \sin(u(t)),
                \end{array}
            \right. \\[1.5em]
            x(t) = (x_1(t), x_2(t)) \in \mathbb{R}^2,
            \quad u(t) \in [- \pi / 2, \pi / 2],
            \quad t \in [t_0, \vartheta],
        \end{array}
    \end{equation}
    with the initial condition $x(t_0) = x_0 = (0, 0)$ and the cost functional
    \begin{equation*}
        J(t_0, x_0, u(\cdot))
        = - x_1(\vartheta),
        \quad u(\cdot) \in \mathcal{U}(t_0, \vartheta).
    \end{equation*}

    According to \cite[Theorem~4.2]{Idczak_Kamocki_2011} (see also \cite{Bourdin_2018,Gomoyunov_2019_FCAA_2}), the fundamental solution matrix $F$ of the differential equation in (\ref{Example_1_system}) is given by
    \begin{equation} \label{Example_1_F}
        F(t, \tau)
        = \sum_{i = 0}^\infty \frac{(t - \tau)^{i \alpha} A^i}{\Gamma((i + 1) \alpha)}
        = \begin{pmatrix}
            \frac{1}{\Gamma(\alpha)} & \frac{(t - \tau)^\alpha}{\Gamma(2 \alpha)} \\
            0 & \frac{1}{\Gamma(\alpha)}
        \end{pmatrix},
        \quad A
        = \begin{pmatrix}
            0 & 1 \\
            0 & 0
        \end{pmatrix},
    \end{equation}
    where $(t, \tau) \in \Omega$ (see (\ref{Omega})).
    Further, due to (\ref{informational_image_Cauchy}), we have
    \begin{equation} \label{Example_1_I}
        \mathcal{I}(t_0, x_0) = (0, 0).
    \end{equation}
    Thus, denoting
    \begin{equation} \label{Example_1_b}
        b_1(t)
        = \frac{1}{\Gamma(\alpha) (\vartheta - t)^{1 - \alpha}},
        \quad
        b_2(t)
        = \frac{1}{\Gamma(2 \alpha) (\vartheta - t)^{1 - 2 \alpha}},
        \quad t \in [t_0, \vartheta),
    \end{equation}
    we come to the auxiliary optimal control problem for the system
    \begin{equation*}
        \begin{array}{c}
            \displaystyle
            \dot{z}(t)
            = b_1(t) \cos(p(t)) + b_2(t) \sin(p(t)), \\[0.5em]
            \displaystyle
            z(t) \in \mathbb{R}, \quad p(t) \in [- \pi / 2, \pi / 2], \quad t \in [t_0, \vartheta),
        \end{array}
    \end{equation*}
    with the initial condition $z(t_0) = z_0 = 0$ and the cost functional
    \begin{equation*}
        J_\ast(t_0, z_0, p(\cdot))
        = - z(\vartheta),
        \quad p(\cdot) \in \mathcal{U}(t_0, \vartheta).
    \end{equation*}

    By direct calculations, we obtain
    \begin{multline*}
        \rho_\ast(t_0, z_0)
        = - \sup_{p(\cdot) \in \mathcal{U}(t_0, \vartheta)} \int_{t_0}^{\vartheta}
        \big( b_1(t) \cos(p(t)) + b_2(t) \sin(p(t)) \big) \, \mathrm{d} t \\
        = - \int_{t_0}^{\vartheta}
        \big( b_1(t) \cos(p^\circ(t)) + b_2(t) \sin(p^\circ(t)) \big) \, \mathrm{d} t
    \end{multline*}
    for the control $p^\circ(t) = \arctan( \Gamma(\alpha) (\vartheta - t)^\alpha / \Gamma(2 \alpha))$, $t \in [t_0, \vartheta]$.
    Then, this control $p^\circ(\cdot)$ is optimal in the auxiliary problem.
    Hence, by Theorem~\ref{Theorem_OCP} and Remark~\ref{Remark_Reduction}, we conclude that $p^\circ(\cdot)$ is an optimal control in the original problem, too.
    Let us note that this result agrees with \cite[Sect.~4.2]{Bergounioux_Bourdin_2019}.
\end{example}

\subsection{Auxiliary Optimal Control Problem with Parameter}
\label{subsection_OCP_Auxiliary_Optimal_Control_Problem_with_Parameter}

Let us fix $\eta \in (0, \vartheta - t_0)$, put $\vartheta_\eta = \vartheta - \eta \in (t_0, \vartheta)$, and consider the auxiliary optimal control problem for the system
\begin{equation} \label{system_z_eta}
    \dot{z}(t) = f_\ast(t, p_\eta(t)),
    \quad z(t) \in \mathbb{R}^n, \quad p_\eta(t) \in \mathbb{U}, \quad t \in [t_0, \vartheta_\eta],
\end{equation}
with the initial condition
\begin{equation} \label{initial_condition_z_eta}
    z(t_0)
    = z_0
    = \mathcal{I}(t_0, x_0)
\end{equation}
and the cost functional
\begin{equation} \label{cost_functional_z_eta}
    J_\eta(t_0, z_0, p_\eta(\cdot))
    = \sigma(z(\vartheta_\eta))
    + \int_{t_0}^{\vartheta_\eta} \chi (t, p_\eta(t)) \, \mathrm{d} t,
    \quad p_\eta(\cdot) \in \mathcal{U}(t_0, \vartheta_\eta).
\end{equation}
The value of the optimal result in this problem is defined by
\begin{equation} \label{Valaeta}
    \rho_\eta(t_0, z_0)
    = \inf_{p_\eta(\cdot) \in \mathcal{U}(t_0, \vartheta_\eta)} J_\eta(t_0, z_0, p_\eta(\cdot)),
\end{equation}
and a control $p^\circ_\eta(\cdot) \in \mathcal{U}(t_0, \vartheta_\eta)$ is called optimal if $J_\eta(t_0, z_0, p_\eta^\circ(\cdot)) = \rho_\eta(t_0, z_0)$.
Thus, the only difference from the auxiliary optimal control problem (\ref{system_z}), (\ref{initial_condition_z}), and (\ref{cost_functional_z}) is that now the time $\vartheta_\eta$ is treated as the terminal one.

Let us note that, since the function (see (\ref{f_ast}))
\begin{equation} \label{f_ast_continuity}
    [t_0, \vartheta_\eta] \times \mathbb{U} \ni (t, u) \mapsto f_\ast(t, u) \in \mathbb{R}^n
\end{equation}
is continuous (and, therefore, bounded), we obtain that, compared to system (\ref{system_z}), system (\ref{system_z_eta}) meets the typical assumptions from the optimal control theory.
On the other hand, in contrast to Theorem~\ref{Theorem_OCP}, due to the presence of the parameter $\eta$, one can not expect that an optimal control $p^\circ_\eta(\cdot)$ in the auxiliary problem (\ref{system_z_eta})--(\ref{cost_functional_z_eta}) determines some optimal control $u^\circ(\cdot)$ in the original problem (\ref{system})--(\ref{cost_functional}).

In this connection, let us consider a notion of $\varepsilon$-optimal controls.
Namely, for a number $\varepsilon > 0$, a control $u^\ast(\cdot) \in \mathcal{U}(t_0, \vartheta)$ is called $\varepsilon$-optimal in the original problem if $J(t_0, x_0, u^\ast(\cdot)) \leq \rho(t_0, x_0) + \varepsilon$.
Respectively, $p^\ast_\eta(\cdot) \in \mathcal{U}(t_0, \vartheta_\eta)$ is $\varepsilon$-optimal in the auxiliary problem if $J_\eta(t_0, z_0, p^\ast_\eta(\cdot)) \leq \rho_\eta(t_0, z_0) + \varepsilon$.

\begin{theorem} \label{Theorem_OCP_varepsilon}
    For any $\varepsilon > 0$, there exist $\eta^\ast = \eta^\ast(\varepsilon) \in (0, \vartheta - t_0)$ and $\varepsilon^\ast = \varepsilon^\ast(\varepsilon) > 0$ such that, for any $\eta \in (0, \eta^\ast]$ and any $x_0 \in B(R_x)$, the following statements hold:
    \begin{description}
        \item[$i)$]
            If a control $p^\ast_\eta(\cdot)$ is $\varepsilon^\ast$-optimal in the auxiliary problem $(\ref{system_z_eta})$--$(\ref{cost_functional_z_eta})$ corresponding to the chosen $\eta$, then the control
            \begin{equation} \label{Theorem_OCP_varepsilon_main}
                u^\ast(t)
                = \begin{cases}
                    p^\ast_\eta(t), & \mbox{if } t \in [t_0, \vartheta_\eta), \\
                    \bar{u}, & \mbox{if } t \in (\vartheta_\eta, \vartheta),
                  \end{cases}
                \quad \bar{u} \in \mathbb{U},
            \end{equation}
            is $\varepsilon$-optimal in the problem $(\ref{system})$--$(\ref{cost_functional})$.
        \item[$ii)$]
            The optimal results in the original and auxiliary problems satisfy the inequality $|\rho(t_0, x_0) - \rho_\eta(t_0, z_0)| \leq \varepsilon$.
    \end{description}
\end{theorem}

The proof of the theorem is given in Appendix.

Theorem~\ref{Theorem_OCP_varepsilon} allows us to apply, via the auxiliary problem (\ref{system_z_eta})--(\ref{cost_functional_z_eta}), the methods and results of the optimal control theory for first-order systems to finding $\varepsilon$-optimal controls in the problem (\ref{system})--(\ref{cost_functional}).
Let us note that a remark similar to Remark~\ref{Remark_Reduction} can also be made in relation to Theorem~\ref{Theorem_OCP_varepsilon}.
Let us consider an example.

\begin{example}
    Let the optimal control problem be described by the system
    \begin{equation} \label{Example_2_system}
        \begin{array}{c}
            \left\{
                \begin{array}{l}
                    (^C D^\alpha x_1)(t) = x_2(t), \\[0.5em]
                    (^C D^\alpha x_2)(t) = u(t),
                \end{array}
            \right. \\[1.5em]
            x(t) = (x_1(t), x_2(t)) \in \mathbb{R}^2,
            \quad u(t) \in [- 1, 1],
            \quad t \in [t_0, \vartheta],
        \end{array}
    \end{equation}
    with the initial condition $x(t_0) = x_0 = (0, 0)$ and the cost functional
    \begin{equation*}
        J(t_0, x_0, u(\cdot))
        = (x_1(\vartheta) - c_1)^2 + \int_{t_0}^{\vartheta} u^2(t) \, \mathrm{d} t,
        \quad u(\cdot) \in \mathcal{U}(t_0, \vartheta),
    \end{equation*}
    where $c_1 \in \mathbb{R}$ is a given number.

    For every $\eta \in (0, \vartheta - t_0)$, taking (\ref{Example_1_F})--(\ref{Example_1_b}) into account, we come to the auxiliary optimal control problem for the system
    \begin{equation} \label{Example_2_system_auxiliary}
        \dot{z}(t)
        = b_2(t) p_\eta(t),
        \quad z(t) \in \mathbb{R}, \quad p_\eta(t) \in [- 1, 1], \quad t \in [t_0, \vartheta_\eta],
    \end{equation}
    with the initial condition $z(t_0) = z_0 = - c_1$ and the cost functional
    \begin{equation*}
        J_\eta(t_0, z_0, p_\eta(\cdot))
        = z^2(\vartheta_\eta) + \int_{t_0}^{\vartheta_\eta} p_\eta^2(t) \, \mathrm{d} t,
        \quad p_\eta(\cdot) \in \mathcal{U}(t_0, \vartheta_\eta).
    \end{equation*}
    Applying the Pontryagin maximum principle to this auxiliary problem, we obtain (see, e.g., \cite[Ch.~3, Theorem~14]{Lee_Markus_1967} and also \cite[Statement~1]{Shaburov_2017}) that the unique optimal control is given by $p^\circ_\eta(t) = b_2(t) \lambda_\eta / S(b_2(t) |\lambda_\eta|)$, $t \in [t_0, \vartheta_\eta)$, where $S(\zeta) = 2$ for $\zeta \in [0, 2]$ and $S(\zeta) = \zeta$ for $\zeta > 2$, and $\lambda_\eta \in \mathbb{R}$ is the unique solution to the equation
    \begin{equation*}
        \lambda \int_{t_0}^{\vartheta_\eta} \frac{b_2^{\, 2}(t)}{S(b_2(t) |\lambda|)} \, \mathrm{d} t + \frac{\lambda}{2}
        = c_1.
    \end{equation*}

    Consequently, by Theorem~\ref{Theorem_OCP_varepsilon} and Remark~\ref{Remark_Reduction}, for every $\varepsilon > 0$, on the basis of the control $p^\circ_\eta(\cdot)$ for a sufficiently small $\eta \in (0, \vartheta - t_0)$, we can determine an $\varepsilon$-optimal in the original problem control $u^\ast(\cdot)$ according to (\ref{Theorem_OCP_varepsilon_main}).
\end{example}

\section{Optimal Feedback Controls}
\label{section_OFC}

In the previous section, we study the question of finding optimal and $\varepsilon$-optimal (open-loop) controls $u(\cdot) \in \mathcal{U}(t_0, \vartheta)$ in the optimal control problem (\ref{system})--(\ref{cost_functional}).
However, it is often more convenient to use feedback (closed-loop) controls, since, in many cases, they are easier to construct, and, moreover, they do not depend on a particular choice of an initial value $x_0 \in B(R_x)$.
In this section, we develop the proposed above reduction of the original problem to the auxiliary problem (\ref{system_z_eta})--(\ref{cost_functional_z_eta}) for obtaining optimal feedback controls.

\subsection{Positional Control Strategies}

We consider a formalization of feedback controls within the framework of positional control strategies \cite{Krasovskii_Subbotin_1988,Krasovskii_Krasovskii_1995} (see also \cite{Gomoyunov_2019_Trudy_Eng,Gomoyunov_2019_SIAM}).
By a (positional) control strategy, we mean an arbitrary function
\begin{equation*}
    G \times (0, \vartheta - t_0) \ni (t, w(\cdot), \eta)
    \mapsto U(t, w(\cdot), \eta) \in \mathbb{U},
\end{equation*}
where $\eta$ plays a role of some accuracy parameter.

Let us fix $\eta \in (0, \vartheta - t_0)$ and a partition $\Delta$ of the time interval $[t_0, \vartheta]$:
\begin{equation} \label{Delta}
    \Delta = \{\tau_j\}_{j = 1}^{k + 1},
    \quad \tau_1 = t_0, \quad \tau_{j+1} > \tau_j, \quad j \in \overline{1, k}, \quad \tau_{k+1} = \vartheta, \quad k \in \mathbb{N}.
\end{equation}
The triple $\{U, \eta, \Delta\}$ is called a control law.
This control law forms in system (\ref{system}) a piecewise constant control $u(\cdot) \in \mathcal{U}(t_0, \vartheta)$ by the following feedback rule:
\begin{equation} \label{Control_law}
    u(t) = U(\tau_j, x_{\tau_j}(\cdot), \eta),
    \quad t \in [\tau_j, \tau_{j + 1}), \quad j \in \overline{1, k},
\end{equation}
where, as usual, we denote $x_{\tau_j}(t) = x(t)$, $t \in [t_0, \tau_j]$.
Let us note that, for every initial value $x_0 \in B(R_x)$, the control law $\{U, \eta, \Delta\}$ determines the control $u(\cdot) = u(\cdot \mid t_0, x_0, \vartheta, U, \eta, \Delta)$ and the corresponding motion $x(\cdot)$ uniquely.
Let us note also that, in accordance with (\ref{Val}), we have
\begin{equation*}
    J\big( t_0, x_0, u(\cdot \mid t_0, x_0, \vartheta, U, \eta, \Delta) \big)
    \geq \rho(t_0, x_0).
\end{equation*}

Taking this into account, we call a control strategy $U^\circ$ optimal (uniformly with respect to initial values $x_0 \in B(R_x)$), if the following statement holds.
For any $\varepsilon > 0$, there exist
\begin{equation*}
    \eta^\circ = \eta^\circ(\varepsilon) \in (0, \vartheta - t_0),
    \quad (0, \eta^\circ] \ni \eta \mapsto \delta^\circ(\eta) = \delta^\circ(\varepsilon, \eta) \in (0, \infty)
\end{equation*}
such that, for any $\eta \in (0, \eta^\circ]$, any partition $\Delta$ (\ref{Delta}) with the diameter $\operatorname{diam}(\Delta) = \max_{j \in \overline{1, k}} (\tau_{j+1} - \tau_j) \leq \delta^\circ(\eta)$, and any $x_0 \in B(R_x)$, the inequality
\begin{equation} \label{optimal_control_strategy}
    J\big( t_0, x_0, u(\cdot \mid t_0, x_0, \vartheta, U^\circ, \eta, \Delta) \big)
    \leq \rho(t_0, x_0) + \varepsilon
\end{equation}
is valid, i.e., the control $u(\cdot \mid t_0, x_0, \vartheta, U^\circ, \eta, \Delta)$ is $\varepsilon$-optimal.

Following the ideas from Sect.~\ref{subsection_OCP_Auxiliary_Optimal_Control_Problem_with_Parameter}, let us construct such an optimal control strategy $U^\circ$ on the basis of optimal control strategies $P^\circ_\eta$ in the auxiliary problems (\ref{system_z_eta})--(\ref{cost_functional_z_eta}) for $\eta \in (0, \vartheta - t_0)$.

\subsection{Positional Control Strategies in the Auxiliary Problem}

Let us fix $\eta \in (0, \vartheta - t_0)$, consider the auxiliary optimal control problem (\ref{system_z_eta})--(\ref{cost_functional_z_eta}), and define the set $G_\eta = [t_0, \vartheta_\eta] \times \mathbb{R}^n$ of all positions of system (\ref{system_z_eta}).
A (positional) control strategy is a function
\begin{equation*}
    G_\eta \times (0, \vartheta - t_0) \ni (t, z, \varkappa) \mapsto P_\eta(t, z, \varkappa) \in \mathbb{U}.
\end{equation*}
where $\varkappa$ is treated as an accuracy parameter.

Let $\varkappa \in (0, \vartheta - t_0)$, and let $\Delta_\eta$ be a partition of the time interval $[t_0, \vartheta_\eta]$:
\begin{equation} \label{Delta_eta}
    \Delta_\eta = \{\tau_j\}_{j = 1}^{k + 1},
    \quad \tau_1 = t_0, \quad \tau_{j+1} > \tau_j, \quad j \in \overline{1, k}, \quad \tau_{k + 1} = \vartheta_\eta, \quad k \in \mathbb{N}.
\end{equation}
The control law $\{P_\eta, \varkappa, \Delta_\eta\}$ forms in the auxiliary system a piecewise constant control $p_\eta(\cdot) \in \mathcal{U}(t_0, \vartheta_\eta)$ by the following feedback rule:
\begin{equation} \label{Control_law_eta}
    p_\eta(t)
    = P_\eta(\tau_j, z(\tau_j), \varkappa),
    \quad t \in [\tau_j, \tau_{j + 1}), \quad j \in \overline{1, k}.
\end{equation}
For every $z_0 \in \mathbb{R}^n$, the control law $\{P_\eta, \varkappa, \Delta_\eta\}$ determines the control $p_\eta(\cdot) = p_\eta(\cdot \mid t_0, z_0, \vartheta_\eta, P_\eta, \varkappa, \Delta_\eta)$ and the corresponding motion $z(\cdot)$ uniquely.

In accordance with condition $(A.1)$, taking $M_F$ from (\ref{M_F_M_f}), let us denote
\begin{equation} \label{R_z}
    M_A
    = \esssup{t \in [t_0, \vartheta]} \|A(t)\|,
    \quad R_z
    = \big(1 + M_F M_A (\vartheta - t_0)^\alpha / \alpha \big) R_x.
\end{equation}
Then, due to (\ref{informational_image_Cauchy}), for every initial value $x_0 \in B(R_x)$ from (\ref{initial_condition}), we have
\begin{equation} \label{I_x_0_R_z}
    \|\mathcal{I}(t_0, x_0)\|
    \leq \Big\| \Big( \operatorname{Id}_n + \int_{t_0}^{\vartheta} \frac{F(\vartheta, \tau) A(\tau)}{(\vartheta - \tau)^{1 - \alpha}} \, \mathrm{d} \tau \Big) x_0 \Big\|
    \leq R_z,
\end{equation}
Therefore, in view of (\ref{initial_condition_z_eta}), in the auxiliary problem, we can restrict ourselves to initial values $z_0 \in B(R_z)$.

Thus, we call a control strategy $P_\eta^\circ$ optimal (uniformly with respect to $z_0 \in B(R_z)$), if the following statement holds.
For any $\varepsilon > 0$, there exist
\begin{equation} \label{kappa_delta_eta_a}
    \varkappa_\eta^\circ = \varkappa_\eta^\circ(\varepsilon) \in (0, \vartheta - t_0),
    \quad (0, \varkappa_\eta^\circ] \ni \varkappa \mapsto \omega_\eta^\circ(\varkappa) = \omega_\eta^\circ(\varepsilon, \varkappa) \in (0, \infty)
\end{equation}
such that, for any $\varkappa \in (0, \varkappa_\eta^\circ]$, any partition $\Delta_\eta$ (\ref{Delta_eta}) with the diameter $\operatorname{diam}(\Delta_\eta) \leq \omega_\eta^\circ(\varkappa)$, and any $z_0 \in B(R_z)$, the inequality below is valid:
\begin{equation} \label{optimal_control_strategy_a}
    J_\eta\big( t_0, z_0, p_\eta(\cdot \mid t_0, z_0, \vartheta_\eta, P_\eta^\circ, \varkappa, \Delta_\eta) \big)
    \leq \rho_\eta(t_0, z_0) + \varepsilon,
\end{equation}
i.e., $p_\eta(\cdot \mid t_0, z_0, \vartheta_\eta, P_\eta^\circ, \varkappa, \Delta_\eta)$ is an $\varepsilon$-optimal control in the auxiliary problem.

Let us note that, due to continuity of $f_\ast$ (see (\ref{f_ast_continuity})) and conditions $(A.3)$ and $(A.4)$, according to, e.g., \cite[Theorems~9.2 and~22.1]{Krasovskii_Krasovskii_1995} (see also \cite[Sect.~6]{Lukoyanov_Gomoyunov_2019_DGAA} and the references therein), such an optimal control strategy $P^\circ_\eta$ exists.
Moreover, it can be constructed, for example, by the method of extremal shift to accompanying points, which is shortly described below.

In the auxiliary problem, we consider the value function (see (\ref{cost_functional_z_eta}) and (\ref{Valaeta}))
\begin{equation*}
    \rho_\eta(t, z)
    = \inf_{p_\eta(\cdot) \in \mathcal{U}(t, \vartheta_\eta)}
    \Big( \sigma(z(\vartheta_\eta)) + \int_{t}^{\vartheta_\eta} \chi(\tau, p_\eta(\tau))  \, \mathrm{d} \tau \Big),
    \quad (t, z) \in G_\eta,
\end{equation*}
where $z(\cdot) = z(\cdot \mid t, z, \vartheta_\eta, p_\eta(\cdot))$.
Let $\varkappa \in (0, \vartheta - t_0)$ be fixed.
For every position $(t, z) \in G_\eta$, relying on the value function, we choose the accompanying point
\begin{equation} \label{OFC_extremal_shift_1}
    (z^\circ, z_{n + 1}^\circ)
    \in \argmin{(\bar{z}, \bar{z}_{n + 1})} (\rho_\eta(t, \bar{z}) + \bar{z}_{n + 1}),
\end{equation}
where the minimum is calculated over the pairs $(\bar{z}, \bar{z}_{n + 1}) \in \mathbb{R}^n \times \mathbb{R}$ such that
\begin{equation} \label{OFC_extremal_shift_2}
    \|z - \bar{z}\|^2 + \bar{z}_{n + 1}^2
    \leq r^2(t, \varkappa),
    \quad r^2(t, \varkappa)
    = \varkappa + (t - t_0) \varkappa,
\end{equation}
and, after that, we determine
\begin{equation} \label{OFC_extremal_shift_3}
    P^\circ_\eta(t, z, \varkappa)
    \in \argmin{p \in \mathbb{U}} \big( \langle z - z^\circ, f_\ast(t, p) \rangle - z_{n + 1}^\circ \chi(t, p) \big).
\end{equation}

Now, we define a control strategy $U^\ast$ in the original problem as follows.
For every $\eta \in (0, \vartheta - t_0)$, let us consider an optimal control strategy $P^\circ_\eta$ in the auxiliary problem and determine the corresponding number $\varkappa^\circ_\eta(\eta)$ according to (\ref{kappa_delta_eta_a}).
Then, for every $(t, w(\cdot)) \in G$, if $t < \vartheta_\eta$, we put
\begin{equation} \label{strategy_U_ast}
    U^\ast(t, w(\cdot), \eta)
    = P^\circ_\eta \big(t, \mathcal{I}(t, w(\cdot)), \varkappa^\circ_\eta (\eta) \big),
\end{equation}
where $\mathcal{I}(t, w(\cdot))$ is the informational image (see (\ref{informational_image}) or (\ref{informational_image_Cauchy})).
If $t \in [\vartheta_\eta, \vartheta]$, we formally define $U^\ast(t, w(\cdot), \eta) = \bar{u}$ for some fixed $\bar{u} \in \mathbb{U}$.

\begin{theorem} \label{Theorem_OFC}
    The control strategy $U^\ast$ defined by $(\ref{strategy_U_ast})$ is optimal in the problem $(\ref{system})$--$(\ref{cost_functional})$.
\end{theorem}

The proof of the theorem is given in Appendix.

Thus, in order to construct an optimal control strategy $U^\circ$ in the original problem in fractional-order system (\ref{system}), it is sufficient to find for every $\eta \in (0, \vartheta - t_0)$ an optimal control strategy $P^\circ_\eta$ in the auxiliary problem in first-order system (\ref{system_z_eta}).

\begin{remark} \label{Remark_Reduction_FC}
    As in Remark~\ref{Remark_Reduction}, let us suppose that the function $\sigma$ from cost functional (\ref{cost_functional}) can be represented as in (\ref{Remark_Reduction_sigma}).
    Then, for every $\eta \in (0, \vartheta - t_0)$, the auxiliary optimal control problem is described by the system
    \begin{equation*}
        \dot{z}(t) = K f_\ast(t, p_\eta(t)),
        \quad z(t) \in \mathbb{R}^{n_\ast}, \quad p_\eta(t) \in \mathbb{U}, \quad t \in [t_0, \vartheta_\eta],
    \end{equation*}
    and the cost functional
    \begin{equation*}
        J_\eta(t_0, z_0, p_\eta(\cdot))
        = \mu(z(\vartheta))
        + \int_{t_0}^{\vartheta_\eta} \chi(t, p_\eta(t)) \, \mathrm{d} t,
        \quad p_\eta(\cdot) \in \mathcal{U}(t_0, \vartheta_\eta).
    \end{equation*}
    Let $P^\circ_\eta$ be an optimal control strategy in this auxiliary problem.
    In accordance with (\ref{strategy_U_ast}), we put
    \begin{equation} \label{strategy_U_ast_Remark_Reduction}
        U^\ast (t, w(\cdot), \eta)
        = P^\circ_\eta \big(t, K(\mathcal{I}(t, w(\cdot)) - c), \varkappa^\circ_\eta (\eta) \big)
    \end{equation}
    for $(t, w(\cdot)) \in G$ such that $t < \vartheta_\eta$.
    By analogy with Theorem~\ref{Theorem_OFC}, one can prove that such a control strategy $U^\ast$ is optimal in the original problem.
\end{remark}

The following two examples illustrate Theorem~\ref{Theorem_OFC} and Remark~\ref{Remark_Reduction_FC}.

\begin{example}
    Let us consider the optimal control problem for system (\ref{Example_2_system}) and the cost functional
    \begin{equation*}
        J(t_0, x_0, u(\cdot))
        = (x_1(\vartheta) - c_1)^2,
        \quad u(\cdot) \in \mathcal{U}(t_0, \vartheta),
    \end{equation*}
    where $c_1 \in \mathbb{R}$ is a given number.
    Then, for every $\eta \in (0, \vartheta - t_0)$, we come to the auxiliary optimal control problem for system (\ref{Example_2_system_auxiliary}) and the cost functional
    \begin{equation*}
        J(t_0, z_0, p_\eta(\cdot))
        = z^2(\vartheta),
        \quad p_\eta(\cdot) \in \mathcal{U}(t_0, \vartheta_\eta).
    \end{equation*}
    Let us associate this auxiliary problem with the Cauchy problem for the corresponding Hamilton--Jacobi--Bellman equation
    \begin{equation*}
        \frac{\partial \varphi(t, z)}{\partial t} - b_2(t) \Big| \frac{\partial \varphi(t, z)}{\partial z} \Big|
        = 0,
        \quad (t, z) \in (t_0, \vartheta_\eta) \times \mathbb{R},
    \end{equation*}
    with the right-end boundary condition $\varphi(\vartheta_\eta, z) = z^2$, $z \in \mathbb{R}$.
    Let us define
    \begin{equation*}
        \varphi_\eta(t, z)
        = \begin{cases}
            (|z| - \psi_\eta(t))^2, & \mbox{if } |z| > \psi_\eta(t), \\
            0, & \mbox{if } |z| \leq \psi_\eta(t),
          \end{cases}
          \quad (t, z) \in G_\eta = [t_0, \vartheta_\eta] \times \mathbb{R},
    \end{equation*}
    where, for the function $b_2$ from (\ref{Example_1_b}), we denote
    \begin{equation*}
        \psi_\eta(t)
        = \int_{t_0}^{\vartheta_\eta} b_2(t) \, \mathrm{d} t
        = \frac{(\vartheta - t)^{2 \alpha} - \eta^{2 \alpha}}{\Gamma(2 \alpha + 1)},
        \quad t \in [t_0, \vartheta_\eta].
    \end{equation*}
    One can verify that $\varphi_\eta$ is a continuously differentiable solution to the considered Cauchy problem.
    Thus, according to, e.g., \cite[Theorem~4.1.1]{Krasovskii_Subbotin_1988}, we have $\rho_\eta(t, z) = \varphi_\eta(t, z)$, $(t, z) \in G_\eta$, and the optimal control strategy in the auxiliary problem is given by
    \begin{equation*}
        P^\circ_\eta(t, z)
        \in \argmin{p \in [-1, 1]} \Big(p \frac{\partial \varphi_\eta(t, z)}{\partial z}\Big)
        = \begin{cases}
            \{1\}, & \mbox{if } z < - \psi_\eta(t), \\
            [-1, 1], & \mbox{if } |z| \leq \psi_\eta(t), \\
            \{- 1\}, & \mbox{if } z > \psi_\eta(t),
          \end{cases}
    \end{equation*}
    where $(t, z) \in [t_0, \vartheta_\eta) \times \mathbb{R}$.

    Hence, by Theorem~\ref{Theorem_OFC} and Remark~\ref{Remark_Reduction_FC}, on the basis of the found optimal control strategies $P^\circ_\eta$, $\eta \in (0, \vartheta - t_0)$, we can construct an optimal control strategy in the original problem as follows:
    \begin{equation*}
        U^\circ (t, w(\cdot), \eta)
        = \begin{cases}
            P^\circ_\eta \big(t, \mathcal{I}_1(t, w(\cdot)) - c_1\big), & \mbox{if } t < \vartheta_\eta, \\
            \bar{u}, & \mbox{if } t \geq \vartheta_\eta,
        \end{cases}
    \end{equation*}
    where $(t, w(\cdot)) \in G$ and $\bar{u} \in \mathbb{U}$, and $\mathcal{I}_1(t, w(\cdot))$ is the first coordinate of the informational image $\mathcal{I}(t, w(\cdot))$ defined by (\ref{informational_image}) or (\ref{informational_image_Cauchy}).
\end{example}

\begin{example}
    Let us consider the original optimal control problem (\ref{system})--(\ref{cost_functional}) in the case when $\sigma(x) = \|K (x - c)\|$, $x \in \mathbb{R}^n$ (see Remark~\ref{Remark_Reduction_FC}).
    Then, in the auxiliary optimal control problem (\ref{system_z_eta})--(\ref{cost_functional_z_eta}), the cost functional takes the form
    \begin{equation*}
        J_\eta(t_0, z_0, p_\eta(\cdot))
        = \|z(\vartheta_\eta)\|
        + \int_{t_0}^{\vartheta_\eta} \chi (t, p_\eta(t)) \, \mathrm{d} t,
        \quad p_\eta(\cdot) \in \mathcal{U}(t_0, \vartheta_\eta).
    \end{equation*}
    Hence, according to, e.g., \cite[\S~23]{Krasovskii_Krasovskii_1995} (see also \cite[Sect.~7.1]{Lukoyanov_Gomoyunov_2019_DGAA} and the references therein), the value function in the auxiliary problem is given by
    \begin{equation*}
        \rho_\eta (t, z)
        = \max_{l \in B(1)} \big( \langle l, z \rangle + \nu_\eta(t, l) \big),
        \quad (t, z) \in G_\eta = [t_0, \vartheta_\eta] \times \mathbb{R}^{n_\ast},
    \end{equation*}
    where $B(1) = \{l \in \mathbb{R}^{n_\ast}: \|l\| \leq 1\}$ and
    \begin{equation*}
        \nu_\eta(t, l)
        = \int_{t}^{\vartheta_\eta} \min_{p \in \mathbb{U}}
        \big( \langle l, K f_\ast(\tau, p) \rangle + \chi(\tau, p) \big) \, \mathrm{d} \tau,
        \quad t \in [t_0, \vartheta_\eta], \quad l \in B(1).
    \end{equation*}
    Moreover, applying the method of extremal shift to accompanying points (\ref{OFC_extremal_shift_1})--(\ref{OFC_extremal_shift_3}), one can construct an optimal control strategy as follows:
    \begin{equation*}
        P^\circ_\eta(t, z, \varkappa)
        \in \argmin{p \in \mathbb{U}} \big( \langle l^\circ, K f_\ast(t, p) \rangle + \chi(t, p) \big),
        \quad (t, z) \in G_\eta, \quad \varkappa \in (0, \vartheta - t_0),
    \end{equation*}
    where
    \begin{equation*}
        l^\circ
        \in \argmax{l \in B(1)} \big( \langle l, z \rangle + \nu_\eta(t, l) - r(t, \varkappa) \sqrt{1 + \|l\|^2} \big).
    \end{equation*}
    Thus, based on these formulas and (\ref{strategy_U_ast_Remark_Reduction}), we can effectively calculate an optimal control strategy in the original problem.
\end{example}

\begin{acknowledgements}
    This work was supported by RSF, project no. 19-11-00105.
\end{acknowledgements}


\section{Introduction}

The celebrated \textsc{Feller} equation was introduced in two seminal papers published by William \textsc{Feller} (1951/1952) in Annals of Mathematics \cite{Feller1951, Feller1952}. These publications studied mathematically (and, henceforth, gave the name) to the following equation\footnote{To be more accurate, W.~\textsc{Feller} studied the following equation \cite{Feller1951}: \begin{equation*}
  p_{\,t}\ =\ \bigl[\,a\,x\,u\,\bigr]_{\,x\,x}\ -\ \bigl[\,(c\ +\ b\,x)\,u\,\bigr]_{\,x}\,,
\end{equation*} where $a\ >\ 0$ and $0\ <\ x\ <\ +\,\infty\,$.}:
\begin{equation}\label{eq:feller}
  p_{\,t}\ +\ \mathfrak{F}_{\,x}\ =\ 0\,, \qquad \mathfrak{F}\,(p,\,x,\,t)\ \mathop{\stackrel{\,\mathrm{def}}{:=}\,}\ -x\cdot\bigl(\gamma\,\,p\ +\ \eta\,p_{\,x}\bigr)\,,
\end{equation}
where the subscripts $t\,$, $x$ denote the partial derivatives, \emph{i.e.}\xspace $(\cdot)_{\,t}\ \mathop{\stackrel{\,\mathrm{def}}{:=}\,}\ \PD{(\cdot)}{t}\,$, $(\cdot)_{\,x}\ \mathop{\stackrel{\,\mathrm{def}}{:=}\,}\ \PD{(\cdot)}{x}\,$. Two parameters $\gamma$ and $\eta\ >\ 0$ can be time dependent in some physical applications, even if in this study we assume they are constants, for the sake of simplicity\footnote{The numerical method we are going to propose can be straightforwardly generalized for this case when $\gamma\ =\ \gamma\,(t)$ and $\eta\ =\ \eta\,(t)\,$. Moreover, \textsc{Feller}'s processes with time-varying coefficients were studied recently in \cite{Masoliver2016}.}. Equation~\eqref{eq:feller} can be seen as the \textsc{Fokker}--\textsc{Planck} (or the forward \textsc{Kolmogorov}) equation with $\gamma\,x$ being the drift and $\eta\,x$ being the diffusion coefficients. One can notice also that Equation~\eqref{eq:feller} becomes singular at $x\ =\ 0$ and $x\ =\ +\,\infty\,$. We remind that practically important solutions to \textsc{Feller}'s equation might be unbounded near $x\ =\ 0\,$. In order to attempt at solving Equation~\eqref{eq:feller}, one has to prescribe an initial condition $p\,(x,\,0)\ =\ p_{\,0}\,(x)$ presumably with a boundary condition at $x\ =\ 0\,$. A popular choice is to prescribe the homogeneous boundary condition $p\,(0,\,t)\ \equiv\ 0\,$. For this choice of the boundary condition it is not difficult to show that the \textsc{Feller} equation dynamics would preserve solution positivity provided that $p_{\,0}\,(x)\ \geq\ 0$ (see Appendix~\ref{app:pos} for a proof). The solution norm is also preserved (see Appendix~\ref{app:norm}). Moreover, \textsc{Feller} using the \textsc{Laplace} transform techniques has shown in \cite{Feller1951} that the initial condition $p_{\,0}\,(x)$ determines uniquely the solution. In other words, \emph{no} boundary condition at $x\ =\ 0$ should be prescribed. This conclusion might appear, perhaps, to be counter-intuitive.

The great interest in \textsc{Feller}'s equation can be explained by its connection to \textsc{Feller}'s processes, which can be described by the following stochastic differential \textsc{Langevin} equation\footnote{The stochastic differential equations is understood in the sense of \textsc{It\={o}}.}:
\begin{equation*}
  \mathrm{d} X_{\,t}\ =\ -\,\gamma\,X_{\,t}\,\mathrm{d} t\ +\ \sqrt{2\,\eta\,X_{\,t}}\,\mathrm{d}\mathcal{W}_{\,t}\,,
\end{equation*}
where $\mathcal{W}_{\,t}$ is the standard \textsc{Wiener} process, \emph{i.e.}\xspace $\xi\,(t)\ \mathop{\stackrel{\,\mathrm{def}}{:=}\,}\ \odd{\mathcal{W}_{\,t}}{t}$ is zero-mean \textsc{Gaussian} white noise, \emph{i.e.}\xspace
\begin{equation*}
  \langle \xi\,(t)\rangle\ =\ 0\,, \qquad
  \langle \xi\,(t)\,\xi\,(s)\rangle\ =\ \delta\,(t\ -\ s)\,,
\end{equation*}
where the brackets $\langle \cdot\rangle$ denote an ensemble averaging operator. Then, the Probability Density Function (PDF) $p\,(x,\,t;\,x_{\,0})$ of the process $X\,(t)\,$, \emph{i.e.}\xspace
\begin{equation*}
  \mathbb{P}\,\bigl\{x\ <\ X\,(t)\ <\ x\ +\ \mathrm{d} x\ \vert\ X\,(0)\ =\ x_{\,0}\bigr\}\ \equiv\ p\,(x,\,t;\,x_{\,0})\,\mathrm{d} x
\end{equation*}
satisfies Equation~\eqref{eq:feller} with the following initial condition \cite{Gardiner2004}:
\begin{equation*}
  p_{\,0}\,(x)\ =\ \delta\,(x\ -\ x_{\,0})\,, \qquad x_{\,0}\ \in\ \mathds{R}^{\,+}\,.
\end{equation*}
The point $x\ =\ 0$ is a singular boundary that the process $X\,(t)$ cannot cross. The \textsc{Feller} process is a continuous representation of branching and birth-death processes, which never attains negative values. This property makes it an ideal model not only in physical, but also in biological and social sciences \cite{Murray2007, Gan2015, Masoliver2016}.

As a general comprehensive reference on generalized \textsc{Feller}'s equations we can mention the book \cite{Lehnikg1993}. Since at least a couple of years there is again a growing interest for studying equation \eqref{eq:feller}. Some singular solutions to \textsc{Feller}'s equation with constant coefficients were constructed in \cite{Gan2015} via spectral decompositions. \textsc{Feller}'s equation and \textsc{Feller}'s processes with time-varying coefficients were studied analytically (always using the \textsc{Laplace} transform) and asymptotically in \cite{Masoliver2016}.

Recently, the \textsc{Feller} equation was derived in the context of the weakly interacting random waves dominated by four-wave interactions \cite{Choi2005}. Wave Turbulence\footnote{We could define the Wave Turbulence (WT) as a physical and mathematical study of systems where random and coherent waves co-exist and interact \cite{Zakharov1992}.} (WT) is a common name for such processes \cite{Zakharov1992}. In WT the \textsc{Feller} equation governs the PDF of squared \textsc{Fourier} wave amplitudes, \emph{i.e.}\xspace $x\ \sim\ \abs{a}^{\,2}\,$. In \cite{Choi2005} some steady solutions to this equation with finite flux in the amplitude space were constructed\footnote{There is a misprint in \cite[p.~366]{Choi2005}. To obtain mathematically correct solutions one has to define $n_{\,k}\ \mathop{\stackrel{\,\mathrm{def}}{:=}\,}\ \dfrac{\eta}{\gamma}$ on the line below Equation (14).}. See also \cite[Chapter~11]{Nazarenko2011} for a detailed discussion and interpretations.

The behaviour of solutions $p\,(x,\,t)$ for large $x$ describes the appearance probability of extreme waves. In the context of ocean waves, these extreme events are known as \emph{rogue} (or \emph{freak}) waves \cite{Dysthe2008}. In the WT literature, any noticeable deviation from the \textsc{Rayleigh} distribution for $x\ \gg\ 1$ is referred to as the \emph{anomalous probability distribution} of large amplitude waves \cite{Choi2005}. For \textsc{Gaussian} wave fields all statistical properties can be derived from the spectrum. However, the PDFs and other higher order moments are compulsory tools to study such deviations.

The present study focuses on the numerical discretization and simulation of \textsc{Feller} equation. The naive approach to solve this equation numerically encounters notorious difficulties. The first question, which arises is what is the (numerical) boundary condition to be imposed at $x\ =\ 0\,$? Moreover, one can notice that Equation \eqref{eq:feller} is posed on a semi-infinite domain. There are three main strategies to tackle this difficulty:
\begin{enumerate}
  \item Map $\mathds{R}^{\,+}$ on a finite interval $[\,0,\,\ell\,]$
  \item Use spectral expansions on $\mathds{R}^{\,+}$ (\emph{e.g.}\xspace \textsc{Laguerre} or associated \textsc{Laguerre} polynomials)
  \item Replace (truncate) $\mathds{R}^{\,+}$ to $[\,0,\,L\,]\,$, with $L\ \gg\ 1\,$.
\end{enumerate}
In most studies the latter option is retained by imposing some appropriate boundary conditions at the artificial boundary $x\ =\ L\,$. In our study we shall propose a method, which is able to handle the semi-infinite domain $\mathds{R}^{\,+}$ without any truncations or simplifications. Finally, the diffusion coefficient in the \textsc{Feller} equation \eqref{eq:feller} is unbounded. If the domain is truncated at $x\ =\ L\,$, then the diffusion coefficient takes the maximal value $\nu_{\,\max}\ \mathop{\stackrel{\,\mathrm{def}}{:=}\,}\ \eta\,L\ \gg\ 1\,$, which depends on the truncation limit $L$ and it can become very large in practice. We remind also that explicit schemes for diffusion equations are subject to the so-called \textsc{Courant}--\textsc{Friedrichs}--\textsc{Lewy} (CFL) stability conditions \cite{Courant1928}:
\begin{equation*}
  \Delta t\ \leq\ \frac{\Delta x^{\,2}}{2\,\nu_{\,\max}}\,.
\end{equation*}
Taking into account the fact that $\nu_{\,\max}$ can be arbitrarily large, no explicit scheme can be usable with \textsc{Feller} equation in practice. Moreover, the dynamics of the \textsc{Feller} equation spreads over the space $\mathds{R}^{\,+}$ even localized initial conditions. In general, one can show that the support of $p\,(x,\,t)\,$, $t\ >\ 0$ is strictly larger\footnote{Using modern analytical techniques it is possible to show even sharper results on the solution support, see \emph{e.g.}\xspace \cite{Carrillo2004}.} than the one of $p\,(x,\,0)\,$. It is the so-called \emph{retention property}. Thus, longer simulation times require larger domains. For all these reasons, it becomes clear that numerical discretization of the \textsc{Feller} equation requires special care.

In this study we demonstrate how to overcome this assertion as well. The main idea behind our study is to bring together PDEs and Fluid Mechanics. First, we observe that the classical \textsc{Eulerian} description is not suitable for this equation, even if the problem is initially formulated in the \textsc{Eulerian} setting. Consequently, the \textsc{Feller} equation will be recast in special \emph{material} or the so-called \textsc{Lagrangian} variables\footnote{It is known that both \textsc{Eulerian} and \textsc{Lagrangian} descriptions were proposed by the same person --- Leonhard \textsc{Euler}.}, which make the resolution easier and naturally adaptive \cite[Chapter~7]{Gosse2013}.

The present manuscript is organized as follows. The symmetry analysis of Equation~\eqref{eq:feller} is performed in Section~\ref{sec:sym}. Then, the governing equation is reformulated in \textsc{Lagrangian} variables in Section~\ref{sec:reform}. The numerical results are presented in Section~\ref{sec:num}. Finally, the main conclusions and perspectives are outlined in Section~\ref{sec:disc}.

\section{Symmetry analysis}
\label{sec:sym}

In general, a linear PDE admits an infinity of conservation laws, with integrating multipliers being solutions to the adjoint PDE \cite{Bluman2010}. Here we provide an interesting conservation law, which was found using the \textsc{GeM} \textsc{Maple} package \cite{Cheviakov2007}:
\begin{equation*}
  \Bigl(\E_{\,1}\,\bigl(-\frac{\gamma\,x}{\eta}\bigr)\,p\Bigr)_{\,t}\ +\ \mathfrak{G}_{\,x}\ =\ 0\,,
\end{equation*}
where $\E_{\,1}\,(z)\ \mathop{\stackrel{\,\mathrm{def}}{:=}\,}\ \int_{\,1}^{\,+\infty}\,\frac{\mathrm{e}^{-\,t\,z}}{t}\;\mathrm{d} t$ is the so-called \emph{exponential integral} function \cite{Abramowitz1965} and the flux $\mathfrak{G}$ is defined as
\begin{equation*}
  \mathfrak{G}\,(x,\,p)\ \mathop{\stackrel{\,\mathrm{def}}{:=}\,}\ -\eta\,\mathrm{e}^{\,\frac{\gamma\,x}{\eta}}\,p\ -\ x\,\Bigl(\gamma\,\E_{\,1}\,\bigl(-\frac{\gamma\,x}{\eta}\bigr)\,p\ +\ \eta\,\E_{\,1}\,\bigl(-\frac{\gamma\,x}{\eta}\bigr)\,p_{\,x}\Bigr)\,.
\end{equation*}

The symmetry group of point transformations can be computed using \textsc{GeM} package as well. The infinitesimal generators are given below:
\begin{align*}
  \text{\textxi}_{\,1}\ &=\ \mathscr{D}_{\,t}\,, \\
  \text{\textxi}_{\,2}\ &=\ p\,\mathscr{D}_{\,p}\,, \\
  \text{\textxi}_{\,3}\ &=\ -\,\frac{\mathrm{e}^{\,-\,\gamma\,t}}{\gamma}\;\mathscr{D}_{\,t}\ +\ \mathrm{e}^{\,-\,\gamma\,t}\,x\,\mathscr{D}_{\,x}\ -\ \mathrm{e}^{\,-\,\gamma\,t}\,p\,\mathscr{D}_{\,p}\,, \\
  \text{\textxi}_{\,4}\ &=\ \frac{\mathrm{e}^{\,\gamma\,t}}{\gamma}\;\mathscr{D}_{\,t}\ +\ \mathrm{e}^{\,\gamma\,t}\,x\,\mathscr{D}_{\,x}\ -\ \frac{\gamma\,\mathrm{e}^{\,\gamma\,t}}{\eta}\;x\,p\,\mathscr{D}_{\,p}\,, \\
  \text{\textxi}_{\,5}\ &=\ \mathrm{e}^{\,(\gamma\ +\ c)\,t}\,\M\Bigl(1\ +\ \frac{c}{\gamma},\,1,\,\frac{\gamma\,x}{\eta}\Bigr)\,\mathrm{e}^{\,-\,\frac{\gamma\,x}{\eta}}\;\mathscr{D}_{\,p}\,, \\
  \text{\textxi}_{\,6}\ &=\ \mathrm{e}^{\,(\gamma\ +\ c)\,t}\,\U\Bigl(1\ +\ \frac{c}{\gamma},\,1,\,\frac{\gamma\,x}{\eta}\Bigr)\,\mathrm{e}^{\,-\,\frac{\gamma\,x}{\eta}}\;\mathscr{D}_{\,p}\,,
\end{align*}
where $c\ \in\ \mathds{R}\,$, $\M(a,\,b,\,z)$ and $\U(a,\,b,\,z)$ are \textsc{Kummer} special functions \cite{Kummer1837, Abramowitz1965} (see also Appendix~\ref{app:kum}). The corresponding point transformations, which map solutions of \eqref{eq:feller} into other solutions, can be readily obtained by integrating several ODE systems (we do not provide integration details here):
\begin{equation*}
  \bigl(\,t,\,x,\,p\,\bigr)\ \mapsto\ \bigl(\,t\ +\ \upvarepsilon_{\,1},\,x,\,p\,\bigr)\,,
\end{equation*}
\begin{equation*}
  \bigl(\,t,\,x,\,p\,\bigr)\ \mapsto\ \bigl(\,t,\,x,\,\mathrm{e}^{\,\upvarepsilon_{\,2}}\,p\,\bigr)\,,
\end{equation*}
\begin{equation*}
  \bigl(\,t,\,x,\,p\,\bigr)\ \mapsto\ \Bigl(\frac{1}{\gamma}\;\ln\bigl(\upvarepsilon_{\,3}\,\gamma\ +\ \mathrm{e}^{\,\gamma\,t}\bigr),\,\frac{\mathrm{e}^{\,\gamma\,t}}{\upvarepsilon_{\,3}\,\gamma\ +\ \mathrm{e}^{\,\gamma\,t}}\;x,\,\bigl(1\ +\ \upvarepsilon_{\,3}\,\gamma\,\mathrm{e}^{\,-\,\gamma\,t}\bigr)\,p\Bigr)\,,
\end{equation*}
\begin{equation*}
  \bigl(\,t,\,x,\,p\,\bigr)\ \mapsto\ \Bigl(t\ -\ \frac{1}{\gamma}\;\ln\bigl(1\ -\ \upvarepsilon_{\,4}\,\gamma\,\mathrm{e}^{\gamma\,t}\bigr)\,, \frac{x}{1\ -\ \upvarepsilon_{\,4}\,\gamma\,\mathrm{e}^{\,\gamma\,t}},\,\mathrm{e}^{-\,\dfrac{\upvarepsilon_{\,4}\,\gamma^{\,2}\,x\,\mathrm{e}^{\,\gamma\,t}}{\eta\,(1\ -\ \upvarepsilon_{\,4}\,\gamma\,\mathrm{e}^{\,\gamma\,t})}}\cdot p\Bigr)\,,
\end{equation*}
\begin{equation*}
  \bigl(\,t,\,x,\,p\,\bigr)\ \mapsto\ \Bigl(t,\,x,\,p\ +\ \upvarepsilon_{\,5}\,\M\bigl(1\ +\ \frac{c}{\gamma},\,1,\,\frac{\gamma\,x}{\eta}\bigr)\,\mathrm{e}^{-\,\frac{\gamma\,x}{\eta}\ +\ (\gamma\ +\ c)\,t}\Bigr)\,,
\end{equation*}
\begin{equation*}
  \bigl(\,t,\,x,\,p\,\bigr)\ \mapsto\ \Bigl(t,\,x,\,p\ +\ \upvarepsilon_{\,6}\,\U\bigl(1\ +\ \frac{c}{\gamma},\,1,\,\frac{\gamma\,x}{\eta}\bigr)\,\mathrm{e}^{-\,\frac{\gamma\,x}{\eta}\ +\ (\gamma\ +\ c)\,t}\Bigr)\,.
\end{equation*}
The first symmetry is the time translation. The second one is the scaling of the dependent variable (the governing equation is linear). Symmetries $3$ and $4$ are exponential scalings. Two last symmetries express the fact that we can always add to the solution a particular solution to the homogeneous equation to obtain another solution. For instance, the solutions invariant under time translations ($\text{\textxi}_{\,1}$) are steady states and their general form is the following:
\begin{equation}\label{eq:steady}
  p\,(x)\ =\ \mathrm{e}^{\,-\,\frac{\gamma\,x}{\eta}}\Bigl(\mathscr{C}_{\,1}\,\E_{\,1}\,\bigl(-\frac{\gamma\,x}{\eta}\bigr)\ +\ \mathscr{C}_{\,2}\Bigr)\,,
\end{equation}
where $\mathscr{C}_{\,1,\,2}$ are `arbitrary' constants, which have to be determined from imposed conditions. Of course, they should be chosen so that the resulting steady solution is a valid probability distribution. It is not difficult to check that the imposed flux $\mathfrak{F}$ on the steady state solution is equal to $\mathscr{C}_{\,1}\,\eta\,$. Please, notice also an important property of the exponential integral function, which is useful in manipulating its values for negative arguments:
\begin{equation*}
  \lim_{\delta\ \to\ 0} \E_{\,1}\,(x\ \pm\ \mathrm{i}\,\delta)\ =\ \E_{\,1}\,(x)\ \mp\ \mathrm{i}\,\pi\ \equiv\ -\,\Ei(-x)\ \mp\ \mathrm{i}\,\pi\,, \qquad
  x\ <\ 0\,,
\end{equation*}
where $\Ei(z)$ is the following exponential integral:
\begin{equation*}
  \Ei(z)\ \mathop{\stackrel{\,\mathrm{def}}{:=}\,}\ -\,\int_{\,-\infty}^{\,z}\,\frac{\mathrm{e}^{\,t}}{t}\;\mathrm{d} t\,.
\end{equation*}
Thus, $\Ei(x)\ \equiv\ -\E_{\,1}\,(-x)$ for $x\ <\ 0\,$. The last relation can be also extended to the entire complex plain:
\begin{equation*}
  \Ei(z)\ \equiv\ -\E_{\,1}\,(-z)\ +\ \frac{1}{2}\;\ln z\ -\ \frac{1}{2}\;\ln\Bigl(\,\frac{1}{z}\,\Bigr)\ -\ \ln(-z)\,, \qquad z\ \in\ \mathds{C}\,.
\end{equation*}
We provide here also the general solutions invariant under the symmetry $(\text{\textxi}_{\,3})\,$:
\begin{equation*}
  p\,(x,\,t)\ =\ \Bigl(\mathscr{C}_{\,2}\ -\ \mathscr{C}_{\,1}\,t\ +\ \frac{\mathscr{C}_{\,1}}{\gamma}\;\ln x\Bigr)\,\mathrm{e}^{-\,\frac{\gamma\,x}{\eta}}
\end{equation*}
and under symmetry $(\text{\textxi}_{\,4})\,$:
\begin{equation*}
  p\,(x,\,t)\ =\ \Bigl(\mathscr{C}_{\,1}\ +\ \mathscr{C}_{\,2}\,t\ +\ \frac{\mathscr{C}_{\,2}}{\gamma}\;\ln x\Bigr)\,\mathrm{e}^{\,\gamma\,t}\,.
\end{equation*}
These solutions might be used, for example, to validate numerical codes.

\section{Reformulation}
\label{sec:reform}

By following the lines of \cite[Chapter~7]{Gosse2013}, we are going to rewrite \textsc{Feller}'s Equation~\eqref{eq:feller} with the so-called \textsc{Lagrangian} or \emph{material} variables. The main advantage of this formulation consists in the fact that we can handle infinite domains \emph{without} any truncations, transformations, \emph{etc.}\xspace It becomes possible to carry computations in infinite domains. Our domain is semi-infinite ($x\ \in\ \mathds{R}^{\,+}$) with the left boundary $x\ =\ 0$ being a reflection point.

As the first step, we introduce the distribution function associated to the probability density $p\,(x,\,t)\,$:
\begin{equation}\label{eq:78}
  \P\,(x,\,t)\ \mathop{\stackrel{\,\mathrm{def}}{:=}\,}\ \int_{\,0}^{\,x}\,p\,(\xi,\,t)\,\mathrm{d}\xi\,.
\end{equation}
The same can be done for the initial condition as well:
\begin{equation*}
  \P_{\,0}(x)\ \mathop{\stackrel{\,\mathrm{def}}{:=}\,}\ \int_{\,0}^{\,x}\,p_{\,0}\,(\xi)\,\mathrm{d}\xi\,, \qquad
  p_{\,0}\ \in\ W_{\,loc}^{\,1,\,1}\,(\mathds{R}^{\,+})\,.
\end{equation*}
We notice also two obvious properties of the function $\P\,(x,\,t)\,$:
\begin{equation*}
  \partial_{\,x}\P\,(x,\,t)\ \equiv\ p\,(x,\,t)\,, \qquad
  \lim_{x\ \to\ 0} \P\,(x,\,t)\ =\ 0\,, \qquad
  \lim_{x\ \to\ +\infty} \P\,(x,\,t)\ =\ 1\,.
\end{equation*}
Due to the positivity preservation property (see Appendix~\ref{app:pos}) the function $\P\,(x,\,t)$ is non-decreasing in variable $x\,$. Thus, we can define its pseudo-inverse\footnote{This mapping is sometimes called the \emph{reciprocal mapping} \cite{Gosse2013} or an \emph{order preserving string} \cite{Brenier2004}.}:
\begin{equation*}
  \mathscr{X}\,:\ [\,0,\,1\,]\times\mathds{R}^{\,+}\ \mapsto\ \mathds{R}^{\,+}\,,
\end{equation*}
which is defined as
\begin{equation*}
  \mathscr{X}\,(\bar{\P},\,t)\ \mathop{\stackrel{\,\mathrm{def}}{:=}\,}\ \inf\;\{\,\xi\ \in\ \mathds{R}^{\,+}\ \vert\ \P\,(\xi,\,t)\ =\ \bar{\P}\,\}\,.
\end{equation*}
Similarly, the initial condition does possess a pseudo-inverse as well:
\begin{equation}\label{eq:initp}
  \mathscr{X}_{\,0}\,(\bar{\P})\ \mathop{\stackrel{\,\mathrm{def}}{:=}\,}\ \inf\;\{\,\xi\ \in\ \mathds{R}^{\,+}\ \vert\ \P_{\,0}\,(\xi)\ =\ \bar{\P}\,\}\,,
\end{equation}
such that $\mathscr{X}\,(\bar{\P},\,0)\ \equiv\ \mathscr{X}_{\,0}\,(\bar{\P})\,$.

If \textsc{Feller} equation \eqref{eq:feller} holds in the sense of distributions, then the following equation holds as well:
\begin{equation}\label{eq:shit}
  \P_{\,t}\ -\ x\,\Bigl[\,\gamma\,\P\ +\ \eta\,\P_{\,x}\,\Bigr]_{\,x}\ =\ 0\,.
\end{equation}
along with the initial condition
\begin{equation*}
  \P\,(x,\,0)\ =\ \P_{\,0}(x)\,.
\end{equation*}
Equation \eqref{eq:shit} can be readily obtained by exploiting the obvious property $p\,(x,\,t)\ =\ \partial_{\,x}\,\P\,(x,\,t)\,$. In Appendix~\ref{app:pos} we show that zero value of the solution $p\,(x,\,t)$ is repulsive. Thus, $\partial_{\,x}\,\P\,(x,\,t)\ \equiv\ p\,(x,\,t)\ >\ 0\,$, $\ \forall\, (x,\,t)\ \in\ \bigl(\mathds{R}^{\,+}\bigr)^{\,2}\,$. Thus, the implicit function theorem \cite{Zorich2004, zorich} guarantees the existence of derivatives of the inverse mapping $\mathscr{X}\,(\bar{\P},\,t)\,$. Let us compute them by differentiating with respect to $\bar{\P}$ and $t$ the following obvious identity:
\begin{equation*}
  \P\,(\mathscr{X}\,(\bar{\P},\,t),\,t)\ \equiv\ \bar{\P}\,.
\end{equation*}
Thus, one can easily show that
\begin{equation*}
  \pd{\mathscr{X}}{\P}\ =\ \frac{1}{\partial_{\,x}\,\P}\,, \qquad
  \pd{\mathscr{X}}{t}\ =\ -\,\frac{\partial_{\,t}\,\P}{\partial_{\,x}\,\P}\,.
\end{equation*}
Using these expressions of partial derivatives, we derive the following evolution equation for the inverse mapping $\mathscr{X}\,(\P,\,\cdot)\,$:
\begin{equation}\label{eq:lagr}
  \bigl(\mathrm{e}^{\,\gamma\,t}\,\mathscr{X}\bigr)_{\,t}\ +\ \mathscr{X}\cdot\Bigl[\,\eta\,\mathrm{e}^{\,\gamma\,t}\,\Bigl(\pd{\mathscr{X}}{\P}\Bigr)^{\,-1}\,\Bigr]_{\,\P}\ =\ 0\,.
\end{equation}
The last equation can be rewritten also by introducing a new dynamic variable:
\begin{equation}\label{eq:xy}
  \mathscr{Y}\,(\P,\,t)\ \mathop{\stackrel{\,\mathrm{def}}{:=}\,}\ \mathrm{e}^{\,\gamma\,t}\,\mathscr{X}\,(\P,\,t)\,, \qquad
  \mathscr{Y}\,(\P,\,0)\ \equiv\ \mathscr{X}\,(\P,\,0)\,.
\end{equation}
It is not difficult to see that Equation~\eqref{eq:lagr} becomes:
\begin{equation}\label{eq:y}
  \mathscr{Y}_{\,t}\ +\ \mathscr{Y}\cdot\Bigl[\,\eta\,\mathrm{e}^{\,\gamma\,t}\,\Bigl(\pd{\mathscr{Y}}{\P}\Bigr)^{\,-1}\,\Bigr]_{\,\P}\ =\ 0\,.
\end{equation}
The last equation will be solved numerically in the following Section.

\subsection{Numerical discretization}

Earlier we derived Equation~\eqref{eq:lagr}, which governs the dynamics of the pseudo-inverse mapping $\mathscr{X}\,(\P,\,\cdot)\,$. The initial condition for Equation~\eqref{eq:lagr} is given by the pseudo-inverse \eqref{eq:initp} of the initial condition $\P_{\,0}\,(x)\,$. We discretize Equation~\eqref{eq:lagr} with an explicit discretization in time since it yields the most straightforward implementation.

The first step in our algorithm consists in choosing the initial sampling interval. We make this choice depending on the provided initial condition. Typically, we want to sample only where it is needed. Thus, it seems reasonable to choose the \emph{initial} segment $[\,0,\,\ell_{\,0}\,]$ with $\ell_{\,0}$ being the leftmost location such that
\begin{equation*}
  1\ -\ \P_{\,0}\,(\ell_{\,0})\ <\ \mathsf{tol}\,.
\end{equation*}
In simulations presented below we chose $\mathsf{tol}\ \sim\ \O\,(10^{\,-5})\,$. Then, we choose the initial sampling $\bigl\{\mathscr{X}_{\,k}^{\,0}\bigr\}_{\,k\,=\,0}^{\,N}\ \in\ [\,0,\,\ell_{\,0}\,]\ \subseteq\ \mathds{R}^{\,+}\,$, with $\mathscr{X}_{\,0}^{\,0}\ =\ 0$ and $\mathscr{X}_{\,N}^{\,0}\ =\ \ell_{\,0}\,$. It is desirable that the initial sampling be adapted to the initial condition, since errors made initially cannot be corrected later. We define also $\P_{\,k}\ =\ \P_{\,0}\,(\mathscr{X}_{\,k}^{\,0})\,$. We stress out that $\bigl\{\P_{\,k}\bigr\}_{\,k\,=\,0}^{\,N}$ stand for a discrete cumulative mass variable and, thus, they are time independent.

More generally, we introduce the following notation:
\begin{equation*}
  \mathscr{X}_{\,k}^{\,n}\ \mathop{\stackrel{\,\mathrm{def}}{:=}\,}\ \mathscr{X}\,(\P_{\,k},\,t^{\,n})\,, \qquad k\ =\ 0,\,1,\,\ldots,\,N\,,
\end{equation*}
with $t^{\,n}\ \mathop{\stackrel{\,\mathrm{def}}{:=}\,}\ n\,\Delta t\,$, $n\ \in\ \mathds{N}$ and $\Delta t\ >\ 0$ is a chosen time step\footnote{We present our algorithm with a constant time step for the sake of simplicity. However, in realistic simulations presented in Section~\ref{sec:num} the time step will be chosen adaptively and automatically to meet the stability and accuracy requirements prescribed by the user.}. We introduce also similar notation for the dynamic variable:
\begin{equation*}
  \mathscr{Y}_{\,k}^{\,n}\ \mathop{\stackrel{\,\mathrm{def}}{:=}\,}\ \mathscr{Y}\,(\P_{\,k},\,t^{\,n})\,, \qquad \mathscr{Y}_{\,k}^{\,0}\ \equiv\ \mathscr{X}_{\,k}^{\,0}\,, \qquad k\ =\ 0,\,1,\,\ldots,\,N\,.
\end{equation*}
Now we can state the fully discrete scheme for Equation~\eqref{eq:y}:
\begin{equation}\label{eq:fd}
  \frac{\mathscr{Y}_{\,k}^{\,n\,+\,1}\ -\ \mathscr{Y}_{\,k}^{\,n}}{\Delta t}\ +\ \eta\,\mathrm{e}^{\,\gamma\,t^{\,n}}\;\frac{\mathscr{Y}_{\,k}^{\,n}}{\Delta \P_{\,k}}\;\biggl\{\frac{\Delta\P_{k\,+\,{\textstyle{1\over2}}}}{\mathscr{Y}_{\,k\,+\,1}^{\,n}\ -\ \mathscr{Y}_{\,k}^{\,n}}\ -\ \frac{\Delta\P_{k\,-\,{\textstyle{1\over2}}}}{\mathscr{Y}_{\,k}^{\,n}\ -\ \mathscr{Y}_{\,k\,-\,1}^{\,n}}\biggr\}\,,
\end{equation}
with $n\ \geq\ 0\,$, $k\ =\ 0,\,1,\,\ldots,\,N\,-\,1$ and
\begin{equation*}
  \Delta\P_{k\,+\,{\textstyle{1\over2}}}\ \mathop{\stackrel{\,\mathrm{def}}{:=}\,}\ \P_{\,k\,+\,1}\ -\ \P_{\,k}\,, \qquad
  \Delta\P_{k\,-\,{\textstyle{1\over2}}}\ \mathop{\stackrel{\,\mathrm{def}}{:=}\,}\ \P_{\,k}\ -\ \P_{\,k\,-\,1}\,.
\end{equation*}
The quantity $\Delta \P_{\,k}$ can be defined as the arithmetic or geometric mean of two neighbouring discretization steps $\Delta\P_{k\,\pm\,{\textstyle{1\over2}}}\,$:
\begin{equation*}
  \Delta \P_{\,k}\ \mathop{\stackrel{\,\mathrm{def}}{:=}\,}\ \frac{\Delta\P_{k\,+\,{\textstyle{1\over2}}}\ +\ \Delta\P_{k\,-\,{\textstyle{1\over2}}}}{2}\,, \qquad
  \Delta \P_{\,k}\ \mathop{\stackrel{\,\mathrm{def}}{:=}\,}\ \sqrt{\Delta\P_{k\,+\,{\textstyle{1\over2}}}\cdot\Delta\P_{k\,-\,{\textstyle{1\over2}}}}\,.
\end{equation*}
To be specific, in our code we implemented the arithmetic mean. The fully discrete scheme can be easily rewritten under the form of a discrete dynamical system:
\begin{equation*}
  \mathscr{Y}_{\,k}^{\,n\,+\,1}\ =\ \mathscr{Y}_{\,k}^{\,n}\ -\ \eta\,\Delta t\,\mathrm{e}^{\,\gamma\,t^{\,n}}\;\frac{\mathscr{Y}_{\,k}^{\,n}}{\Delta \P_{\,k}}\;\biggl\{\frac{\Delta\P_{k\,+\,{\textstyle{1\over2}}}}{\mathscr{Y}_{\,k\,+\,1}^{\,n}\ -\ \mathscr{Y}_{\,k}^{\,n}}\ -\ \frac{\Delta\P_{k\,-\,{\textstyle{1\over2}}}}{\mathscr{Y}_{\,k}^{\,n}\ -\ \mathscr{Y}_{\,k\,-\,1}^{\,n}}\biggr\}\,, \quad n\ \geq\ 0\,.
\end{equation*}

\begin{remark}
We would like to say a few words about the implementation of boundary conditions. First of all, no boundary condition is required on the left side, where $\mathscr{X}_{\,0}^{\,n}\ =\ \mathscr{Y}_{\,0}^{\,n}\ \equiv\ 0\,$. On the right boundary we prefer to impose the homogeneous \textsc{Neumann}-type boundary condition, which yields the exact `mass' conservation at the discrete level as well. Namely, at the rightmost cell we have the following fully discrete scheme:
\begin{equation*}
  \mathscr{Y}_{\,N}^{\,n\,+\,1}\ =\ \mathscr{Y}_{\,N}^{\,n}\ +\ \eta\,\Delta t\,\mathrm{e}^{\,\gamma\,t^{\,n}}\;\frac{\mathscr{Y}_{\,N}^{\,n}}{\Delta \P_{\,N}}\cdot\frac{\Delta\P_{N\,-\,{\textstyle{1\over2}}}}{\mathscr{Y}_{\,N}^{\,n}\ -\ \mathscr{Y}_{\,N\,-\,1}^{\,n}}\,, \qquad n\ \geq\ 0\,,
\end{equation*}
with $\Delta \P_{\,N}\ \mathop{\stackrel{\,\mathrm{def}}{:=}\,}\ \dfrac{\P_{\,N}\ -\ \P_{\,N\,-\,2}}{2}\,$. As a result, we obtain the exact conservation of `mass' at the discrete level:
\begin{equation*}
  \sum_{k\,=\,0}^{\,N}\,\Delta \P_{\,k}\,\mathscr{X}_{\,k}^{\,n}\ \equiv\ \sum_{k\,=\,0}^{\,N}\,\Delta \P_{\,k}\,\mathscr{X}_{\,k}^{\,0}\,, \qquad \forall n\ \in\ \mathds{N}\,.
\end{equation*}
\end{remark}

To summarize, our numerical strategy consists in:
\begin{enumerate}
  \item We compute the pseudo-inverse of the initial data $p_{\,0}\,(x)$ to obtain $\mathscr{X}\,(\P,\,0)\ \equiv\ \mathscr{Y}\,(\P,\,0)\,$.
  \item This initial condition $\mathscr{Y}\,(\P,\,0)$ is evolved in (discrete) time using an \emph{explicit} marching scheme in order to obtain numerical approximation to $\mathscr{Y}\,(\P,\,t)\,$, $t\ >\ 0\,$.
  \item The variable $\mathscr{X}\,(\P,\,t)$ is recovered by inverting \eqref{eq:xy}, \emph{i.e.}\xspace
  \begin{equation*}
    \mathscr{X}\,(\P,\,t)\ =\ \mathrm{e}^{\,-\gamma\,t}\,\mathscr{Y}\,(\P,\,t)\,.
  \end{equation*}
  \item Thanks to \eqref{eq:78} we can deduce the values of $p\,\bigl(\mathscr{X}(\P,\,t),\,t\bigr)\ \in\ [\,0,\,1\,]\,$ by applying a favourite finite difference formula\footnote{In our code we employed the simplest forward finite differences and it lead satisfactory results. This point can be easily improved when necessary.}.
\end{enumerate}
Working with the pseudo-inverse allows to overcome the issue of the retention phenomenon, which manifests as the expanding support of $p\,(x,\,t)$ for positive (and possibly large) times $t\ >\ 0\,$, $t\ \gg\ 1\,$, since the computational domain was transformed to $[\,0,\,1\,]\,$. This method is the \textsc{Lagrangian} counterpart of the moving mesh technique in the \textsc{Eulerian} setting \cite{Khakimzyanov2015a, Khakimzyanov2015b}.

A simple \textsc{Matlab} code, which implements the scheme we described hereinabove, is freely available for reader's convenience at this URL address:
\smallskip

\url{https://github.com/dutykh/Feller/}

\section{Numerical results}
\label{sec:num}

In this Section we validate and illustrate the application of the proposed algorithm on several examples. However, first we begin with a straightforward validation test. The only difference with the proposed above algorithm is that we are using a higher order adaptive time stepping for our practical simulations. The explicit first order scheme was used to simplify the presentation. In practice, much more sophisticated time steppers can be used. For instance, we shall employ the explicit embedded \textsc{Dormand}--\textsc{Prince} \textsc{Runge}--\textsc{Kutta} pair $(4,\,5)$ \cite{Dormand1980} implemented in \textsc{Matlab} under the \texttt{ode45} routine \cite{Shampine1997}. Conceptually, this method is similar to the explicit \textsc{Euler} scheme presented above. It conserved all good properties we showed, but it provides additionally the higher accuracy order and totally automatic adaptivity of the time step, which matches very well with adaptive features of the \textsc{Lagrangian} scheme in space. The values of absolute and relative tolerances used in the time step choice are systematically reported below.

\subsection{Steady state preservation}
\label{sec:wb}

In order to validate the numerical algorithm we have at our disposal a family of steady state solutions \eqref{eq:steady}. Hence, if we take such a solution as an initial condition, normally the algorithm has to keep it intact under the discretized dynamics. The parameters $\eta\,$, $\gamma$ of the equation, those of the steady solution and numerical parameters used in our computation are reported in Table~\ref{tab:params0}. The initial condition at $t\ =\ 0$ along with the final state at $t\ =\ T$ are shown in Figure~\ref{fig:steady}. Up to graphical resolution they coincide completely. We can easily check that during the whole simulation the points barely moved, \emph{i.e.}\xspace
\begin{equation*}
  \norm{\mathscr{X}(\cdot,\,T)\ -\ \mathscr{X}(\cdot,\,0)}_{\,\infty}\ \approx\ 0.008577\ldots
\end{equation*}
We can check also other quantities. For instance, $\P\,(\text{\textxi},\,t)$ is preserved up to the machine precision. If we reconstruct the probability distribution $p\,(x,\,t)\,$, we obtain:
\begin{equation*}
  \norm{p\,(\cdot,\,T)\ -\ p\,(\cdot,\,0)}_{\,\infty}\ \approx\ 0.003051\ldots
\end{equation*}
The last error comes essentially from the fact that we apply simple first order finite difference to reconstruct the variable $p\,(x,\,t)$ from its primitive $\P\,(x,\,t)\,$. We can improve this point, but even this simple reconstruction seems to be consistent with the overall scheme accuracy. Thus, this example shows that our implementation of the proposed algorithm is also practically well-balanced \cite{Gosse2013}, since steady state solutions are well preserved.

\begin{table}
  \begin{tabular}{l|c}
  \hline\hline
  \textit{Parameter} & \textit{Value} \\
  \hline\hline
  Drift coefficient, $\gamma$ & $1.0$ \\
  Diffusion coefficient, $\eta$ & $1.0$ \\
  Integration constant, $\mathscr{C}_{\,1}$ & $0.0$ \\
  Integration constant, $\mathscr{C}_{\,2}$ & $\dfrac{\eta}{\gamma}\ \equiv\ 1.0$ \\
  Final simulation time, $T$ & $10.0$ \\
  Number of discretization points, $N$ & $500$ \\
  Absolute tolerance, $\mathsf{tol}_{\,\mathrm{a}}$ & $10^{\,-5}$ \\
  Relative tolerance, $\mathsf{tol}_{\,\mathrm{r}}$ & $10^{\,-5}$ \\
  \hline\hline
  \end{tabular}
  \bigskip
  \caption{\small\em Numerical parameters used in the steady state computations.}
  \label{tab:params0}
\end{table}

\begin{figure}
  \centering
  \includegraphics[width=0.59\textwidth]{figs/Steady.eps}
  \caption{\small\em Comparison of a steady state solution of class \eqref{eq:steady} at $t\ =\ 0$ and at $t\ =\ T\ =\ 10\,$. They are indistinguishable up to the graphical resolution, which validates the solver.}
  \label{fig:steady}
\end{figure}

\subsection{Transient computations}

In this Section we present a couple of extra truly unsteady computations in order to illustrate the capabilities of our method. Namely, we shall simulate the probability distributions emerging from a family of initial conditions (normalized to have the probability distribution):
\begin{equation*}
  p_{\,0}\,(x)\ =\ \frac{\mathrm{e}^{\,-\frac{x}{\sigma_{\,1}}}\ +\ \mathrm{e}^{\,-\frac{x\ -\ x_{\,0}}{\sigma_{\,2}}}}{\sigma_{\,1}\ +\ \sigma_{\,2}\,\mathrm{e}^{\frac{x_{\,0}}{\sigma_{\,2}}}}\,.
\end{equation*}
The primitive of the last distribution can be easily computed as well:
\begin{equation*}
  \P\,(x)\ =\ 1\ -\ \frac{\sigma_{\,1}\,\mathrm{e}^{\,-\frac{x}{\sigma_{\,1}}}\ +\ \sigma_{\,2}\,\mathrm{e}^{\,-\frac{x\ -\ x_{\,0}}{\sigma_{\,2}}}}{\sigma_{\,1}\ +\ \sigma_{\,2}\,\mathrm{e}^{\frac{x_{\,0}}{\sigma_{\,2}}}}\,.
\end{equation*}
We design two different experiments \emph{in silico} to show two completely different behaviours of solutions to \textsc{Feller} equation \eqref{eq:feller} depending on the sign of the drift coefficient $\gamma\,$. These will constitute additional tests for the proposed numerical method. In both cases, the initial positions of particles are chosen according to the logarithmic distribution (\texttt{logspace} function in \textsc{Matlab}) on the segment $[\,0,\,20\,]\,$. This choice is made to represent more accurately the exponentially decaying initial condition since the errors made in the initial condition cannot be corrected in the dynamics.

\subsubsection{Expanding drift}

If the drift coefficient $\gamma\ <\ 0\,$, this term will cause \textsc{Feller}'s equation \eqref{eq:feller} to transport information in the rightward direction. This situation is quite uncomfortable from the numerical point of view since the initial condition expands quickly towards the (positive) infinity. We submitted our method to this case. All numerical and physical parameters are provided in Table~\ref{tab:params1} (the middle column). The initial condition along with the probability distribution at the end of our simulation are shown simultaneously in Figure~\ref{fig:expa}(\textit{a,b}) in \textsc{Cartesian} and semi-logarithmic coordinates correspondingly. As expected, the support of the initial condition more than triples from $t\ =\ 0$ to $t\ =\ T\ =\ 3.0\,$. We remind that the diffusion and drift coefficients are proportional to $x$ and the scheme handles the growth of these terms automatically. The smooth decay of the solution on the semi-logarithmic plot (see Figure~\ref{fig:expa}(\textit{b})) shows the absence of any numerical instabilities. The trajectory of grid nodes is shown in Figure~\ref{fig:expa}(\textit{c}). One can see that points follow the expansion of the solution. Nevertheless, they concentrate in the areas where the probability takes significant values.

\begin{table}
  \begin{tabular}{l|c|c}
  \hline\hline
  \textit{Parameter} & \textit{Expanding experiment} & \textit{Confining experiment} \\
  \hline\hline
  Drift coefficient, $\gamma$ & $-0.1$ & $0.5$ \\
  Diffusion coefficient, $\eta$ & $1.0$ & $1.0$ \\
  Final simulation time, $T$ & $3.0$ & $12.0$ \\
  Number of discretization points, $N$ & $100$ & $100$ \\
  Initial condition parameter, $\sigma_{\,1}$ & $2.0$ & $2.0$ \\
  Initial condition parameter, $\sigma_{\,2}$ & $1.0$ & $1.0$ \\
  Initial condition parameter, $x_{\,0}$ & $3.0$ & $3.0$ \\
  Absolute tolerance, $\mathsf{tol}_{\,\mathrm{a}}$ & $10^{\,-5}$ & $10^{\,-5}$ \\
  Relative tolerance, $\mathsf{tol}_{\,\mathrm{r}}$ & $10^{\,-5}$ & $10^{\,-5}$ \\
  \hline\hline
  \end{tabular}
  \bigskip
  \caption{\small\em Numerical parameters used in unsteady computations.}
  \label{tab:params1}
\end{table}

\begin{figure}
  \centering
  \subfigure[]{\includegraphics[width=0.48\textwidth]{figs/Expanding.eps}}
  \subfigure[]{\includegraphics[width=0.48\textwidth]{figs/ToLogExpanding.eps}}
  \subfigure[]{\includegraphics[width=0.65\textwidth]{figs/ExpandingTraj.eps}}
  \caption{\small\em Numerical result of the expanding experiment with negative drift $\gamma\ =\ -0.1\ <\ 0\,$: (a) initial and terminal states of the numerical discretized solution; (b) initial and terminal states of the numerical discretized solution in semi-logarithmic coordinates; (c) trajectories of grid nodes. All numerical parameters for this computation are reported in Table~\ref{tab:params1} (middle column).}
  \label{fig:expa}
\end{figure}

\subsubsection{Confining drift}

In the case of the positive drift coefficient $\gamma\ >\ 0\,$, the \textsc{Feller} dynamics gets even more interesting since we have two competing effects:
\begin{enumerate}
  \item Positive drift moving information towards the origin $x\ =\ 0\,$,
  \item Diffusion spreading solution towards the positive infinity.
\end{enumerate}
It might happen that both effects balance each other and result in a steady solution. Such a simulation is presented in this Section. The numerical and physical parameters are given in Table~\ref{tab:params1} (the right column). The initial condition along with the probability distribution at the end of our simulation are shown simultaneously in Figure~\ref{fig:std}(\textit{a,b}) in \textsc{Cartesian} and semi-logarithmic coordinates respectively. The trajectories of grid nodes are shown in Figure~\ref{fig:expa}(\textit{c}). One can see the rapid initial expansion stage, which is followed by a stabilization process, when we almost achieved the expected steady state. Again, the grid nodes trajectories show excellent adaptive features of the numerical scheme: at the end of the simulation the relative density of nodes turns out to be preserved. The semi-logarithmic plot shown in Figure~\ref{fig:std}(\textit{b}) shows that the numerical solution is free of numerical instabilities. The obtained steady solution will be preserved by discrete dynamics thanks to the well-balanced property demonstrated in Section~\ref{sec:wb}.

\begin{figure}
  \centering
  \subfigure[]{\includegraphics[width=0.47\textwidth]{figs/ToSteady.eps}}
  \subfigure[]{\includegraphics[width=0.49\textwidth]{figs/ToLogSteady.eps}}
  \subfigure[]{\includegraphics[width=0.65\textwidth]{figs/SteadyTraj.eps}}
  \caption{\small\em Numerical result of the confining experiment with positive drift $\gamma\ =\ 0.5\ >\ 0\,$: (a) initial and terminal states of the numerical discretized solution; (b) initial and terminal states of the numerical discretized solution in semi-logarithmic coordinates; (c) trajectories of grid nodes. All numerical parameters for this computation are reported in Table~\ref{tab:params1} (right column).}
  \label{fig:std}
\end{figure}

\section{Discussion}
\label{sec:disc}

Above we presented some rationale about the discretization, existence and uniqueness theory for \textsc{Feller}'s equation. The main conclusions and perspectives are outlined below.

\subsection{Conclusions}

The celebrated \textsc{Feller} equation was studied mathematically in two seminal papers published by William \textsc{Feller} (1951/1952) in Annals of Mathematics \cite{Feller1951, Feller1952}. In particular, the uniqueness was shown in \cite{Feller1952} without any boundary condition required at $x\ =\ 0\,$. This result is notable and we use it in our study. The main goal of our work was to present an efficient numerical scheme, which were able to handle the unbounded diffusion inherent to Equation~\eqref{eq:feller}. Moreover, the retention phenomenon causes the solution support to expand all the time. Thus, if we wait sufficiently long time, it will reach the computational domain boundaries (since $\mathds{R}^{\,+}\ \ni\ x$ was truncated to a finite segment to make the simulation \emph{in silico} possible). To overcome this difficulty, we proposed a simple and explicit \textsc{Lagrangian} scheme using the notion of the pseudo-inverse. In this way, we obtained a stable numerical scheme (under not so stringent stability conditions), which turns out to be naturally \emph{adaptive} as well, since particles move together with the growing support (the rightmost particle position defines the support) and they tend to concentrate where it is really needed. The performance of our scheme was illustrated on several examples. We share also the \textsc{Matlab} code so that anybody can check the claims we made hereinabove and use it to solve numerically the \textsc{Feller} equation in their applications:
\smallskip

\url{https://github.com/dutykh/Feller/}

\subsection{Perspectives}

All the numerical schemes and results presented in this paper were in one `spatial' dimension. The \textsc{Feller} equation considered here is $1-$D as well. However, it seems interesting\footnote{Regardless the physical applications.} to consider generalized \textsc{Feller} equations in higher space dimensions and to extend the proposed numerical strategy to the multi-dimensional case as well. Another possible extension direction consists in proposing higher order schemes to capture numerical solutions with better accuracy with the same number of degrees of freedom. On the more theoretical side, we would like to obtain an alternative well-posedness theory for \textsc{Feller} equation by taking a continuous limit in our numerical scheme.

\subsection*{Acknowledgments}
\addcontentsline{toc}{subsection}{Acknowledgments}

The author would like to thank Prof.~Laurent~\textsc{Gosse} (CNR--IAC Rome, Italy) for his invaluable help in understanding his first book. We thank also Prof.~Francesco~\textsc{Fedele} (Georgia Tech, Atlanta, USA) for pointing to us this interesting problem.

\bigskip