# Self-similarity in turbulence and its applications

Koji Ohkitani
Research Institute for Mathematical Sciences,
Kyoto University, Kyoto 606-8502 Japan

###### Abstract

First, we discuss the non-Gaussian type of self-similar solutions to the
Navier-Stokes equations. We revisit a class of self-similar solutions which
was studied in Canonne-Planchon (1996). In order to shed some light on it, we
study self-similar solutions to the 1D Burgers equation in detail, completing
the most general form of similarity profiles that it can possibly possess. In
particular, on top of the well-known source-type solution we identify a kink-
type solution. It is represented by one of the confluent hypergeometric
functions, viz. Kummer’s function $M.$

For the 2D Navier-Stokes equations, on top of the celebrated Burgers vortex we
derive yet another solution to the associated Fokker-Planck equation. This can
be regarded as a ’conjugate’ to the Burgers vortex, just like the kink-type
solution above. Some asymptotic properties of this kind of solution have been
worked out. Implications for the 3D Navier-Stokes equations are suggested.

Second, we address an application of self-similar solutions to explore more
general kind of solutions. In particular, based on the source-type self-
similar solution to the 3D Navier-Stokes equations, we consider what we could
tell about more general solutions.

Burgers equation, Navier-Stokes equations, self-similarity

##  1 Introduction

Self-similarity is an useful concept in handling partial differential
equations particularly arising from fluid mechanics. Our motivation for this
study is as follows.

An initial-boundary-value problem for the 3D Navier-Stokes equations is
studied in [1] postulating that the velocity field has a self-similar form,
from the very initial data to the final state. To this end they introduced
Besov spaces to accommodate a singular initial velocity field which is as
rough as $u(x)\propto 1/x$. The knock-on effect is that we would have
similarity profiles, viz. steady solutions to the scaled Navier-Stokes
equations, which are not well-localised in space. In order to shed some light
on its implications, we will take up a simpler system, the 1D Burgers equation
to consider solutions in such an enlarged function space. In particular, we
will be concerned with the fate of solutions which are not in
$L^{1}$marginally, but in some Besov space nearby.

There are two schematic themes in this paper. One is the study of non-Gaussian
type solutions (loosely called ’kink-type’ solutions) of scaled version of
fluid dynamical equations. The other one is proposal of a new approach of
constructing general solutions on the basis of particular source-type self-
similar solutions.

For the first theme, non-Gaussian solutions are given explicitly with some of
their properties discussed, whereas clarification of their possible
significance as a dynamical system is left for further study. For the second
theme a protocol for (re)building more general solutions out of self-similar
solutions is exemplified using the Burgers and 2D Navier-Stokes equations.
This suggests that applications to the 3D Navier-Stokes equations deserve
further investigation.

This paper is organised as follows. We study self-similar profiles for the
Burgers equation in Section 2, emphasising the properties of the newly
identified kink-type solutions. We consider self-similar solutions to the
Navier-Stokes equations, especially in two dimensions in Section 3. We address
possible applications of those similar solutions, in particular about
obtaining information regarding more general class of solutions in Section 4.
Section 5 is devoted to summary and outlook.

##  2 Self-similar solutions of the Burgers equation

When it is demanded that the initial data themselves are self-similar,
inevitably the initial velocity field goes singular like $u(x)\propto 1/x$.
This has difficulties both at the origin and at infinity: (1) it is singular
and non-integrable at the origin and (2) its decay at far distances is too
slow to be integrable. For those reasons, Besov spaces111They are associated
with a norm defined by finite-differences;
$\|\bm{u}\|_{B^{s}_{pq}}\equiv\left\\{\sum_{j=1}^{\infty}\left(2^{sj}\|\Delta_{j}(\bm{u})\|_{L^{p}}\right)^{q}\right\\}^{1/q},$
where $1\leq p,q\leq\infty,s\in\mathbb{R}$ and $\Delta_{j}(\bm{u})$ represents
the band-filtered velocity at frequency $2^{j}$. were introduced to construct
solutions to the initial-boundary problem in [1].

Assume the initial data is in a Besov space $B^{0}_{3,\infty}(\mathbb{R}^{3})$
but not in $L^{3}(\mathbb{R}^{3})$. The construction is concerned with a class
of solutions for $t>0,$ which are in $B^{0}_{3,\infty}(\mathbb{R}^{3})$ but
not in $L^{3}(\mathbb{R}^{3})$. The ’source-type’ self-similar solutions are
irrelevant here, because they are well-localised with a finite $L^{3}$-norm.
The obvious question is: which functions can the scaled solutions possibly
tend to, if they approach steady solutions at all ?

To shed some light on solutions constructed in Besov spaces we consider a
simpler problem of the 1D Burgers equation. For this purpose, instead of
$B^{0}_{3,\infty}(\mathbb{R}^{3})$ we consider velocity fields in
$B^{0}_{1,\infty}(\mathbb{R}^{1})$ to study an analogous problem in one
spatial dimension. It is readily checked that the theory developed in [1]
holds valid mutatis mutandis for the 1D Burgers equation. For more recent
references on self-similar solutions, including studies with use of BMO-type
spaces, see for example [2, 3, 4].

###  2.1 Cole-Hopf transform as Riccati substitution

We consider the Burgers equation in $\mathbb{R}^{1}$

$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial
x}=\nu\frac{\partial^{2}u}{\partial x^{2}},$ (1)

for an initial data $u(x,0)=u_{0}(x)$. It is well-known that (1) is
linearisable by the so-called Cole-Hopf transform. We recall the basic results
obtained in [5].

  * •

If $\int_{0}^{x}u_{0}(y)dy=o(x^{2})$ for large $|x|,$ unique existence of
solutions is guaranteed for all $t>0$.

  * •

If $\int_{0}^{x}u_{0}(y)dy=O(x^{2})$ for large $|x|,$ the existence of
solutions is guaranteed only for finite time $0\leq t<T$.

An instructive example for the second result is given by
$u(x,t)=\frac{x}{t-T}$ [5].

To discuss forward self-similar solutions we introduce dynamic scaling
transformations

$\xi=\frac{x}{\lambda(t)},\;\tau={\color[rgb]{0,0,0}\frac{1}{2a}\log t},$
$u(x,t)=\frac{1}{\lambda(t)}U(\xi,\tau),$

where $\lambda(t)=\sqrt{2at}$ denotes a scaling length and $a(>0)$ a zoom-in
parameter. Applying them to (1) we find the scaled form of the Burgers
equation

$\frac{\partial U}{\partial\tau}+U\frac{\partial
U}{\partial\xi}=a\frac{\partial}{\partial\xi}\left(\xi
U\right)+\nu\frac{\partial^{2}U}{\partial\xi^{2}}.$ (2)

If a steady solution is established in the self-similar variables as
$\tau\to\infty,$ it satisfies

$\frac{d^{2}U}{d\xi^{2}}+\frac{a}{\nu}\frac{d}{d\xi}\left(\xi
U\right)=\frac{d}{d\xi}\frac{U^{2}}{2\nu}.$

After an integration we find the following form

$\frac{dU}{d\xi}+\frac{a}{\nu}\xi U=\frac{U^{2}}{2\nu}+2C,$ (3)

where the prefactor of $2$ is inserted in front of a constant $C$ for
subsequent convenience. (In this paper $c,C,C_{1},\ldots$ etc. denote
constants which may be different from line to line.)

We now distinguish two cases.

  1. 1.

If $C=0,$ the problem is well-understood. It can be solved by either
introducing an integrating factor $\frac{dU}{d\xi}+\frac{a}{\nu}\xi
U=e^{-\frac{a\xi^{2}}{2\nu}}\frac{d}{d\xi}\left(Ue^{\frac{a\xi^{2}}{2\nu}}\right),$
or linearising with a transformation $V=1/U$ regarding it as a Bernoulli
equation. Either way we find the so-called source-type solution, e.g. [6, 7]

$U(\xi)=\frac{C^{\prime}\displaystyle{\exp\left(-\frac{a\xi^{2}}{2\nu}\right)}}{1-\displaystyle{\frac{C^{\prime}}{2\nu}\int_{0}^{\xi}}\exp\left(-\frac{a\eta^{2}}{2\nu}\right)d\eta},$
(4)

for a constant $C^{\prime}$. The solution is well-localised spatially, see
Figs.2, 2. The name has come from the fact that $\lim_{t\to
0}\frac{1}{\sqrt{2at}}U(\xi)=M\delta(x),$ where $\delta(\cdot)$ denotes the
Dirac mass and $M\equiv\int_{-\infty}^{\infty}u_{0}(x)dx.$ We also recall [6,
7] that for all $u_{0}\in L^{1}(\mathbb{R}^{1})$ and $1\leq p\leq\infty,$

$t^{\frac{1}{2}\left(1-\frac{1}{p}\right)}\left\|u(x,t)-\frac{1}{\sqrt{2at}}U(\xi)\right\|_{L^{p}}\to
0\;\;\mbox{as}\;\;t\to\infty,$

showing some degree of the universality of the profile.

  2. 2.

If $C\neq 0,$ (3) is a Riccati equation and the method above does not work. In
this case, we ought to introduce a Riccati substitution
$U=-2\nu\dfrac{\Psi^{\prime}}{\Psi},$ to reduce it to a linear equation (see
Appendix A)

$\Psi^{\prime\prime}+\frac{a}{\nu}\xi\Psi^{\prime}+\frac{C}{\nu}\Psi=0.$ (5)

It is to be observed that the Cole-Hopf transform arises as a natural course
of solution to the Riccati equation [8]. If the initial data $u_{0}\notin
L^{1}(\mathbb{R}^{1}),$ but $u_{0}\in B^{0}_{1,\infty}(\mathbb{R}^{1}),$ the
other kind of solutions, i.e. the kink-type ones will come into play. Below we
will discuss those solutions to (5) in detail.

###  2.2 Kink-type solution

As noted already, the theory developed in [1] for the 3D Navier-Stokes
equations works for initial velocity $\notin
L^{3}(\mathbb{R}^{{\color[rgb]{0,0,0}3}})$, but $\in
B^{0}_{{\color[rgb]{0,0,0}3},\infty}(\mathbb{R}^{3}).$ We will consider an
analogous 1D problem as an illustration.

In fact, the equation (5) can be solved using confluent hypergeometric
functions, viz. the Kummer’s functions. Consider

$\Psi^{\prime\prime}+A\xi\Psi^{\prime}+B\Psi=0,$ (6)

where $A$ and $B$ are constants (to be set $A=\frac{a}{\nu}$ and
$B=\frac{C}{\nu}$). Actually we have the following

###### Proposition 2.1

The solution to (6) can be written

$\Psi=\xi
e^{-\frac{A}{2}\xi^{2}}w\left(1-\frac{B}{2A},\frac{3}{2},\frac{A}{2}\xi^{2}\right),$

where the function $w\left(1-\frac{B}{2A},\frac{3}{2},z\right)$ satisfies the
confluent hypergeometric equation

$z\frac{d^{2}w}{dz^{2}}+\left(\frac{3}{2}-z\right)\frac{dw}{dz}+\left(\frac{B}{2A}-1\right)w=0.$
(7)

Hence the general solution is given by

$\Psi=C_{1}\xi\exp\left(-\frac{a}{2\nu}\xi^{2}\right)M_{\rm
K}\left(1-\frac{C}{2a},\frac{3}{2},\frac{a}{2\nu}\xi^{2}\right)+C_{2}\xi\exp\left(-\frac{a}{2\nu}\xi^{2}\right)U_{\rm
K}\left(1-\frac{C}{2a},\frac{3}{2},\frac{a}{2\nu}\xi^{2}\right),$

where $M_{\rm K}(\alpha,\gamma,z)$ and $U_{\rm K}(\alpha,\gamma,z)$ denote two
fundamental solutions222Standard notations for Kummer’s functions are
$M(\alpha,\gamma,z)$ and $U(\alpha,\gamma,z).$ We add the subscript K to avoid
confusion with the scaled velocity $U$. of (7).

We note that the real parameter $A$ does not have to be positive in the
general solutions above. This will be important when we consider a backward
self-similar solution below.

Proof
This is done by straightforward calculations. Taking
$A{\color[rgb]{0,0,0}(=a/\nu)}=1$ without loss of generality, consider
$\Psi=\xi e^{-\frac{\xi^{2}}{2}}w\left(1-\dfrac{B}{2},\dfrac{3}{2},z\right),$
with $z=\xi^{2}/2.$ Direct calculations show

$\partial_{\xi}\Psi=(1-\xi^{2})e^{-\frac{\xi^{2}}{2}}w+\xi^{2}e^{-\frac{\xi^{2}}{2}}\partial_{z}w,$

and

$\partial_{\xi\xi}\Psi=(\xi^{3}-3\xi)e^{-\frac{\xi^{2}}{2}}w+(3\xi-2\xi^{3})e^{-\frac{\xi^{2}}{2}}\partial_{z}w+\xi^{3}e^{-\frac{\xi^{2}}{2}}\partial_{zz}w.$

Thus we get

$\displaystyle\partial_{\xi\xi}\Psi+\xi\partial_{\xi}\Psi+B\Psi$
$\displaystyle=$ $\displaystyle
e^{-\frac{\xi^{2}}{2}}\left\\{\xi^{3}\partial_{zz}w+(3\xi-\xi^{3})\partial_{z}w+(B-2)\xi
w\right\\}$ $\displaystyle=$ $\displaystyle
e^{-\frac{\xi^{2}}{2}}\xi\left\\{\xi^{2}\partial_{zz}w+(3-\xi^{2})\partial_{z}w+(B-2)w\right\\}.$

As $\xi^{2}=2z,$ we deduce

$z\frac{d^{2}w}{dz^{2}}+\left(\frac{3}{2}-z\right)\frac{dw}{dz}+\left(\frac{B}{2}-1\right)w.\;\;\square$

We need to check which option, $M_{\rm K}$ or $U_{\rm K},$ is acceptable for
our purpose. First we check $M_{\rm K}$. By the asymptotic formulas for
$|z|\to\infty$ $(\Re z>0)$ (see Appendix B), we have as $|\xi|\to\infty$

$M_{\rm
K}\left(1-\frac{B}{2A},\frac{3}{2},\frac{A\xi^{2}}{2}\right)\approx\frac{\Gamma(\frac{3}{2})}{\Gamma(1-\frac{B}{2})}e^{A\xi^{2}/2}\left(\frac{A\xi^{2}}{2}\right)^{-\frac{1}{2}(1+\frac{B}{A})}\propto
e^{A\xi^{2}/2}\xi^{-(1+\frac{B}{A})}.$

Hence with this choice we obtain

$\Psi\approx\xi^{-\frac{B}{A}}\;\;\mbox{and}\;\;U=-2\nu\partial_{\xi}\log\Psi\propto\xi^{-1},$

which is consistent with the boundary behaviour, i.e. the condition $C\neq 0$
above.

On the other hand, for the other solution $U_{\rm K}$ we have as
$|\xi|\to\infty$

$U_{\rm
K}\left(1-\frac{B}{2A},\frac{3}{2},\frac{A\xi^{2}}{2}\right)\approx\left(\frac{A\xi^{2}}{2}\right)^{\frac{B}{2A}-1}\propto\xi^{\frac{B}{A}-2},$

thus in this case we find

$\Psi\propto
e^{-A\xi^{2}/2}\xi^{\frac{B}{A}-1}\;\mbox{and}\;U\propto-A\xi+\left(\frac{B}{A}-1\right)\xi^{-1}.$

We should discard this option because, when $A\neq 0$, the corresponding $U$
does not even belong to $B^{0}_{1,\infty}(\mathbb{R}^{1}),$ due to the
presence of the linear term in $\xi$.

We conclude that for the kink solution, we should choose

$\Psi=C_{1}\xi\exp\left(-\frac{a}{2\nu}\xi^{2}\right)M_{\rm
K}\left(1-\frac{C}{2a},\frac{3}{2},\frac{a}{2\nu}\xi^{2}\right).$

It may be in order to have a look at a specific example of the class of
solutions. Replacing $z\to iz,$ in the following identity, e.g. [9],

$M_{\rm K}\left(\frac{1}{2}\frac{3}{2},-z^{2}\right)=\frac{\sqrt{\pi}}{2z}{\rm
erf}(z),$

we have

$M_{\rm K}\left(\frac{1}{2}\frac{3}{2},z^{2}\right)=\frac{\sqrt{\pi}}{2iz}{\rm
erf}(iz)=\frac{\sqrt{\pi}}{2z}{\rm erfi}(z).$

Here ${\rm erf}(z)=\frac{2}{\sqrt{\pi}}\int_{0}^{z}e^{-t^{2}}dt$ denotes the
error function and its imaginary version ${\rm
erfi}(z)=\frac{2}{\sqrt{\pi}}\int_{0}^{z}e^{t^{2}}dt.$ In order to make use of
the identity, we take $C/a=1$. With this a typical example can be given

$\Psi={\color[rgb]{0,0,0}C_{1}}\xi\exp\left(-\frac{a}{2\nu}\xi^{2}\right)M_{\rm
K}\left(\frac{1}{2}\frac{3}{2},\frac{a\xi^{2}}{2\nu}\right)={\color[rgb]{0,0,0}C_{1}}\sqrt{\frac{2\nu}{a}}D\left(\sqrt{\frac{a}{2\nu}}\xi\right),$

where $D(x)\equiv
e^{-x^{2}}\int_{0}^{x}e^{t^{2}}dt=\frac{\sqrt{{\color[rgb]{0,0,0}\pi}}}{2}H[e^{-x^{2}}]$
denotes Dawson’s integral and $H[\cdot]$ the Hilbert transform. See Figs. 4,
4, in which we also assume $a=1$ and $\nu=1/2$ for simplicity. Note that the
profile in Fig.4 is a steady solution in scaled space, which develops from the
initial profile in the original variables in Fig.3 under time evolution.

It can be seen that in higher spatial dimensions (6) is generalised to

$\triangle\Psi+\frac{a}{\nu}(\bm{\xi}\cdot\nabla)\Psi+\frac{C}{\nu}\Psi=0.$
(8)

However, representation formulas for solutions to (8) are not known.

See Appendix C for a motivation or rationale for studying two different kinds
of steady solutions.

###  2.3 Backward self-similar solutions

Generally speaking, we talk about forward self-similar solutions to study the
decaying process in the late stage of evolution, whereas talk about backward
self-similar solutions to study whether solutions blow up in finite time.
Needless to mention, no solutions to the Burgers equation blow up if they
start from well-localised smooth initial data under natural boundary
conditions. Nonetheless, because the backward problem differs from the forward
problem only in the sign of the parameter $a,$ if one of them is obtained it
is readily transferable to the other one by flipping the sign of $a,$ provided
that we do not bother their boundary behaviour.

In the spirit of [1] we have studied forward self-similar solutions to the
Burgers equation under the setting of $B^{0}_{1,\infty}(\mathbb{R}^{1}).$ What
would happen if we consider the backward self-similar solutions in the
enlarged function class, or even broader one ?

We can study ’possible’ blow up with backward self-similarity, putting

$u(x.t)=\frac{1}{\lambda(t)}U\left(\frac{x}{\lambda(t)}\right)$

with the length-scale $\lambda(t)=\sqrt{2a(t_{*}-t)},$ where $t_{*}$ denotes
the time of blowup. The steady equation reads

$U\frac{\partial U}{\partial\xi}+a\left(\xi\frac{\partial
U}{\partial\xi}+U\right)=\nu\frac{\partial^{2}U}{\partial\xi^{2}},$

to which the only smooth solution is a trivial one $U\equiv 0$ under natural
boundary conditions and smoothness conditions. Relaxing those conditions, a
non-trivial solution is nonetheless obtained as

$U(\xi)=\frac{C^{\prime}e^{a\xi^{2}/(2\nu)}}{1-\frac{C^{\prime}}{2\nu}\int_{0}^{\xi}e^{a\eta^{2}/(2\nu)}d\eta}.$

This solution is badly-behaved at far distances;333In Appendix D of [12] it
was stated erroneously that
$U(\xi)\to\frac{1}{|\xi|}\;\;\mbox{as}\;\;\xi\to\pm\infty,$ which should be
corrected as above. $U(\xi)\propto\xi$ as $\xi\to\pm\infty.$ Furthermore this
has a singular point (a pole) somewhere, say at $\xi=\xi_{*},$ where the
denominator vanishes. Thus we have $U(\xi)\propto 1/(\xi-\xi_{*})$ around it
and is non-integrable $U\notin L^{1}_{\rm loc}.$ See Fig.5.

It is of interest to have another look at the backward self-similar solution.
We first recast the source-type solution using Kummer’s function. Taking
$\beta=1/2$ in the following identity [9]

$e^{-z}M_{\rm K}(1,\beta+1,z)=\beta z^{-\beta}\gamma(\beta,z),$

where $\gamma(\beta,z)\equiv\int_{0}^{z}t^{\beta-1}e^{-t}dt,$ we have

$e^{-z}M_{\rm
K}\left(1,\frac{3}{2},z\right)=z^{-1/2}\int_{0}^{\sqrt{z}}e^{-s^{2}}ds.$

Putting $z=\frac{a}{2\nu}\xi^{2},$ we find

$\Psi=\xi\exp\left(-\frac{a}{2\nu}\xi^{2}\right)M_{\rm
K}\left(1,\frac{3}{2},\frac{a}{2\nu}\xi^{2}\right)=\sqrt{\frac{\pi\nu}{2a}}{\rm
erf}\left(\sqrt{\frac{a}{2\nu}}\xi\right).$

When we replace $a\to-a,$ the right-hand side is changed as

$\Psi\to\frac{1}{i}\sqrt{\frac{\pi\nu}{2a}}{\rm
erf}\left(i{\color[rgb]{0,0,0}\sqrt{\frac{a}{2\nu}}}\xi\right)=\sqrt{\frac{\pi\nu}{2a}}{\rm
erfi}\left(\sqrt{\frac{a}{2\nu}}\xi\right),$

which agrees with backward self-similar solution obtained above. Note that the
asymptotic formula for $M_{\rm K}(\alpha,\gamma,z)$ for $|z|\to\infty$ does
not hold valid for $\alpha=1,$ because $\Gamma(1)=0.$

It may be interesting to consider the following question. The existence of the
forward self-similar (source-type) solutions $\bm{U}$ to the 3D Navier-Stokes
equations is known. Suppose we make a replacement $a\to-a$ in such solutions,
the backward profile $\bm{U}$ must be singular and/or ill-behaved at far
distances, in view of the non-existence of self-similar blowup [10, 11]. We
would still be interested in what kind of spatial structure the profile
possesses because such a solution may be helpful in putting constraints under
the replacement $a\leftrightarrow-a$ on the sought-after forward source-type
solutions.

Figure 1: The initial data, Dirac delta function $u(x)=\delta(x)$ at $t=0$.
Depicted as
$\delta_{\epsilon}(x)\equiv\exp(-x^{2}/\epsilon)/\sqrt{\pi\epsilon},$ with
$\epsilon=1\times 10^{-4}.$ (The figures in this paper are meant to be
schematic.)

Figure 2: The source-type solution as $\tau\to\infty$;
$U(\xi)=\dfrac{e^{-\xi^{2}}}{1-\int_{0}^{\xi}e^{-\eta^{2}}d\eta.}$

Figure 3: The self-similar initial data $u(x)=\frac{1}{x}$ at $t=0$.

Figure 4: The kink-type solution as $\tau\to\infty$;
$U(\xi)=-\frac{D^{\prime}(\xi)}{D(\xi)}=-\frac{1-2\xi D(\xi)}{D(\xi)}.$

Figure 5: A blowup profile
$U=\dfrac{e^{\xi^{2}}}{1-\int_{0}^{\xi}e^{\eta^{2}}d\eta},$which behaves
$\propto\frac{1}{\xi_{*}-\xi}$ near some $\xi_{*}$ and $\propto\xi$ as
$|\xi|\to\infty.$

##  3 Self-similar solutions of the Navier-Stokes equation

###  3.1 2D Navier-Stokes equations

The 2D Navier-Stokes equations is described by the vorticity equation

$\dfrac{\partial\omega}{\partial
t}+\bm{u}\cdot\nabla\omega=\nu\triangle\omega.$

The dynamically-scaled form of the vorticity equation reads

$\dfrac{\partial\Omega}{\partial\tau}+\bm{U}\cdot\nabla\Omega=\nu\triangle\Omega+a\nabla\cdot(\bm{\xi}\Omega),$
(9)

where $\Omega(\bm{\xi},\tau)=\frac{1}{2at}\omega\left(\bm{\xi},\tau\right)$
denotes the scaled vorticity, $\bm{\xi}=\frac{\bm{x}}{\sqrt{2at}}$ and
$\tau=\frac{1}{2a}\log t$. Under the assumption of radial symmetry, i.e. the
similarity profile $\Omega$ being a function of $\xi=|\bm{\xi}|$ only, it has
a solution

$\Omega(\xi)=\frac{a\Gamma}{2\pi\nu}\exp\left(-\frac{a|\bm{\xi}|^{2}}{2\nu}\right),$
(10)

where $\Gamma=\int_{\mathbb{R}^{2}}\Omega(\bm{\xi})d\bm{\xi}$ denotes the
circulation. This is the celebrated Burgers vortex, which is well-localised in
space and $\in L^{1}(\mathbb{R}^{2})$.

All the radially-symmetric solutions to the 2D Fokker-Planck equations were
worked out explicitly in [13]. One of them is the Burgers vortex above and the
other one is lesser-localised in space. The latter solution reads

$\Omega(\xi)=\frac{a}{4\nu}\exp\left(-\frac{a\xi^{2}}{2\nu}\right)\left({\rm
Ei}\left(\frac{a\xi^{2}}{2\nu}\right)-\log\left(\frac{a}{2\nu}\right)-\gamma\right),$
(11)

where ${\rm Ei}(x)\equiv-{\rm p.v.}\int_{-x}^{\infty}\frac{e^{-t}}{t}dt$
denotes the exponential integral and $\gamma\approx 0.577$ the Euler’s
constant.  Taking $a/2\nu=1$ for simplicity we consider

$\Omega(\xi)=\frac{1}{2}\exp(-\xi^{2})\left({\rm Ei}(\xi^{2})-\gamma\right).$

Recalling ${\rm
Ei}(x)=O\left(\dfrac{e^{x}}{x}\right)\;\;\,\mbox{as}\;\;x\to\infty$ and ${\rm
Ei}(x)\approx\gamma+\log x+x\;\;\,\mbox{as}\;\;x\to 0+,$ we note its
asymptotic properties as follows.

  * •

It decays slowly
$\Omega(\xi)\approx\dfrac{1}{2\xi^{2}}\;\;\mbox{as}\;\;\xi\to\infty.$ Hence
$\Omega(\xi)$ does not belong to $L^{1}(\mathbb{R}^{2})$, but does belong to
the Besov space $B^{0}_{1,\infty}(\mathbb{R}^{2})$.

  * •

Near the origin it shows a mild singularity as
${\color[rgb]{0,0,0}\Omega(\xi)\approx\log\xi}\;\;\mbox{as}\;\;\xi\to 0+,$
that is, it is discontinuous in vorticity at $\xi=0$.

  * •

It can be verified nonetheless that the corresponding velocity is continuous
there. The azimuthal component velocity is given by

${\color[rgb]{0,0,0}U_{\theta}(\xi)=\frac{1}{\xi}\int_{0}^{\xi}\eta\Omega(\eta)d\eta=\frac{1}{4\xi}\left\\{\log\xi^{2}-e^{-\xi^{2}}\left({\rm
Ei}(\xi^{2})-\gamma\right)\right\\}}.$

It has the following asymptotic behaviours:
$U_{\theta}(\xi)\approx\frac{\xi}{2}\log\xi\;\;\mbox{as}\;\;\xi\to
0+\;\mbox{and}\;U_{\theta}(\xi)\approx\frac{1}{2\xi}\log\xi\;\;\mbox{as}\;\;\xi\to\infty.$

Reverting to the original variables from (11) and discarding terms associated
with the Gaussian solution, we obtain a particular solution444 This is same
the procedure by which we get the so-called Lamb-Oseen decaying vortex from
(3.2).

$\omega(r,t)=\frac{1}{2at}\Omega\left(\frac{r}{\sqrt{2at}}\right)=\frac{1}{8\nu
t}\exp\left(-\frac{r^{2}}{4\nu t}\right){\rm Ei}\left(\frac{r^{2}}{4\nu
t}\right).$

We have, for fixed $t$,

$\omega(\bm{x},t)\approx\frac{\log r}{4\nu t}\;\;\mbox{as}\;\;r\to 0+$

by ${\rm Ei}(x)\approx\log x\;\;\mbox{as}\;\;x\to 0+.$ Hence the property of
the vorticity being singular at the origin persists throughout time evolution.
We also have

$\omega(x,t)\approx\frac{1}{2r^{2}}\;\;\mbox{as}\;\;t\to 0+.$

The corresponding initial data is exactly scale-invariant as expected.

The azimuthal velocity field in the original variables is given by

$u_{\theta}(r,t)=\frac{1}{r}\int_{0}^{r}s\omega(s,t)ds$
$=\frac{1}{4r}\left\\{\log\left(\frac{r^{2}}{4\nu
t}\right)-\exp\left(-\frac{r^{2}}{4\nu t}\right){\rm
Ei}\left(\frac{r^{2}}{4\nu t}\right)+\gamma\right\\}.$

Its asymptotic properties are

$u_{\theta}(r,t)\approx\frac{\log r}{2r}\;\mbox{as}\;t\to
0+\;\mbox{and}\;u_{\theta}(r,t)\approx\frac{r}{8\nu t}\log r\;\mbox{as}\;r\to
0+.$

The second non-Gaussian solution may serve as a replacement for the ’kink-
type’ solution to the problem. However, neither realisability nor stability of
(11) is known.

###  3.2 3D Navier-Stokes equations

As already noted, in three dimensions the existence of the forward self-
similar solutions is known, e.g. [1, 14], but the precise functional form of
the solution is not. If the initial condition is well-localised, we know that
it takes a near-Gaussian form if the dependent variable is chosen suitably,
i.e. as the vorticity curl in three-dimensions.

Using the vorticity curl $\bm{\chi}=\nabla\times\bm{\omega}$ the governing
equations read

$\frac{\partial\bm{\chi}}{\partial t}=\triangle(\bm{u}\cdot\nabla\bm{u}+\nabla
p)+\nu\triangle\bm{\chi},$

where $\bm{u}$ denotes the velocity and Under the dynamic scaling
transformations

$\bm{\xi}=\frac{\bm{x}}{\sqrt{2at}},\;\;\tau=\frac{1}{2a}\log
t,\bm{\chi}(\bm{x},t)=\frac{1}{(2at)^{3/2}}\bm{X}(\bm{\xi},\tau),$

we find

$\frac{\partial\bm{X}}{\partial\tau}=\triangle\left(\bm{U}\cdot\nabla\bm{U}+\nabla
P\right)+\nu\triangle\bm{X}+a\nabla\cdot(\bm{\xi}\otimes\bm{X}).$

A perturbative attempt of the determination of source-type solutions can be
found in [13]. The leading-order approximation solution is explicitly given by
a Gaussian function modulo incompressibility and the corrections due to the
nonlinear term is estimated to be small. The conjugate solution at leading-
order is also worked out explicitly.

##  4 Lifting of self-similar solution to more general ones

We discuss what can be learnt from studying self-similar solutions and
identify open problems in this regard.

###  4.1 1D Burgers equation

When we recast (4) as

$U(\xi)=-2\nu\frac{\partial}{\partial\xi}\log\left({1-\displaystyle{\frac{C^{\prime}}{2\nu}\int_{0}^{\xi}}\exp\left(-\frac{a\eta^{2}}{2\nu}\right)d\eta}\right),$
(12)

it is reminiscent of the structure of the celebrated Cole-Hopf transform. In
other words, the source-type solution encodes the vital information as to how
we may linearise this nonlinear equation.

We don’t intend to add anything new to the understanding of the Burgers
equation, rather we will show how we can recover the Cole-Hopf transformation
on the basis of a particular self-similar solution (viz. the source type
solution).

Let us pretend that we don’t know the Cole-Hopf transformation. It’s still
straightforward to derive the source-type solution, which is slightly non-
Gaussian. It is then possible to recover the Cole-Hopf linearisation by
replacing the self-similar heat flow with a general heat flow, which we call
the lifting procedure.

For illustration we will have a look at the details of the procedure. Assume
that the velocity potential $\Phi(\xi)$ is given by

$\Phi(\xi)=f(\widehat{\Psi}),$

where $f$ is a function of the self-similar source type solution of the heat
equation
$\widehat{\Psi}=\int_{0}^{\xi}\exp\left(-\frac{a\eta^{2}}{2\nu}\right)d\eta$.
Then we have

$U=\partial_{\xi}\Phi=f^{\prime}(\widehat{\Psi})\partial_{\xi}\widehat{\Psi}.$

Writing

$U(\xi)=F(\widehat{\Psi};\partial_{\xi}\widehat{\Psi}),$

we have trivially

$F(x;y)=yf^{\prime}(x),$

and the function $F(x;y)$ satisfies the following scaling property

$F(x;\alpha y)=\alpha F(x;y),\;\mbox{for}\;\forall\;\alpha>0.$ (13)

It should be noted that hereafter the arguments $x,y,$ etc. in the function
$F$ do not represent spatial coordinates.

The procedure of spotting (or, recovering, in this case) more general
solutions can be formalised in the following steps.

Step 1. Assume a self-similar profile is obtained in the form

$U(\xi)=F(\widehat{\Psi}(\xi);\partial_{\xi}\widehat{\Psi}(\xi))\equiv\frac{\partial_{\xi}\widehat{\Psi}}{1-\frac{1}{2\nu}\widehat{\Psi}}.$
(14)

The profile $U(\xi)$ is a near-identity transformation of the last argument of
the function, i.e. that of the Gaussian function
$\partial_{\xi}\widehat{\Psi}(\xi)$. Luckily for the Burgers equation
$F(\cdot;\cdot)$ is known explicitly $F(x;y)=\frac{y}{1-x/(2\nu)},$ as
indicated by the symbol $\equiv$ above. By definition, a particular solution
is obtained by reverting to the original variables:

$u(x,t)=\frac{1}{\sqrt{2at}}F\left(\widehat{\Psi}\left(\frac{x}{\sqrt{2at}}\right);\partial_{\xi}\widehat{\Psi}\left(\frac{x}{\sqrt{2at}}\right)\right)\equiv\frac{1}{\sqrt{2at}}\frac{\partial_{\xi}\widehat{\Psi}(\frac{x}{\sqrt{2at}})}{1-\frac{1}{2\nu}\widehat{\Psi}(\frac{x}{\sqrt{2at}})}.$

Step 2. Replacing the self-similar heat solution with the general heat flow,
we obtain a more general class of solutions

$\displaystyle u(x,t)$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2at}}F\left(\widehat{\psi}\left(x,t\right);\sqrt{2at}\partial_{x}\widehat{\psi}\left(x,t\right)\right)$
$\displaystyle=$ $\displaystyle
F\left(\widehat{\psi}\left(x,t\right);\partial_{x}\widehat{\psi}\left(x,t\right)\right)\equiv\frac{\partial_{x}\widehat{\psi}(x,t)}{1-\frac{1}{2\nu}\widehat{\psi}(x,t)},$

where

$\widehat{\psi}(x,t)=\frac{1}{\sqrt{4\pi\nu
t}}\int_{-\infty}^{\infty}\widehat{\psi}(y,0)\exp\left(-\frac{(x-y)^{2}}{4\nu
t}\right)dy$

denotes the general heat flow. Note that the length scales $\sqrt{2at}$ cancel
out because of (13). In the case of the Burgers equation, this last form
provides the general solution.

Before closing this subsection it may be in order to emphasise the following
fact. In terms of the scaled velocity potential $\Phi(\xi)$ the exact solution
is given by $\Phi=-2\nu\log\left(1-\frac{1}{2\nu}\widehat{\Psi}\right),$ whose
leading-order approximation agrees with the scaled heat flow
$\Phi(\xi)\approx\widehat{\Psi}(\xi).$ This confirms the near-identity nature
of $\Phi(\xi)$.

###  4.2 2D Navier-Stokes equations

In 2D incompressible flows the source-type solution is given by the Burgers
vortex, which is not only near-Gaussian, but also exactly Gaussian. Because of
this peculiarity, the final steady state does not contain any useful
information regarding the nonlinear terms. Hence it is impossible to lift (or
generalise) the final state to find a more general class of solutions and we
end up with obtaining the linearised solution only. Curiously enough the
better understood 2D Navier-Stokes equations defy the current approach to gain
some information about their solutions.

Regarding this, it is in order to include a bit more detailed description
based on kinematic relationship. The stream function corresponding to (10),
with $\Gamma=1$, is given by $\widehat{\Psi}(\xi)=\dfrac{1}{4\pi}{\rm
Ei}\left(-\frac{a|\bm{\xi}|^{2}}{2\nu}\right),$ which agrees with the scaled
heat flow in two dimensions. Assume that the scaled stream function
$\Psi(\bm{\xi})$ is given by $\Psi(\bm{\xi})=f(\widehat{\Psi}(\xi)).$ It is
readily derived

$\Omega(\bm{\xi})=-\triangle\Psi=-\triangle\widehat{\Psi}f^{\prime}(\cdot)-|\nabla\widehat{\Psi}(\xi)|^{2}f^{\prime\prime}(\cdot).$

Writing
$\Omega(\bm{\xi})=F\left(\widehat{\Psi}(\xi),\nabla\times\widehat{\Psi}(\xi),-\triangle\widehat{\Psi}(\xi)\right),$
we find

$F(x,y;z)=f^{\prime}(\cdot)z-f^{\prime\prime}(\cdot)y^{2}.$

The function $F$ satisfies the following scaling

$F(x,\alpha
y;\alpha^{2}z)=\alpha^{2}F(x,y;z),\;\mbox{for}\;\forall\;\alpha>0.$

The only solution we know is the identity, that is, $f(x)=x.$

###  4.3 3D Navier-Stokes equations

We take the vorticity curl, $\bm{\chi}=\nabla\times\bm{\omega}$ as the basic
dependent variable, whose dynamically-scaled version is denoted by
$\bm{X}(\bm{\xi}).$ With this choice, the linearised equations have the
Fokker-Planck operator In view of the critical scale-invariance of type 2 [13]
it is most convenient for our analysis, as the leading-order approximation is
basically given by the Gaussian function, i.e. the Gaussian function modulo
incompressibility.

The steady version of the dynamically-scaled Navier-Stokes equations reads
[13]

$\triangle^{*}\bm{X}\equiv\triangle\bm{X}+\frac{a}{\nu}\nabla\cdot(\bm{\xi}\otimes\bm{X})=-\frac{1}{\nu}\triangle\mathbb{P}\left(\triangle^{-1}\bm{X}\cdot\nabla\triangle^{-1}\bm{X}\right).$
(15)

By a formal analysis we show the following

###### Proposition 4.1

The successive approximations to the solution $\bm{X}$ of the equations (15)
are given by a functional 555It is a functional rather than a function as we
need to take into account nonlocal interactions due to the incompressible
condition. of

$\\{(\nabla\times)^{k}\widehat{\bm{\Psi}}(\bm{\xi})\;|\;k=0,1,2,3\\}.$

Here $(\nabla\times)^{3}\widehat{\bm{\Psi}}(\bm{\xi})=\mathbb{P}\bm{M}G,$
defined with $\bm{M}=\int\bm{X}(\bm{\xi})d\bm{\xi},$
$G=\left(\frac{a}{2\pi\nu}\right)^{3/2}\exp(-\frac{a}{2\nu}|\bm{\xi}|^{2})$
and $\mathbb{P}$ solenoidal projection such that
$\nabla\cdot\widehat{\bm{\Psi}}(\bm{\xi})=0$.

Note that $\widehat{\bm{\Psi}}(\bm{\xi})$ denotes the scaled heat flow,
whereas $\bm{\Psi}(\bm{\xi})$ the scaled vector potential.

Proof
We consider the following expansion for small $\epsilon>0,$

$\bm{X}(\bm{\xi})=\epsilon\bm{X}_{1}+\epsilon^{2}\bm{X}_{2}+\epsilon^{3}\bm{X}_{3}+\ldots,$

deferring the justification of smallness of $\epsilon$ by its identification
as the Reynolds number. We prove by mathematical induction that each
$\bm{X}_{n}$ can be represented by a combination of functionals of
$\widehat{\bm{\Psi}}(\bm{\xi}),$ including e.g.
$\triangle^{-1}\widehat{\bm{\Psi}}$ and the derivatives of
$\widehat{\bm{\Psi}}(\bm{\xi})$. Equating the terms with the same powers in
$\epsilon,$ we derive equations for $\bm{X}_{n}$ for $n\geq 1.$

(i) We will first confirm that this is the case for $n=1$.
To leading order at $O(\epsilon)$ we have

$\triangle^{*}\bm{X}_{1}=0,$

from which it follows that $\bm{X}_{1}=\mathbb{P}\bm{M}G.$ Indeed $\bm{X}_{1}$
is a functional of the desired form. By definition
$\bm{X}_{1}=(\nabla\times)^{3}\bm{\Psi}_{1}(\bm{\xi}),$ we also note that the
leading-order approximation satisfies
$\bm{\Psi}_{1}(\bm{\xi})=\widehat{\bm{\Psi}}(\bm{\xi}).$

(ii) Assuming that the statement holds up to step $k(\leq n),$ we will deduce
that the it also holds for step $(n+1).$ For illustration let us first take a
look at, for example, $O(\epsilon^{2})$ and $O(\epsilon^{3})$. To next-to-
leading order at $O(\epsilon^{2})$ we have

$\triangle^{*}\bm{X}_{2}=-\frac{1}{\nu}\triangle\mathbb{P}(\triangle^{-1}\bm{X}_{1}\cdot\nabla\triangle^{-1}\bm{X}_{1})$

and at the third order $O(\epsilon^{3})$

$\triangle^{*}\bm{X}_{3}=-\frac{1}{\nu}\triangle\mathbb{P}(\triangle^{-1}\bm{X}_{1}\cdot\nabla\triangle^{-1}\bm{X}_{2}+\triangle^{-1}\bm{X}_{2}\cdot\nabla\triangle^{-1}\bm{X}_{1}).$

Likewise, at $O(\epsilon^{n+1})$ for $\forall n\in\mathbb{N},$ we have

$\triangle^{*}\bm{X}_{n+1}=-\frac{1}{\nu}\triangle\mathbb{P}\sum_{l=1}^{n}(\triangle^{-1}\bm{X}_{l}\cdot\nabla\triangle^{-1}\bm{X}_{n+1-l}),$

or

$\bm{X}_{n+1}=-(\nu\triangle^{*})^{-1}\triangle\mathbb{P}\sum_{l=1}^{n}(\triangle^{-1}\bm{X}_{l}\cdot\nabla\triangle^{-1}\bm{X}_{n+1-l}).$
(16)

Here use has been made of the inverse Fokker-Planck operator is given by [13]

$(\nu\triangle^{*})^{-1}\equiv-\int_{0}^{\infty}dse^{\nu s\triangle^{*}}=\int
d\bm{\eta}g(\bm{\xi},\bm{\eta})$

with

$g(\bm{\xi},\bm{\eta})\equiv\frac{-1}{(2\pi\nu)^{3/2}}{\rm
f.p.}\int_{\sqrt{a}}^{\infty}\frac{\sigma^{2}d\sigma}{\sigma^{2}-a}e^{-\frac{1}{2\nu}|\sigma\bm{\xi}-\bm{\eta}\sqrt{\sigma^{2}-a}|^{2}},$

where f.p. denotes Hadamard’s finite part. Noting that
$(\triangle^{*})^{-1}\triangle\mathbb{P}$ represents an integral operator of
order zero, we see that the right-hand side of (16) is a functional of
$\bm{X}_{k}$’s $(1\leq k\leq n).$ Because they are all represented as
functionals of $\widehat{\bm{\Psi}}(\bm{\xi})$ and its derivatives, so is
$\bm{X}_{n+1}$. Hence we deduce that $\bm{X}_{n}$ can be expressed as a
functional of $\widehat{\bm{\Psi}}(\bm{\xi})$ and its derivatives for all
$n\in\mathbb{N}.$ $\;\;\square$

We note that that the nonlinear contribution $\bm{X}_{n}\;(n\geq 2)$ has the
second order derivative of $\bm{\Psi}_{1}(\bm{\xi})$ at most and the
$(5-2n)$-th order derivative at least. At $n=3$ already we have $5-2n<0$ and
this means that it may involve integrals of $\bm{\Psi}_{1}(\bm{\xi})$, hence
$F$ may be a functional, rather than a function. We also note that when the
$\epsilon$-expansion is uniformly convergent, $\bm{X}$ itself is a functional
of $\widehat{\bm{\Psi}}(\bm{\xi})$ and its derivatives. To see how we can make
$\epsilon$ arbitrarily small, take for example, the next-to-leading order
approximation (iteration) [13]

$\bm{X}=\bm{X}_{1}-\frac{1}{\nu}\mathbb{P}\left(\triangle^{-1}\bm{X}_{1}\cdot\nabla\triangle^{-1}\bm{X}_{1}\right)$
$=\mathbb{P}\bm{M}G-\frac{1}{\nu}\mathbb{P}\left(\triangle^{-1}\mathbb{P}\bm{M}G\cdot\nabla\triangle^{-1}\mathbb{P}\bm{M}G\right).$
(17)

Introducing the following variables for non-dimensionalisation
$\widetilde{\bm{X}}=\bm{X}/\nu\,\mbox{and}\,\widetilde{\bm{M}}=\bm{M}/\nu,$ we
find

$\widetilde{\bm{X}}=\mathbb{P}\widetilde{\bm{M}}G-\mathbb{P}\left(\triangle^{-1}\mathbb{P}\widetilde{\bm{M}}G\cdot\nabla\triangle^{-1}\mathbb{P}\widetilde{\bm{M}}G\right),$

where $Re=|\widetilde{\bm{M}}|=|\bm{M}|/\nu$ denotes the Reynolds number.
Identifying $Re=\epsilon,$ this corresponds to the $O(\epsilon^{2})$
approximation. We can argue similarly for the $O(\epsilon^{n})$ approximations
as well for $\forall n\in\mathbb{N}$.

On the basis of Proposition 4.1, we consider a self-similar profile in the
following form666Care should be taken of the notations; it means that $\bm{F}$
is a functional of $\widehat{\bm{\Psi}}$ and its derivatives with
$(\nabla\times)^{3}\bm{\widehat{\Psi}}$ at the highest.

$\bm{X}(\bm{\xi})=\bm{F}\left(\bm{\widehat{\Psi}}(\bm{\xi}),\nabla\times\bm{\widehat{\Psi}}(\bm{\xi}),(\nabla\times)^{2}\bm{\widehat{\Psi}}(\bm{\xi});(\nabla\times)^{3}\bm{\widehat{\Psi}}(\bm{\xi})\right)$
(18)

for some functional $\bm{F},$ where $\bm{\widehat{\Psi}}$ denotes the scaled
heat flow,
$(\nabla\times)^{3}\bm{\widehat{\Psi}}=\mathbb{P}\bm{K}\widehat{G},$ $\bm{K}$
is a function of $\bm{M}=\int\bm{X}(\bm{\xi})d\bm{\xi}$ and
$\widehat{G}=\exp(-\frac{a}{2\nu}|\bm{\xi}|^{2}).$ The functional
$F(\bm{x},\bm{y},\bm{z};\bm{w})$ is a near-identity transformation of the last
argument $\bm{w}$ (the solenoidal Gaussian function) and satisfies

$\bm{F}(\bm{x},\alpha\bm{y},\alpha^{2}\bm{z};\alpha^{3}\bm{w})=\alpha^{3}\bm{F}(\bm{x},\bm{y},\bm{z};\bm{w})\;\mbox{for}\;\forall\alpha>0.$
(19)

The profile (18) is a rather strong assumption, even though we take it
approximately. If $\bm{F}$ is available, at least part of solutions to the
Navier-Stokes equations is reducible to those of the heat equations.
Unfortunately, the precise form of $\bm{F}$ is not known at the moment, but we
know that it is close to the Gaussian and how close it is [13]. Assuming that
$F$ is obtained, the procedure goes as follows.

Step 1. When the self-similar profile $\bm{F}$ is given explicitly, a
particular self-similar solution is obtained as

$\bm{\chi}(\bm{x},t)=\frac{1}{(2at)^{3/2}}\bm{F}\left(\bm{\widehat{\Psi}}\left(\frac{\bm{x}}{\sqrt{2at}}\right),\nabla\times\bm{\widehat{\Psi}}\left(\frac{\bm{x}}{\sqrt{2at}}\right),(\nabla\times)^{2}\bm{\widehat{\Psi}}\left(\frac{\bm{x}}{\sqrt{2at}}\right);(\nabla\times)^{3}\bm{\widehat{\Psi}}\left(\frac{\bm{x}}{\sqrt{2at}}\right)\right).$

Step 2. We seek a more general class of solution on this basis

$\bm{\chi}(\bm{x},t)=\frac{1}{(2at)^{3/2}}\bm{F}\left(\bm{\widehat{\psi}}(\bm{x},t),(2at)^{1/2}\,\nabla\times\bm{\widehat{\psi}}(\bm{x},t),2at\,(\nabla\times)^{2}\bm{\widehat{\psi}}(\bm{x},t);(2at)^{3/2}\,(\nabla\times)^{3}\bm{\widehat{\psi}}(\bm{x},t)\right).$

The first three arguments give rise to a factor of $(2at)^{3/2}$ in total, as
a result of scaling (19) and hence we have

$\bm{\chi}(\bm{x},t)=\bm{F}\left(\bm{\widehat{\psi}}(\bm{x},t),\nabla\times\bm{\widehat{\psi}}(\bm{x},t),(\nabla\times)^{2}\bm{\widehat{\psi}}(\bm{x},t);(\nabla\times)^{3}\bm{\widehat{\psi}}(\bm{x},t)\right).$

It is not known how general a class of functions such a construction can
cover.

To see how cancellations take place, it is helpful to consider the following
example based on the next-to-leading order approximation. On the right-most
side of (17), the first term gives rise to $(2at)^{3/2}$ so does the second
$(2at)^{1/2}\cdot(2at)=(2at)^{3/2}.$

Remark
In three dimensions, under the assumption of
$\bm{\Psi}=\bm{f}(\bm{\widehat{\Psi}}),$ a simpler form of $\bm{F}$ can be
obtained using spatial derivatives rather than curls. Write

$\bm{X}=\widetilde{\bm{F}}\left(\bm{\widehat{\Psi}},D\bm{\widehat{\Psi}},D^{2}\bm{\widehat{\Psi}};D^{3}\bm{\widehat{\Psi}}\right),$

where $D_{i}=\partial_{\xi_{i}}\;\;i=1,2,3$ stands for any spatial
derivatives. After some algebra we find explicitly

$X_{i}=-\epsilon_{ikq}\left\\{\frac{\partial
s_{j}}{\partial\xi_{k}}\frac{\partial s_{m}}{\partial\xi_{l}}\frac{\partial
s_{n}}{\partial\xi_{l}}\frac{\partial^{3}f_{q}}{\partial s_{j}\partial
s_{m}\partial s_{n}}+\left(\frac{\partial s_{j}}{\partial\xi_{k}}\triangle
s_{m}+2\frac{\partial^{2}s_{j}}{\partial\xi_{k}\partial\xi_{l}}\frac{\partial
s_{m}}{\partial\xi_{l}}\right)\frac{\partial^{2}f_{q}}{\partial s_{j}\partial
s_{m}}+\frac{\partial\triangle s_{j}}{\partial\xi_{k}}\frac{\partial
f_{q}}{\partial s_{j}}\right\\}$

where we have put $s_{i}=\widehat{\Psi}_{i}(\bm{\xi})$ for simplicity. Observe
that the above expression reduces to
$\bm{X}=-\nabla\times\triangle\widehat{\bm{\Psi}}.$ for the identity
transformation $\frac{\partial f_{q}}{\partial s_{j}}=\delta_{qj}.$

##  5 Summary and outlook

In this paper the following things are done. As part of forward self-similar
solutions to the Burgers equation we have identified a class of kink-type
solutions, on top of the well-known source-type solutions. They are given by
Kummer’s confluent hypergeometric function $M_{\rm K}$. We also noted ill-
behaved ’blow-up’ profiles by flipping the sign of the parameter $a.$

We discussed the Navier-Stokes equations; in two dimensions we discussed a
self-similar profile which may be regarded as a conjugate to the Burgers
vortex. Some of its asymptotic properties have been studied, whereas
clarification of its significance it may play, e.g. its stability as a
dynamical system, requires further investigation.

We also discussed applications of the self-similar solution in three
dimensions. Some properties of the self-similar profile have been analysed
formally and possible lifting to more general class of solutions, at least
approximately, is suggested. For the final topic, it may be worthwhile to try
computing approximate solutions by numerical methods. This is also left for
future study.

##  Appendix A Riccati equation

Consider the Riccati’s ordinary differential equation

$\frac{dy}{dx}=a(x)y^{2}+b(x)y+c(x),$

where $a(x),b(x)$ and $c(x)$ are given functions. It is known that a
substitution of the following form

$y=-\frac{u^{\prime}(x)}{a(x)u(x)}$

reduces the above to a linear second-order homogeneous equation

$\frac{d^{2}u}{dx^{2}}-\left(\frac{a^{\prime}(x)}{a(x)}+b(x)\right)\frac{du}{dx}+a(x)c(x)u=0.$

##  Appendix B Confluent hypergeometric equation

Kummer’s confluent hypergeometric equation [9]

$z\frac{d^{2}w}{dz^{2}}+\left(\gamma-z\right)\frac{dw}{dz}-\alpha w=0$

has two fundamental solutions $w(z)=M_{\rm K}(\alpha,\gamma,z)$ and $U_{\rm
K}(\alpha,\gamma,z)$. Their asymptotic behaviours are given as follows.

For $\Re(z)>0,$ as $|z|\to\infty$ we have

$\left\\{\begin{array}[]{l}M_{\rm
K}(\alpha,\gamma,z)\approx\dfrac{\Gamma(\gamma)}{\Gamma(\alpha)}e^{z}z^{\alpha-\gamma},(\alpha\neq
1,\gamma\neq-n),\\\ \vskip 5.69046pt\cr U_{\rm K}(\alpha,\gamma,z)\approx
z^{-\alpha},\end{array}\right.$

where $\Gamma(\cdot)$ denotes the gamma function.

Also, as $|z|\to 0$ we have

$\left\\{\begin{array}[]{ll}M_{\rm K}(\alpha,\gamma,z)=1,&(\gamma\neq-n),\\\
\vskip 5.69046pt\cr U_{\rm
K}(\alpha,\gamma,z)=\dfrac{\Gamma(\gamma-1)}{\Gamma(\alpha)}z^{1-\gamma}+O(1),&(\alpha\neq
1,1<\Re(\gamma)<2),\end{array}\right.$

where $n\in\mathbb{N}$.

##  Appendix C Source-kink duality

The PDE, Burgers equation, allows at least in two different kinds of
approximations based on ODEs via particle systems. Those particle pictures are
well-known, but best stated here for motivation.

One is
(1) Propagation of Wiener’s chaos, e.g.[16]

$dX_{i}=\frac{1}{2N}\sum_{j\neq i}b(X_{j}-X_{i})dt+dW_{i},$

where $dW$ denotes Brownian motion, the drift velocity $b=\delta$ for the
Burgers equation. This is related with the source-type solution.

The other one is
(2) Pole decomposition, e.g.[15]

$u(x,t)=-2\nu\sum_{j=1}^{N}\frac{1}{x-z_{j}(t)},$
$\frac{dz_{j}}{dt}=-2\nu\sum_{k=1,k\neq j}^{N}\frac{1}{z_{j}-z_{k}},$

where $z_{j}$ denotes the locations of poles in the complex plane. This is
related with the kink-type solution, represented by Kummer’s $M_{\rm K}$.

A crude explanation why we have two different views is as follows. There are
two fundamental solutions to the Fokker-Planck equation; the Gaussian function
and the Dawson’s integral, which are related by the Hilbert transform to each
other. One of them converges to the Dirac mass $\delta$ and the other to a
Cauchy kernel $1/z,$ respectively in suitable limits.

## References

  * [1] Cannone M and Planchon F. 1996. Self-similar solutions for Navier-Stokes equations in $\mathbb{R}^{3}$. Commun. in P.D.E. 21, 179–193.
  * [2] Lemarié-Rieusset PG. 2002 Recent developments in the Navier-Stokes problem CRC Press.
  * [3] Germain P, Pavlović N and Staffilani G. 2007 Regularity of solutions to the Navier-Stokes equations evolving from small data in BMO-1. Intern. Math. Res. Notices 2007.9, rnm087.
  * [4] Bradshaw Z and Tsai T-P. 2018 Self-similar solutions to the Navier-Stokes equations: A survey of recent results. preprint arXiv:1802.00038.
  * [5] Hopf E. 1950. The partial differential equation $u_{t}+uu_{x}=\mu u_{xx}$. Commun. Pure and Appl. Math. 3, 201–230.
  * [6] Escobedo M and Zuazua E. 1991. Large time behavior for convection-diffusion equations in $R^{N}$. J. Func. Anal. 100, 119–161.
  * [7] Biler P, Karch G and Woyczyński WA. 1999. Asymptotics for multifractal conservation laws. Studia Mathematica 135, 231-252
  * [8] Drazin PG and Johnson RS. 1989. Solitons: an introduction Cambridge: Cambridge University Press.
  * [9] Olver FW, Lozier DW, Boisvert RF and Clark CW. (eds.) 2010 NIST Handbook of Mathematical Functions Cambridge: Cambridge University Press.
  * [10] Nec̆as J, Růz̆ic̆ka M and S̆verák V. 1996 On Leray’s self-similar solutions of the Navier-Stokes equations. Acta Math. 176, 283–294.
  * [11] Tsai T-P. 1998 On Leray’s Self-Similar Solutions of the Navier-Stokes Equations Satisfying Local Energy Estimates Arch. Rat. Mech. Anal. 143, 29–51.
  * [12] Ohkitani K. 2020. Study of the Hopf functional equation for turbulence: Duhamel principle and dynamical scaling. Phys. Rev. E 101, 013104-1–15.
  * [13] Ohkitani K and Vanon R. 2021 Self-similar source-type solutions to the three-dimensional Navier-Stokes equations. preprint arXiv:2107.02952.
  * [14] Jia H and Šverák V. 2014. Local-in-space estimates near initial time for weak solutions of the Navier-Stokes equations and forward self-similar solutions. Invent. Math. 196, 233–265.
  * [15] Aref H. 1983 Integrable, chaotic, and turbulent vortex motion in two-dimensional flows. Ann. Rev. Fluid Mech. 15, 345–389.
  * [16] Sznitman AS. 1991 Topics in propagation of chaos. in Ecole d’été de probabilités de Saint-Flour XIX-1989, 165–251, Berlin, Heidelberg: Springer.


# Consistent Discretization of a Class of Predefined-Time Stable Systems

Esteban Jiménez-Rodríguez  Rodrigo Aldana-López  Juan D. Sánchez-Torres  David
Gómez-Gutiérrez  Alexander G. Loukianov  Department of Electrical Engineering,
Cinvestav-Guadalajara, Jalisco, 45019 México (e-mail: {ejimenezr,
<EMAIL_ADDRESS> Department of Computer Science and Systems
Engineering, University of Zaragoza, Zaragoza, 50009 España, (e-mail:
<EMAIL_ADDRESS> Research Laboratory on Optimal Design,
Devices and Advanced Materials -OPTIMA-, Department of Mathematics and
Physics, ITESO, Jalisco, 45604 México (e-mail<EMAIL_ADDRESS> Multi-agent
autonomous systems lab, Intel Labs, Intel Tecnología de México, Jalisco, 45019
Mexico (e-mail<EMAIL_ADDRESS> Tecnologico de Monterrey, Escuela de
Ingeniería y Ciencias, Jalisco, 45138 Mexico

###### Abstract

As the main contribution, this document provides a consistent discretization
of a class of fixed-time stable systems, namely predefined-time stable
systems. In the unperturbed case, the proposed approach allows obtaining not
only a consistent but exact discretization of the considered class of
predefined-time stable systems, whereas in the perturbed case, the consistent
discretization preserves the predefined-time stability property. All the
results are validated through simulations and compared with the conventional
explicit Euler scheme, highlighting the advantages of this proposal.

###### keywords:

Predefined-time stability, Discrete-time systems, Digital implementation,
Stability of nonlinear systems, Fixed-time stability

††thanks: This work has been submitted to IFAC for possible publication

##  1 Introduction

Given a system with tunable parameters whose origin is fixed-time stable,
generally, it is not straightforward, and sometimes it is even impossible to
achieve any desired upper bound of the settling time through the selection of
the parameters of the system (Jiménez-Rodríguez et al., 2019, Example 1). To
overcome this drawback, a class of dynamical systems that exhibit the property
of predefined-time stability has been studied within the last six years
(Sánchez-Torres et al., 2018).

Predefined-time stability refers to the property that exhibits a particular
class of fixed-time stable systems whose solutions converge to the origin
before an arbitrary user-prescribed time, which is assigned through an
appropriate selection of the tunable parameters of the system (Sánchez-Torres
et al., 2018). Due to its remarkable features, the study of the properties and
applications of the predefined-time stability notion has attracted much
attention, mainly in continuous time (Aldana-López et al., 2019, 2019; Muñoz-
Vázquez et al., 2019; Sánchez-Torres et al., 2019).

The numerical simulation examples conducted in the above works were done
applying the conventional explicit (forward) Euler discretization, with a tiny
step size ($1\times 10^{-4}$ or $1\times 10^{-5}$ time units). However, the
right side of the ordinary differential equations exhibiting predefined-time
stability does not satisfy the Lipschitz condition at the origin. In this
case, the explicit Euler discretization scheme does not guarantee that
predefined-time stability property will be preserved (Levant, 2013; Huber et
al., 2016; Polyakov et al., 2019). In other words, the discrete-time equation
solutions may be inconsistent with the solutions of the continuous-time one.

In this sense, this work concerns with the discretization of a class of
predefined-time stable systems. The discretization process uses the idea of
(Polyakov et al., 2019, Figure 1.1), but different from the contribution of
the mentioned paper, which achieves consistent discretizations of homogeneous
finite-time and practically fixed-time stable systems, the proposed approach
allows the consistent discretization of the class of predefined-time stable
systems proposed by Aldana-López et al. (2019); Jiménez-Rodríguez et al.
(2019). Hence, this document provides a consistent discretization of a class
of fixed-time stable systems.

The rest of this paper is organized as follows: $\mathcal{K}^{1}$ functions
and the considered class of predefined-time stable systems are introduced in
Section 2. Then, the exact discretization of the considered class of
predefined-time stable systems is developed in Section 3, in two different
ways. Finally, Section 4 considers a consistent discretization of perturbed
systems with discontinuous predefined-time control.

##  2 Preliminaries

###  2.1 Class $\mathcal{K}^{1}$ functions

The following definition was inspired on the comparison class $\mathcal{K}$
functions (Kellett, 2014, Definition 1)

######  Definition 1 ($\mathcal{K}^{1}$ functions)

A scalar continuous function $\kappa:\mathbb{R}_{\geq 0}\to\left[0,1\right)$
is said to belong to class $\mathcal{K}^{1}$, denoted as
$\kappa\in\mathcal{K}^{1}$, if it is strictly increasing, $\kappa(0)=0$ and
$\lim_{r\to\infty}\kappa(r)=1$.

######  Remark 2 (Invertibility of class $\mathcal{K}^{1}$ functions)

Let $\kappa\in\mathcal{K}^{1}$.

  * (i)

$\kappa$ is injective since it is strictly increasing;

  * (ii)

$\kappa$ is onto since it is continuous, $\kappa(0)=0$ and
$\lim_{r\to\infty}\kappa(r)=1$.

Hence, $\kappa$ is bijective and its inverse exist.

######  Remark 3 (A connection with probability theory)

Let $\kappa\in\mathcal{K}^{1}$ be continuously differentiable on its domain.
In this case, there exists a continuous function $\Phi:\mathbb{R}_{\geq
0}\to\mathbb{R}_{+}$ such that $\kappa(r)=\int_{0}^{r}\Phi(z)\text{d}z$.
Moreover, since $\lim_{r\to\infty}\kappa(r)=1$, the function $\Phi$ is
required to satisfy $\int_{0}^{\infty}\Phi(z)\text{d}z=1$, i.e., functions
$\Phi$ and $\kappa$ can be viewed as probability density functions and
cumulative distribution functions, respectively, of positive random variables.

###### Example 4

Examples of class $\mathcal{K}^{1}$ functions are:

  * •

$\kappa(r)=\frac{2}{\pi}\arctan(ar)$, with $a>0$;

  * •

$\kappa(r)=\frac{r}{r+a}$, with $a>0$;

  * •

$\kappa(r)=1-a^{-r}$, with $a>1$;

  * •

$\kappa(r)=P(a,r)$, with $a>0$, where $P(a,r)=\frac{\gamma(a,r)}{\Gamma(a)}$
is the regularized Incomplete Gamma Function,
$\gamma(a,r)=\int_{0}^{r}t^{a-1}e^{-t}\text{d}t$ is the Incomplete Gamma
Function and $\Gamma(a)=\int_{0}^{\infty}t^{a-1}e^{-t}\text{d}t$ is the Gamma
Function (Abramowitz and Stegun, 1965);

  * •

$\kappa(r)=I(a_{1},a_{2},\frac{r}{r+1})$, with $a_{1},a_{2}>0$, where
$I(a_{1},a_{2},r)=\frac{b(a_{1},a_{2},r)}{\mathcal{B}(a_{1},a_{2})}$ is the
regularized Incomplete Beta Function,
$b(a_{1},a_{2},r)=\int_{0}^{r}t^{a_{1}-1}(1-t)^{a_{2}-1}\text{d}t$ is the
Incomplete Beta Function and
$\mathcal{B}(a_{1},a_{2})=\int_{0}^{1}t^{a_{1}-1}(1-t)^{a_{2}-1}\text{d}t$ is
the Beta Function (Abramowitz and Stegun, 1965).

###  2.2 A class of predefined-time stable systems

Consider (the scalar form of) the class of systems introduced by Aldana-López
et al. (2019) and Jiménez-Rodríguez et al. (2019)

$\dot{x}=-\frac{1}{\rho_{1}(1-\rho_{2})}\frac{\kappa(\left\lvert
x\right\rvert)^{\rho_{2}}}{\kappa^{\prime}(\left\lvert
x\right\rvert)}\text{sign}\left(x\right),\ \ x(0)=x^{0},$ (1)

where $x:\mathbb{R}_{\geq 0}\to\mathbb{R}$ is the state of the system,
$\rho_{1}>0$ and $0\leq\rho_{2}<1$ are tunable parameters,
$\kappa\in\mathcal{K}^{1}$ is continuously differentiable, and
$\text{sign}\left(\cdot\right)$ stands for the signum function given by

$\text{sign}\left(x\right)\in\left\\{\begin{array}[]{ccc}\\{1\\}&\text{if}&x>0\\\
\\{-1\\}&\text{if}&x<0\\\ {[-1,1]}&\text{if}&x=0.\end{array}\right.$

From the theory of ordinary differential equations, (1) is a separable first-
order equation, whose solution $x(t)$ satisfies

$\int_{x^{0}}^{x(t)}(1-\rho_{2})\frac{\kappa^{\prime}(\left\lvert\xi\right\rvert)}{\kappa(\left\lvert\xi\right\rvert)^{\rho_{2}}}\text{sign}\left(\xi\right)\mathrm{d}\xi=-\int_{0}^{t}\frac{1}{\rho_{1}}\mathrm{d}\tau.$

Applying the change of variable
$\xi\to\kappa(\left\lvert\xi\right\rvert)^{1-\rho_{2}}\text{sign}\left(\xi\right)$
in the left side integral and integrating both sides of the above equation, it
yields

$\kappa(\left\lvert x(t)\right\rvert)^{1-\rho_{2}}=\kappa(\left\lvert
x^{0}\right\rvert)^{1-\rho_{2}}-\frac{t}{\rho_{1}},$

for $\kappa(\left\lvert x^{0}\right\rvert)^{1-\rho_{2}}-\frac{t}{\rho_{1}}\geq
0$. From the above, the solution of (1) is

$x(t)=\\\
\small{\left\\{\begin{array}[]{l}\kappa^{-1}\left(\left[\kappa\left(\left\lvert
x^{0}\right\rvert\right)^{1-\rho_{2}}-\frac{t}{\rho_{1}}\right]^{\frac{1}{1-\rho_{2}}}\right)\text{sign}\left(x^{0}\right)\\\
\qquad\qquad\qquad\qquad\qquad\qquad\qquad\text{ if }0\leq
t\leq\rho_{1}\kappa\left(\left\lvert x^{0}\right\rvert\right)^{1-\rho_{2}}\\\
0\\\ \qquad\qquad\qquad\qquad\qquad\qquad\qquad\text{ if
}t>\rho_{1}\kappa\left(\left\lvert
x^{0}\right\rvert\right)^{1-\rho_{2}},\end{array}\right.}$ (2)

where $\kappa^{-1}$ stands for the inverse of the function
$\kappa\in\mathcal{K}$.

Hence, one can conclude that the origin $x=0$ of (1) is finite-time stable and
the settling-time function is $T(x^{0})=\rho_{1}\kappa\left(\left\lvert
x^{0}\right\rvert\right)^{1-\rho_{2}}$. But the settling-time function
$T(x^{0})=\rho_{1}\kappa\left(\left\lvert
x^{0}\right\rvert\right)^{1-\rho_{2}}<\rho_{1}$ is bounded, then the origin
$x=0$ of (1) is, in fact, fixed-time stable.

Moreover, for any desired upper bound of the settling-time function $T_{c}$,
there exists an appropriate selection of the parameter $\rho_{1}=T_{c}$ such
that the convergence time will always be less than $T_{c}$ for any initial
condition $x^{0}$. This property is known as predefined-time stability (see
Appendix A).

###### Remark 5

It is worth noticing that the right-hand side of (1) is continuous and non-
Lipschitz if $0<\rho_{2}<1$, and discontinuous if $\rho_{2}=0$. In any case,
all the above derivation remains true.

##  3 Exact discretization of the considered class of predefined-time stable
systems

In this section, we develop a discrete-time system whose solutions are equal
to the solution of (1) in the sampling instants. Hence, we provide the exact
form of simulating the class of predefined-time stable systems represented by
(1).

###  3.1 Direct construction

Consider the continuous-time system (1) and let $x_{k}=x(kh)$, where $h>0$ is
the sampling step size, and $k\in\mathbb{N}$.

To obtain the exact discretization of system (1), one can apply the same
procedure as in Subsection 2.2, integrating the left side between the samples
$x_{k}$ and $x_{k+1}$, and the right side between $kh$ and $(k+1)h$,
respectively, obtaining

$\int_{x_{k}}^{x_{k+1}}(1-\rho_{2})\frac{\kappa^{\prime}(\left\lvert\xi\right\rvert)}{\kappa(\left\lvert\xi\right\rvert)^{\rho_{2}}}\text{sign}\left(\xi\right)\mathrm{d}\xi=-\int_{kh}^{(k+1)h}\frac{1}{\rho_{1}}\mathrm{d}\tau.$
(3)

Recall the transformation

$y_{\rho_{2}}\colon\begin{array}[]{>{\displaystyle}r @{} >{{}}c<{{}} @{}
>{\displaystyle}l}\mathbb{R}&\rightarrow&(-1,1)\\\
x&\mapsto&y_{\rho_{2}}(x)=\kappa(\left\lvert
x\right\rvert)^{1-\rho_{2}}\text{sign}\left(x\right)\end{array}$ (4)

used in Subsection 2.2 to find the solutions of (1). This transformation is
invertible, and its inverse $y_{\rho_{2}}^{-1}:(-1,1)\to\mathbb{R}$ is given
by

$y_{\rho_{2}}^{-1}(x)=\kappa^{-1}\left(\left\lvert
x\right\rvert^{\frac{1}{1-\rho_{2}}}\right)\text{sign}\left(x\right).$

Using (4) in the left side integral of (3), one obtains

$\int_{\kappa(\left\lvert
x_{k}\right\rvert)^{1-\rho_{2}}\text{sign}\left(x_{k}\right)}^{\kappa(\left\lvert
x_{k+1}\right\rvert)^{1-\rho_{2}}\text{sign}\left(x_{k+1}\right)}\text{sign}\left(y\right)\mathrm{d}y=-\int_{kh}^{(k+1)h}\frac{1}{\rho_{1}}\mathrm{d}\tau.$

Hence, integrating both sides of the above expression and having in mind that
$\text{sign}\left(x_{k+1}\right)=\text{sign}\left(x_{k}\right)$ for the first-
order system (1), the above equation turns into

$\kappa(\left\lvert
x_{k+1}\right\rvert)^{1-\rho_{2}}\text{sign}\left(x_{k+1}\right)\\\
=\kappa(\left\lvert
x_{k}\right\rvert)^{1-\rho_{2}}\text{sign}\left(x_{k}\right)-\frac{h}{\rho_{1}}\text{sign}\left(x_{k+1}\right).$
(5)

Finally, let $\alpha\in\mathbb{R}_{+}$ and consider the discontinuous function
$\theta_{\rho_{2},\alpha}:\mathbb{R}\to(-\alpha-1,\alpha+1)$ defined by
$\theta_{\rho_{2},\alpha}(x)=\left(\kappa(\left\lvert
x\right\rvert)^{1-\rho_{2}}+\alpha\right)\text{sign}\left(x\right)$. The
inverse of $\theta_{\rho_{2},\alpha}$ is the continuous function
$\theta_{\rho_{2},\alpha}^{-1}:(-\alpha-1,\alpha+1)\to\mathbb{R}$ given by

$\theta_{\rho_{2},\alpha}^{-1}(x)=\kappa^{-1}\left(\left(\max\\{\left\lvert
x\right\rvert-\alpha,0\\}\right)^{\frac{1}{1-\rho_{2}}}\right)\text{sign}\left(x\right)$
(6)

Thus, replacing (6) in (5), it yields

$\displaystyle\begin{split}x_{k+1}&=\theta_{\rho_{2},h/\rho_{1}}^{-1}\left(\kappa(\left\lvert
x_{k}\right\rvert)^{1-\rho_{2}}\text{sign}\left(x_{k}\right)\right),\\\
&=\kappa^{-1}\left(\max\left\\{\kappa(\left\lvert
x_{k}\right\rvert)^{1-\rho_{2}}-\frac{h}{\rho_{1}},0\right\\}^{\frac{1}{1-\rho_{2}}}\right)\text{sign}\left(x_{k}\right),\end{split}$
(7)

with $x_{0}=x^{0}$.

This following theorem summarizes the previous constructive result:

###### Theorem 6

The solution of the discrete-time system (7) is equivalent to the solution (2)
of the continuous-time system (1) at the sample instants. This is,
$x_{k}=x(kh)$, for every $k\in\mathbb{N}$.

Theorem 6 has the following attractive corollary:

###### Corollary 7

For any $x^{0}\in\mathbb{R}$, the solution of the discrete-time system (7)
satisfies $x_{k}=0,\text{ for
}k\geq\left\lceil\frac{\rho_{1}}{h}\right\rceil,$ where
$\left\lceil\cdot\right\rceil$ stands for the ceiling function.

The solution of the continuous-time system (1) satisfies

$x(t)=0,\text{ for }t\geq\rho_{1}\kappa\left(\left\lvert
x^{0}\right\rvert\right)^{1-\rho_{2}}.$

Moreover, from Theorem 6, we can establish for the solution of the discrete-
time system (7) that

$x_{k}=x(kh)=0,\text{ for }kh\geq\rho_{1}\kappa\left(\left\lvert
x^{0}\right\rvert\right)^{1-\rho_{2}}.$

which also holds for $kh\geq\rho_{1}$, given that $\kappa\in\mathcal{K}$.

Hence, the result follows. $\blacksquare$

###### Example 8

Consider system (1) with the particular selections of
$\kappa(r)=\frac{2}{\pi}\arctan(r)$, $\rho_{1}=1$, $\rho_{2}=0.5$ and
$x^{0}=10$.

Figures 1 and 2 present the comparison of the exact discretization (7) with
the conventional explicit Euler discretization for two different time steps.
In both cases, the solution of the discrete-time system (7) accurately
recovers the exact solution (2) at the sampling instants. For $h=0.02$ (Fig.
1), one can note small differences between the exact solution and the explicit
Euler discretization, whereas, for $h=0.0651318636497$ (Fig. 2), the explicit
Euler discretization becomes unstable.

Figure 1: Comparison of the exact solution (solid line), the discrete-time
equivalent (star points) and the explicit Euler discretization (round points)
with a step size of $h=0.02$. Figure 2: Comparison of the exact solution
(solid line), the discrete-time equivalent (star points) and the explicit
Euler discretization (round points) with a step size of $h=0.0651318636497$.

###  3.2 Parallel with implicit Euler discretization

Consider the nonlinear transformation given by (4), and define
$w(t)=y_{\rho_{2}}(x(t))$. Its time-derivative along the trajectories of
system (1) is

$\dot{w}=-\frac{1}{\rho_{1}}\text{sign}\left(w\right),\ \
w(0)=y_{\rho_{2}}(x^{0}),$ (8)

since $\text{sign}\left(w\right)=\text{sign}\left(x\right)$.

###### Remark 9

If the explicit or forward Euler procedure with time step $h>0$ were used to
simulate system (8), it would lead to oscillations with amplitude
$\mathcal{O}(h)$ in the neighborhood of the origin $w=0$ (Drakunov and Utkin,
1989; Utkin, 1994), and these oscillations would also be translated into
oscillations in the original variable $x$. Those mentioned oscillations, known
as numerical chattering, are an undesired effect of an inadequate
discretization, since the solution (2) does not oscillate.

Using the implicit or backward Euler method (as in (Drakunov and Utkin, 1989;
Utkin, 1994)) with time step $h>0$ for the discretization of the transformed
system (8), it yields

$w_{k+1}=w_{k}-\frac{h}{\rho_{1}}\text{sign}\left(w_{k+1}\right),$

where $w_{k}=y_{\rho_{2}}(x_{k})$ and $x_{k}=x(kh)$. Rewriting the above in
terms of $x_{k}$, one obtains (3), as before.

##  4 Consistent discretization of predefined-time control of perturbed
systems

###  4.1 Continuous-time predefined-time stabilization of first order
perturbed systems

Consider the perturbed control system

$\dot{x}=u+\Delta(t,x),\,\,x(0)=x^{0},$ (9)

where $x:\mathbb{R}_{\geq 0}\to\mathbb{R}$ is the state of the system,
$u\in\mathbb{R}$ is the control input signal, and $\Delta:\mathbb{R}_{\geq
0}\times\mathbb{R}\to\mathbb{R}$ is an unknown perturbation term which is
assumed to be bounded of the form $\sup_{(t,x)\in\mathbb{R}_{\geq
0}\times\mathbb{R}}\left\lvert\Delta(t,x)\right\rvert\leq\delta$, with
$\delta\in\mathbb{R}_{+}$ known.

It is well known that, in continuous time, the perturbation term $\Delta(t,x)$
can only be entirely rejected by a discontinuous control term, since no
conditions of smoothness, Lipschitz continuity neither continuity are assumed
(Jiménez-Rodríguez et al., 2019).

According to the above, a suitable feedback predefined-time controller design
is

$u=-\left(\frac{1}{\rho_{1}}+\rho_{3}\kappa^{\prime}(0)\right)\frac{1}{\kappa^{\prime}(\left\lvert
x\right\rvert)}\text{sign}\left(x\right),$ (10)

where $\rho_{1}>0$, $\rho_{3}\geq\delta$, and $\kappa\in\mathcal{K}^{1}$ is
continuously differentiable and such that
$\kappa^{\prime}(0):=\kappa^{\prime}(0^{+})\geq\kappa^{\prime}(r)$ for all
$r\in\mathbb{R}_{\geq 0}$. Note that (10) is an augmented gain version of the
right side of (1) with $\rho_{2}=0$.

To verify the predefined-time convergence of the closed-loop system (9)-(10),
note that

$\displaystyle\frac{d\left\lvert x\right\rvert}{dt}$
$\displaystyle=-\left(\frac{1}{\rho_{1}}+\rho_{3}\kappa^{\prime}(0)\right)\frac{1}{\kappa^{\prime}(\left\lvert
x\right\rvert)}+\Delta(t,x)\text{sign}\left(x\right)$
$\displaystyle\leq-\frac{1}{\rho_{1}\kappa^{\prime}(\left\lvert
x\right\rvert)}-\rho_{3}\frac{\kappa^{\prime}(0)}{\kappa^{\prime}(\left\lvert
x\right\rvert)}+\left\lvert\Delta(t,x)\right\rvert$
$\displaystyle\leq-\frac{1}{\rho_{1}\kappa^{\prime}(\left\lvert
x\right\rvert)}-(\rho_{3}-\delta)$
$\displaystyle\leq-\frac{1}{\rho_{1}\kappa^{\prime}(\left\lvert
x\right\rvert)},$

for $x\neq 0$. Reasoning as in Subsection 2.2, and using the comparison lemma
in the above differential inequality, one can easily conclude that
$\left\lvert x(t)\right\rvert=x(t)=0$ for $t\geq\rho_{1}$.

###  4.2 Consistent discretization/simulation

Consider transformation (4), with $\rho_{2}=0$, given by
$z(t)=y_{0}(x(t))=\kappa(\left\lvert
x(t)\right\rvert)\text{sign}\left(x(t)\right)$. Its derivative along the
trajectories of the closed-loop system (9)-(10) is

$\dot{z}=-\beta\text{sign}\left(z\right)+f(t,x),\,\,z(0)=y_{0}(x^{0}),$ (11)

where $\beta:=\frac{1}{\rho_{1}}+\rho_{3}\kappa^{\prime}(0)$, and
$f(t,x)=\kappa^{\prime}(\left\lvert x\right\rvert)\Delta(t,x)$ is the
perturbation term in the transformed coordinate $z$, which complies to
$\sup_{(t,x)\in\mathbb{R}_{\geq 0}\times\mathbb{R}}\left\lvert
f(t,x)\right\rvert\leq\kappa^{\prime}(0)\delta$.

Again, using the implicit Euler method with time step $h>0$ for the
discretization of the system (11), it yields

$z_{k+1}=z_{k}-h\beta\text{sign}\left(z_{k+1}\right)+hf_{k},$ (12)

where $f_{k}=f(kh,x_{k})$ is the perturbation term at the sample instants.

Hence, replacing the function $\theta_{0,h\beta}^{-1}$ (6) into (12), one
obtains

$x_{k+1}=\\\ \kappa^{-1}\left(\max\left\\{\left\lvert\kappa(\left\lvert
x_{k}\right\rvert)\text{sign}\left(x_{k}\right)+hf_{k}\right\rvert-h\beta,0\right\\}\right)\times\\\
\text{sign}\left(\kappa(\left\lvert
x_{k}\right\rvert)\text{sign}\left(x_{k}\right)+hf_{k}\right),$ (13)

with $x_{0}=x^{0}$.

The discrete-time system (13) provides a way of simulating the closed-loop
continuous-time system (9)-(10). Moreover, it preserves the predefined-time
stability behavior, as stated in the following proposition:

###### Proposition 10

For any $x^{0}\in\mathbb{R}$, the solution of the discrete-time system (13)
satisfies $x_{k}=0,\text{ for
}k\geq\left\lceil\frac{\rho_{1}}{h}\right\rceil$.

Taking absolute value in both sides of (13) and using the triangle inequality
one gets,

$\displaystyle\left\lvert x_{k+1}\right\rvert$
$\displaystyle=\kappa^{-1}\left(\max\left\\{\left\lvert\kappa(\left\lvert
x_{k}\right\rvert)\text{sign}\left(x_{k}\right)+hf_{k}\right\rvert-h\beta,0\right\\}\right)$
$\displaystyle\leq\kappa^{-1}\left(\max\left\\{\kappa(\left\lvert
x_{k}\right\rvert)+h\left\lvert f_{k}\right\rvert-h\beta,0\right\\}\right).$

On the other hand, from
$\beta:=\frac{1}{\rho_{1}}+\rho_{3}\kappa^{\prime}(0)$, $\rho_{3}\geq\delta$,
$\left\lvert f_{k}\right\rvert\leq\kappa^{\prime}(0)\delta$ and the above, the
inequality

$\left\lvert
x_{k+1}\right\rvert\leq\kappa^{-1}\left(\max\left\\{\kappa(\left\lvert
x_{k}\right\rvert)-\frac{h}{\rho_{1}},0\right\\}\right)$

is obtained.

Finally, noticing that the right side of the above inequality is equal to the
absolute value of the right side of (7) with $\rho_{2}=0$, and that it is a
non-decreasing function of $x_{k}$, we get

$\left\lvert x_{k}\right\rvert\leq 0,\,\,\forall
k\geq\left\lceil\frac{\rho_{1}}{h}\right\rceil$

applying Lemma 13. $\blacksquare$

###### Example 11

Consider the closed-loop system (9)-(10) with the particular selections of
$\kappa(r)=\frac{2}{\pi}\arctan(r)$, $\rho_{1}=1$, $\rho_{3}=1.1$ and
$x^{0}=10$. It is assumed that $\Delta(t,x)=\sin{10\pi t}$.

Figure 3 present the comparison of the consistent discretization (13) with the
conventional explicit Euler discretization for $h=0.05$. One can see that
$x_{k}=0$ for $k\geq\left\lceil\frac{\rho_{1}}{h}\right\rceil$ for the
consistent discretization (13), whereas the explicit Euler discretization
induce undesired oscillations around the origin.

Figure 3: Comparison of the the discrete-time equivalent (star points) and the
explicit Euler discretization (round points) with a step size of $h=0.05$ for
the predefined-time control of perturbed systems.

##  5 Conclusion

This paper presented the development of a consistent discretization for the
class of predefined-time stable systems introduced in Aldana-López et al.
(2019); Jiménez-Rodríguez et al. (2019).

The proposed approach allowed the exact discretization of the considered class
of systems when no perturbations were assumed. In the perturbed case, the
developed consistent discretization preserved the predefined-time stability
property.

All the results were confirmed through numerical simulations and compared with
the conventional explicit Euler scheme. Even with relatively large time steps,
the proposed discretization worked well, as expected, whereas the explicit
Euler discretization produced unstable oscillations.

Esteban Jiménez acknowledges to CONACYT–México for the D.Sc. scholarship
number 481467 and the project 252405.

## References

  * Abramowitz and Stegun (1965) Abramowitz, M. and Stegun, I.A. (1965).  _Handbook of Mathematical Functions: With Formulas, Graphs, and Mathematical Tables_.  Dover Publications, 9 edition.
  * Aldana-López et al. (2019) Aldana-López, R., Gómez-Gutiérrez, D., Jiménez-Rodríguez, E., Sánchez-Torres, J.D., and Defoort, M. (2019).  Enhancing the settling time estimation of a class of fixed-time stable systems.  _International Journal of Robust and Nonlinear Control_ , 29(12), 4135–4148.
  * Aldana-López et al. (2019) Aldana-López, R., Gómez-Gutiérrez, D., Jiménez-Rodríguez, E., Sánchez-Torres, J.D., and Defoort, M. (2019).  On the design of new classes of fixed-time stable systems with predefined upper bound for the settling time.
  * Bitsoris and Gravalou (1995) Bitsoris, G. and Gravalou, E. (1995).  Comparison principle, positive invariance and constrained regulation of nonlinear systems.  _Automatica_ , 31(2), 217 – 222.  https://doi.org/10.1016/0005-1098(94)E0044-I.
  * Drakunov and Utkin (1989) Drakunov, S. and Utkin, V. (1989).  On discrete-time sliding modes.  _IFAC Proceedings Volumes_ , 22(3), 273 – 278.  https://doi.org/10.1016/S1474-6670(17)53647-2.
  * Filippov (1988) Filippov, A.F. (1988).  _Differential equations with discontinuous righthand sides_.  Kluwer Academic Publishers Group, Dordrecht.
  * Huber et al. (2016) Huber, O., Acary, V., and Brogliato, B. (2016).  Lyapunov stability and performance analysis of the implicit discrete sliding mode control.  _IEEE Transactions on Automatic Control_ , 61(10), 3016–3030.  10.1109/TAC.2015.2506991.
  * Jiménez-Rodríguez et al. (2019) Jiménez-Rodríguez, E., Muñoz Vázquez, A.J., Sánchez-Torres, J.D., Defoort, M., and Loukianov, A.G. (2019).  A Lyapunov-like Characterization of Predefined-Time Stability.  _ArXiv e-prints_.  URL https://arxiv.org/abs/1910.14604.
  * Kellett (2014) Kellett, C.M. (2014).  A compendium of comparison function results.  _Mathematics of Control, Signals, and Systems_ , 26(3), 339–374.  10.1007/s00498-014-0128-8.
  * Levant (2013) Levant, A. (2013).  On fixed and finite time stability in sliding mode control.  In _52nd IEEE Conference on Decision and Control_ , 4260–4265.
  * Muñoz-Vázquez et al. (2019) Muñoz-Vázquez, A.J., Sánchez-Torres, J.D., Jiménez-Rodríguez, E., and Loukianov, A. (2019).  Predefined-time robust stabilization of robotic manipulators.  _IEEE/ASME Transactions on Mechatronics_ , 1–1.  10.1109/TMECH.2019.2906289.
  * Polyakov et al. (2019) Polyakov, A., Efimov, D., and Brogliato, B. (2019).  Consistent discretization of finite-time and fixed-time stable systems.  _SIAM Journal on Control and Optimization_ , 57(1), 78–103.  10.1137/18M1197345.
  * Sánchez-Torres et al. (2019) Sánchez-Torres, J.D., Defoort, M., and Muñoz-Vázquez, A.J. (2019).  Predefined-time stabilisation of a class of nonholonomic systems.  _International Journal of Control_ , 1–8.  10.1080/00207179.2019.1569262.
  * Sánchez-Torres et al. (2018) Sánchez-Torres, J.D., Gómez-Gutiérrez, D., López, E., and Loukianov, A.G. (2018).  A class of predefined-time stable dynamical systems.  _IMA Journal of Mathematical Control and Information_ , 35(Suppl 1), i1–i29.  10.1093/imamci/dnx004.
  * Utkin (1994) Utkin, V.I. (1994).  Sliding mode control in discrete-time and difference systems.  In A. Zinober (ed.), _Variable Structure and Lyapunov Control_ , chapter 5, 87–107. Springer, Berlin Heidelberg.

##  Appendix A Predefined-time stability

Predefined-time stability refers to the property that exhibits a particular
class of fixed-time stable systems with tunable parameters, for which an upper
bound of the settling-time function can be arbitrarily chosen through a
suitable selection of the parameters (Jiménez-Rodríguez et al., 2019,
Defintion 2).

This notion is formally defined considering an autonomous system of the form

$\dot{\bm{x}}=\bm{f}(\bm{x};\bm{\rho}),\ \ \bm{x}(0)=\bm{x}_{0},$ (14)

with $\bm{x}:\mathbb{R}_{\geq 0}\to\mathbb{R}^{n}$ the system state,
$\bm{\rho}\in\mathbb{R}^{l}$ the tunable parameters of (14), and
$\bm{f}:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ a nonlinear function.

For the origin $\bm{x}=\bm{0}$ to be a finite-, in particular, a predefined-
time stable equilibrium of (14), the function $\bm{f}$ must be a non Lipschitz
(maybe discontinuous) function of $\bm{x}$. Then, $\bm{f}$ is assumed to be
such that the solutions of (14) exist and are unique in the sense of Filippov
(Filippov, 1988).

###### Definition 12

The origin of (14) is said to be predefined-time stable if it is fixed-time
stable and for any $T_{c}\in\mathbb{R}_{+}$, there exists some
$\bm{\rho}\in\mathbb{R}^{l}$ such that the settling-time function of (14)
satisfies

$\sup_{\bm{x}_{0}\in\mathbb{R}^{n}}T(\bm{x}_{0})\leq T_{c}.$

##  Appendix B A comparison lemma for discrete-time systems

The following is a particular case of (Bitsoris and Gravalou, 1995,
Proposition 1).

###### Lemma 13

Let $u_{k}$ be the solution of

$u_{k+1}=f(u_{k}),\,\,u_{0}=u^{0},\,\,u\in\mathbb{R},$

where $f:\mathbb{R}\to\mathbb{R}$ is continuous and non-decreasing, and let
$v_{k}$ be such that

$v_{k+1}\leq f(v_{k}),\,\,v_{0}\leq u^{0},\,\,v\in\mathbb{R}.$

Then, $v_{k}\leq u_{k}$, for all $k\in\mathbb{N}$.

We proceed by induction. The base case follows from the hypothesis

$v_{0}\leq u^{0}=u_{0}.$

Now, assume that $v_{k}\leq u_{k}$. Then,

$v_{k+1}\leq f(v_{k})\leq f(u_{k})=u_{k+1},$

since $f$ is non-decreasing. The result follows. $\blacksquare$