# Self-regulated biological transportation structures with general entropy
dissipations, part I: the 1D case

Clarissa Astuto  King Abdullah University of Science and Technology (KAUST),
4700, Thuwal, Saudi Arabia  Jan Haskovec  King Abdullah University of Science
and Technology (KAUST), 4700, Thuwal, Saudi Arabia  Peter Markowich  King
Abdullah University of Science and Technology (KAUST), 4700, Thuwal, Saudi
Arabia  Department of Mathematics, University of Vienna, Oskar-Morgenstern-
Platz 1, 1090, Austria  Simone Portaro  King Abdullah University of Science
and Technology (KAUST), 4700, Thuwal, Saudi Arabia

###### Abstract

We study self-regulating processes modeling biological transportation networks
as presented in [15]. In particular, we focus on the 1D setting for Dirichlet
and Neumann boundary conditions. We prove an existence and uniqueness result
under the assumption of positivity of the diffusivity $D$. We explore
systematically various scenarios and gain insights into the behavior of $D$
and its impact on the studied system. This involves analyzing the system with
a signed measure distribution of sources and sinks. Finally, we perform
several numerical tests in which the solution $D$ touches zero, confirming the
previous hints of local existence in particular cases.

##  1 Introduction

In recent years, there has been growing interest in studying principles and
mechanisms of formation of biological transportation structures. This line of
research has numerous applications beyond biology, in particular in fields
such as medicine, chemistry and engineering. This paper contributes to the
discussion by studying the spatially one–dimensional version of self-regulated
processes modeling biological transportation structures. Such structures, in
particular organization of leaf venation networks, vascular and neural
networks, have been in the focus of biophysical and mathematical research
communities, see, e.g., the recent works [2, 3, 4, 7, 8, 12, 16]. One is
typically interested in geometric, topological and statistical properties of
the optimal transportation structures, where the optimality criteria typically
involve a trade-off between cost and transportation efficiency [6]. Biological
transportation structures develop without centralized control [19] and
therefore they can be considered as emergent structures resulting from self-
regulating processes.

In the recent work [15] a class of self-regulating processes has been
introduced, which is governed by minimization of the entropy dissipation

(1.1)

The functional $E=E[D]$ is to be minimized with respect to the symmetric and
positive semidefinite diffusivity tensor field $D=D(x)$. Here
$\Phi:\mathbb{R}\rightarrow\mathbb{R}$ is the convex entropy generating
function and $u=u(x)$ denotes the concentration of the transported medium,
e.g., chemical species, ions or nutrients. Moreover,
$\Omega\subset\mathbb{R}^{d}$, $d\geq 1$ is a domain with $C^{1}$ boundary.

The functional (1.1) is coupled to the mass conservation equation

$-\mathrm{div}(D\nabla u)=S\quad\mathrm{in}\,\,\Omega,$

where $S=S(x)$ represents the prescribed, time-independent distribution of
sources and sinks of the system. We equip (1) with the Dirichlet boundary
condition $u=c\quad\mathrm{on}\,\,\partial\Omega.$ Here $c\in\mathbb{R}$ is a
constant representing the equilibrium of the system. Indeed, we assume $\Phi$
to have a critical point at $c$, i.e., $\Phi^{\prime}(c)=0.$

The formal $L^{2}$-gradient flow of the energy functional (1.1) constrained by
(1) is of the form $\frac{\partial D}{\partial t}=\Phi^{\prime\prime}(u)\nabla
u\otimes\nabla u+\frac{\nabla\sigma\otimes\nabla u+\nabla
u\otimes\nabla\sigma}{2}\qquad\mbox{in }(0,\infty)\times\Omega,$ see Lemma 1
of [15] for details. Here $t\geq 0$ is the time-like variable induced by the
gradient flow and $\sigma=\sigma(x,t)$ is the solution of the boundary value
problem $-\mathrm{div}(D\nabla\sigma)=\Phi^{\prime\prime\prime}(u)\nabla
u\cdot D\nabla u\qquad\mathrm{in}\,\,(0,\infty)\times\Omega,$ subject to the
homogeneous Dirichlet boundary condition $\sigma=0$ on $\partial\Omega$.

In this paper we focus on the spatially one–dimensional setting of (1–1). In
particular, we shall formulate sufficient conditions for global existence of
solutions. Moreover, we provide insights into global and local well-posedness
propeties of the system by discussing several special cases and constructing
relevant examples. The last part of the paper is dedicated to an appropriate
discretization of (1–1) presentation of several numerical examples.

While investigating a numerical scheme that solves the issue of different
scales in time, a solution has been proposed that treats a specific component
of the system implicitly, while keeping the rest explicit. By adopting this
approach, it becomes possible to create a time discretization wherein the
implicit portion of the system can be inverted relatively easily, often
bypassing the need for solving nonlinear systems. At the same time, the
explicit part maintains its accuracy and capability to capture the stiffness.
The numerical scheme that addresses these issues is the IMEX Runge–Kutta
scheme [9, 17, 18] and it becomes possible to develop high-order linearly
implicit schemes using the conventional Runge–Kutta methods. They are
generally associated to a Butcher’s tableau, that is composed by a matrix of
coefficients $A=\\{a_{i,j}\\}$, a vector of weights $b_{i}$, and a vector of
nodes $c_{i}$.

In this paper we shall employ this methodology to obtain efficient
semi–implicit discretizations that offer high order accuracy in the framework
of biological transportation structures. IMEX schemes are proven to be very
effective for many applications, e.g., for Navier-Stokes equations [10], Euler
equations [5] and for generic high order PDEs [11]. In this paper we make use
of a third–order semi–implicit scheme in time for the spatially
one–dimensional version of (1–1).

##  2 The spatially 1D setting

Without loss of generality, let $\Omega=(0,1)$ and the equilibrium $c=0$. The
one–dimensional version of (1), (1) and (1) reads

$\displaystyle-(Du_{x})_{x}=S,$ (2.1a)
$\displaystyle-(D\sigma_{x})_{x}=\Phi^{\prime\prime\prime}(u)Du_{x}^{2},$
(2.1b) $\displaystyle D_{t}=\Phi^{\prime\prime}(u)u_{x}^{2}+u_{x}\sigma_{x},$
(2.1c)

for all $(x,t)\in(0,1)\times(0,\infty)$, where $u_{x}$ denotes the derivative
of $u$ with respect to the spatial variable. The system is equipped with
homogeneous Dirichlet boundary conditions for $u$ and $\sigma$ and with
uniformly positive and bounded initial condition
$D(0,x)=D_{I}(x)\geq\gamma>0\quad\forall x\in[0,1],\quad D_{I}\in
L^{\infty}(0,1).$ To simplify the notation, we will denote the explicit
dependence of the quantities of interest on either the spatial coordinate $x$
or temporal coordinate $t$ only when the quantity in question is dependent
solely on one of these variables. Otherwise, we will omit it.

System (2.1a–2.1c) can be formulated in an integro-differential form. By
integrating (2.1a) with respect to $x$ we get

$\displaystyle u_{x}=\frac{\alpha(t)-R(x)}{D},$ (2.2)

where $\alpha=\alpha(t)$ is constant in space and
$R(x):=\int_{0}^{x}S(y)\,\mathrm{d}y.$ Another subsequent integration of (2.2)
in $\Omega$, exploiting the homogeneous Dirichlet boundary condition for $u$,
leads to

$\displaystyle\alpha(t)=\frac{\bigintsss_{0}^{1}\frac{R(y)}{D}\,\mathrm{d}y}{\bigintsss_{0}^{1}\frac{\,\mathrm{d}y}{D}}.$

Moreover, restricting the integration to $(0,x)$, we get the explicit
expression for $u$,

$\displaystyle
u=\alpha(t)\int_{0}^{x}\frac{\,\mathrm{d}y}{D(y)}-\int_{0}^{x}\frac{R(y)}{D}\,\mathrm{d}y.$
(2.3)

Performing a similar reasoning for (2.1b), we obtain

$\displaystyle-D\sigma_{x}+\beta(t)$
$\displaystyle=\Phi^{\prime\prime}(u)Du_{x}-\int_{0}^{x}\Phi^{\prime\prime}(u)(Du_{x})_{x}\,\mathrm{d}y$
$\displaystyle=\Phi^{\prime\prime}(u)Du_{x}+\int_{0}^{x}\Phi^{\prime\prime}(u)S(y)\,\mathrm{d}y,$

which can be reformulated in the more convenient way
$\sigma_{x}=\frac{\beta(t)}{D}-\Phi^{\prime\prime}(u)u_{x}-\frac{1}{D}\int_{0}^{x}\Phi^{\prime\prime}(u)S(y)\,\mathrm{d}y,$
with

$\displaystyle\beta(t)=\frac{\bigintsss_{0}^{1}\frac{1}{D}\bigintsss_{0}^{y}\Phi^{\prime\prime}(u)S(z)\,\mathrm{d}z\,\mathrm{d}y}{\bigintsss_{0}^{1}\frac{\,\mathrm{d}y}{D}},$

obtained through an integration of (2) and employing (1). Finally,
substitution of (2.2) and (2) into (2.1c) gives

$\displaystyle D_{t}$
$\displaystyle=-\frac{u_{x}}{D}\left(\int_{0}^{x}\Phi^{\prime\prime}(u)S(y)\,\mathrm{d}y-\beta(t)\right)$
$\displaystyle=\left(\frac{R(x)-\alpha(t)}{D^{2}}\right)\left(\int_{0}^{x}\Phi^{\prime\prime}(u)S(y)\,\mathrm{d}y-\beta(t)\right),$

which we recast as

$\displaystyle D^{2}D_{t}$
$\displaystyle=\left(\frac{\bigintsss_{0}^{1}\frac{R(x)-R(y)}{D}\,\mathrm{d}y}{\bigintsss_{0}^{1}\frac{\,\mathrm{d}y}{D}}\right)\left(\frac{\bigintsss_{0}^{1}\frac{\bigintsss_{0}^{x}\Phi^{\prime\prime}(u)S(z)\,\mathrm{d}z-\bigintsss_{0}^{y}\Phi^{\prime\prime}(u)S(z)\,\mathrm{d}z}{D}\,\mathrm{d}y}{\bigintsss_{0}^{1}\frac{\,\mathrm{d}y}{D}}\right)$
$\displaystyle=\left(\frac{\bigintsss_{0}^{1}\frac{\bigintsss_{y}^{x}S(z)\,\mathrm{d}z}{D}\,\mathrm{d}y}{\bigintsss_{0}^{1}\frac{\,\mathrm{d}y}{D}}\right)\left(\frac{\bigintsss_{0}^{1}\frac{\bigintsss_{y}^{x}\Phi^{\prime\prime}(u)S(z)\,\mathrm{d}z}{D}\,\mathrm{d}y}{\bigintsss_{0}^{1}\frac{\,\mathrm{d}y}{D}}\right).$
(2.4)

In order to give a meaning to the elliptic equations (2.1a–2.1b) we look for a
solution $D>0$.

It is immediate to notice that $u$ is bounded in terms of the $L^{1}$ norm of
$\frac{1}{D}$:

######  Proposition 1.

Let $u$ be given by equation (2.3). Then for all $t\geq 0$,

$\displaystyle\|u(\cdot,t)\|_{L^{\infty}(0,1)}\leq\|S\|_{L^{1}(0,1)}\bigg{\|}\frac{1}{D(\cdot,t)}\bigg{\|}_{L^{1}(0,1)}.$
(2.5)

###### Proof.

Since $u=0$ for $x\in\\{0,1\\}$, it follows that the function $|u|$ attains
its maximum at some point $x_{m}\in(0,1)$. Then, due to the differentiability
of $u$, we have $u_{x}=0$ at $x=x_{m}$, so that

$\displaystyle
R(x_{m})=\frac{\bigintsss_{0}^{1}\frac{R(y)}{D}\,\mathrm{d}y}{\bigintsss_{0}^{1}\frac{\,\mathrm{d}y}{D}}.$

Then, for some $\xi(t)\in(0,1)$ we get the following estimate

$\displaystyle|u(x_{m})|$
$\displaystyle=\bigg{|}R(x_{m})\bigintsss_{0}^{x_{m}}\frac{\,\mathrm{d}y}{D}-\bigintsss_{0}^{x_{m}}\frac{R(y)}{D}\,\mathrm{d}y\bigg{|}$
$\displaystyle=\bigg{|}\left(R(x_{m})-R(\xi(t)\right)\bigintsss_{0}^{x_{m}}\frac{\,\mathrm{d}y}{D}\bigg{|}$
$\displaystyle\leq\|S\|_{L^{1}(0,1)}\bigintsss_{0}^{1}\frac{1}{|D|}\,\mathrm{d}y,$

which concludes the proof.

Furthermore, if $\Phi^{\prime\prime}$ is bounded on $\mathbb{R}$, it is
possible to obtain a uniform bound (on every finite time interval) for
$\Phi^{\prime}(u)$ which is independent of
$\int_{0}^{1}\frac{\,\mathrm{d}y}{D(y)}$.

######  Proposition 2.

Let $S\in L^{1}(0,1)$, $D_{I}\in L^{1}_{+}(0,1)$ with (2) and
$\Phi^{\prime\prime}\in L^{\infty}(\mathbb{R})$. If $D(x,t)\geq 0$ for
$(x,t)\in[0,1]\times[0,T]$, with $0<T<\infty$, then there exists $C(T)>0$ such
that $\|\Phi^{\prime}(u(\cdot,t))\|_{L^{\infty}(0,1)}\leq C(T)\quad\forall
t\in[0,T].$

###### Proof.

Fix $\varepsilon>0$. Using (2.2) and the convexity of $\Phi$, we can recast
the energy functional (1.1) as

$\displaystyle E[D(t)]$
$\displaystyle=\int_{0}^{1}Du^{2}_{x}\Phi^{\prime\prime}(u)\,\mathrm{d}x$
$\displaystyle=\int_{0}^{1}|Du_{x}||u_{x}\Phi^{\prime\prime}(u)|\,\mathrm{d}x$
$\displaystyle=\int_{0}^{1}|\alpha(t)-R(x)|\big{|}\Phi^{\prime}(u)_{x}\big{|}\,\mathrm{d}x$
$\displaystyle=\int_{|\alpha(t)-R(x)|\leq\varepsilon}|\alpha(t)-R(x)|\big{|}\Phi^{\prime}(u)_{x}\big{|}\,\mathrm{d}x+\int_{|\alpha(t)-R(x)|>\varepsilon}|\alpha(t)-R(x)|\big{|}\Phi^{\prime}(u)_{x}\big{|}\,\mathrm{d}x$
$\displaystyle=:I_{\varepsilon_{-}}+I_{\varepsilon_{+}}.$

Since $D=D(t)$ represents a gradient flow of the energy $E$, we have
$E[D(t)]\leq E[D_{I}]<\infty$ for all $t\geq 0$, and the equality holds if and
only if $D\equiv D_{I}$. Therefore,

$\displaystyle
I_{\varepsilon_{+}}:=\int_{|\alpha(t)-R(x)|>\varepsilon}|\alpha(t)-R(x)|\big{|}\Phi^{\prime}(u)_{x}\big{|}\,\mathrm{d}x<E[D_{I}],$

and, consequently,

$\displaystyle\int_{|\alpha(t)-R(x)|>\varepsilon}\big{|}\Phi^{\prime}(u)_{x}\big{|}\,\mathrm{d}x<\frac{E[D_{I}]}{\varepsilon}.$
(2.6)

Now, take $x\in\\{z\in[0,1];\,|\alpha(t)-R(z)|\leq\varepsilon\\}$. Then, from
(2) we obtain the estimate

$\displaystyle\frac{1}{3}\big{|}(D^{3})_{t}\big{|}$
$\displaystyle=\bigg{|}\left[R(x)-\alpha(t)\right]\left(\int_{0}^{x}\Phi^{\prime\prime}\left(u\right)S(y)\,\mathrm{d}y-\beta(t)\right)\bigg{|}$
$\displaystyle\leq
2\varepsilon\|\Phi^{\prime\prime}\|_{L^{\infty}(0,1)}\|S\|_{L^{1}(0,1)}.$

Integrating the latter expression in time between $(0,t)$ and choosing
$\varepsilon=\frac{7}{48}\frac{\gamma^{3}}{\|\Phi^{\prime\prime}\|_{L^{\infty}(0,1)}\|S\|_{L^{1}(0,1)}T}$
we get the following lower bound for $D(x,t)$

$\displaystyle D^{3}(x,t)\geq
D_{I}(x)-6\varepsilon\|\Phi^{\prime\prime}(u)\|_{L^{\infty}(0,1)}\|S\|_{L^{1}(0,1)}t\geq\frac{1}{8}\gamma^{3}$

for all
$(x,t)\in\\{z\in[0,1];\,|\alpha(t)-R(z)|\leq\varepsilon\\}\times[0,T]$.
Therefore, from the latter inequality and (2.2) we find the following bound
for $u_{x}$,

$\displaystyle|u_{x}|\leq\frac{2\varepsilon}{\gamma}\qquad\mbox{for all
}(x,t)\in\\{z\in[0,1];\,|\alpha(t)-R(z)|\leq\varepsilon\\}\times[0,T],$

and

$\displaystyle\int_{|\alpha(t)-R(x)|\leq\varepsilon}\big{|}\Phi^{\prime}(u)_{x}\big{|}\,\mathrm{d}x\leq\frac{2\varepsilon\|\Phi^{\prime\prime}\|_{L^{\infty}(0,1)}}{\gamma}\qquad\mbox{for
all }t\in[0,T].$ (2.7)

A combination of (2.6) and (2.7) yields

$\displaystyle\int_{0}^{1}\big{|}\Phi^{\prime}(u)_{x}\big{|}\,\mathrm{d}x<\frac{E[D_{I}]}{\varepsilon}+\frac{2\varepsilon\|\Phi^{\prime\prime}\|_{L^{\infty}(0,1)}}{\gamma}\qquad\mbox{for
all }t\in[0,T].$

Finally, since $\Phi^{\prime}\left(u(x,t)\right)=\Phi^{\prime}(0)=0$ for
$x\in\\{0,1\\}$ we conclude that

$\displaystyle\|\Phi^{\prime}\left(u(\cdot,t)\right)\|_{W^{1,1}(0,1)}\leq\tilde{C}(T)\qquad\mbox{for
all }t\in[0,T],$

so that

$\displaystyle\|\Phi^{\prime}\left(u(\cdot,t)\right)\|_{L^{\infty}(0,1)}\leq
C(T)\qquad\mbox{for all }t\in[0,T].$

In particular, Proposition 2 implies that, if $\Phi^{\prime\prime}\in
L^{\infty}(\mathbb{R})$, then $\|u(\cdot,t)\|_{L^{\infty}(0,1)}$ is locally
bounded in $t$ if $\Phi^{\prime}(\infty)=\infty$ and
$\Phi^{\prime}(-\infty)=-\infty$.

The following theorem provides existence and uniqueness of a solution of the
integro-differential system (2.3–2).

######  Theorem 1.

The system (2.3–2) complemented with the initial condition (2) has a unique
positive solution $D\in C^{1}\left([0,T);C[0,1]\right)$. The maximal time of
existence $T>0$ is finite if and only if $\liminf_{t\rightarrow
T^{-}}min_{x\in[0,1]}D(x,t)=0$, otherwise the solution is global with
$T=+\infty$.

The proof of the theorem is a direct consequence of the Cauchy-Lipschitz-
Picard theorem. First, a bound for $u$ depends only on the minimum of $D$ from
Proposition 1. Employing Picard iterations (or Banach contraction) in the
space of continuous functions for a sufficiently small time $T$, such that $D$
remains bounded from below by a positive constant, which implies existence of
a local solution. The solution is global in time if $D=D(x,t)$ remains
positive for all $t\geq 0$.

To facilitate computations and gain deeper insight into the system, we now
replace the homogeneous Dirichlet boundary conditions with mixed Neumann-
Dirichlet,

$\displaystyle u_{x}(t,0)=0,$ $\displaystyle\quad u(t,1)=0,$
$\displaystyle\sigma_{x}(t,0)=0,$ $\displaystyle\quad\sigma(t,1)=0.$

This modification does not alter neither the gradient flow computations nor
Proposition 1 and Theorem 1. Moreover, the resulting integro-differential
system is a straightforward modification of (2), (2.3):

$\displaystyle
D^{2}D_{t}=R(x)\int_{0}^{x}\Phi^{\prime\prime}(u)S(y)\,\mathrm{d}y,$ (2.8a)
$\displaystyle u=\int_{x}^{1}\frac{R(y)}{D}\,\mathrm{d}y,$ (2.8b)

complemented with the initial condition (2), and the energy functional (1.1)
can be written as
$E[D]=\int_{0}^{1}\frac{R(x)^{2}}{D}\Phi^{\prime\prime}(u)\,\mathrm{d}x.$

We now focus on identifying conditions under which the diffusivity $D(x,t)$
remains positive for all $t>0$, implying global well-posedness of the problem
by Theorem 1. The following proposition provides some examples.

######  Proposition 3.

Let $D$ be the unique solution of (2.8a), (2.8b) with the initial condition
(2). If any of the following conditions hold:

  1. 1.

$S(x)\geq 0$ or $S(x)\leq 0$ almost everywhere on $(0,1)$,

  2. 2.

$R(x)\geq 0$ for $x\in[0,1]$ and $\Phi^{\prime\prime\prime}(u)\geq 0$, or
$R(x)\leq 0$ for $x\in[0,1]$ and $\Phi^{\prime\prime\prime}(u)\leq 0$.

  3. 3.

Define $S^{+}(x):=\max\\{S(x),0\\}$ and, resp., $S^{-}(x):=\max\\{-S(x),0\\}$
the positive and, resp., negative parts of $S=S(x)$. If there exists constants
$K_{+},K_{-}>0$ such that $K_{-}\leq\Phi^{\prime\prime}(u)\leq
K_{+}\,\,\forall u\in\mathbb{R}$ and
$\int_{0}^{x}S^{+}(y)\,\mathrm{d}y\geq\frac{K_{+}}{K_{-}}\int_{0}^{x}S^{-}(y)\,\mathrm{d}y,\quad\forall
x\in(0,1)$ or
$\int_{0}^{x}S^{+}(y)\,\mathrm{d}y\leq\frac{K_{-}}{K_{+}}\int_{0}^{x}S^{-}(y)\,\mathrm{d}y,\quad\forall
x\in(0,1),$

then there exists a $\gamma>0$ such that $D(x,t)\geq\gamma$ for all
$(x,t)\in[0,1]\times[0,\infty)$.

###### Proof.

If (1) holds, then trivially $D_{t}\geq 0$ and the result follows.
Assume now (2). Integrating (2.8a) by parts and using $u_{x}=-\frac{R(x)}{D}$,
we obtain

$\displaystyle D^{2}D_{t}$
$\displaystyle=R(x)\int_{0}^{x}\Phi^{\prime\prime}(u)R^{\prime}(y)\,\mathrm{d}y$
$\displaystyle=R(x)\left(\Phi^{\prime\prime}(u)R(x)-\int_{0}^{x}\Phi^{\prime\prime\prime}(u)u_{x}R(y)\,\mathrm{d}y\right)$
$\displaystyle=R^{2}(x)\Phi^{\prime\prime}(u)+R(x)\int_{0}^{x}\Phi^{\prime\prime\prime}(u)\frac{R(y)^{2}}{D}\,\mathrm{d}y\geq
0,$

which yields the statement. Finally, let (3) hold. Then,

$\displaystyle R(x)$
$\displaystyle=\int_{0}^{x}\big{(}S^{+}(y)-S^{-}(y)\big{)}\,\mathrm{d}y$
$\displaystyle\geq\left(\frac{K_{+}}{K_{-}}-1\right)\int_{0}^{x}S^{-}(y)\,\mathrm{d}y\geq
0\quad\forall x\in(0,1),$

and, consequently

$\displaystyle D^{2}D_{t}\geq
R(x)\int_{0}^{x}\big{(}K_{-}S^{+}(y)-K_{+}S^{-}(y)\big{)}\,\mathrm{d}y\geq
0\quad\forall x\in(0,1),$

which proves the proposition. An analogous argument can be made in the case
where (3) holds.

Let us remark that by a simple modification of the proof of Proposition 2, one
obtains $\|u(\cdot,t)\|_{L^{\infty}(0,1)}$ locally bounded in $t$ also for the
case of mixed boundary conditions, whenever $\Phi^{\prime\prime}\in
L^{\infty}(\mathbb{R})$ and $\Phi^{\prime}(\infty)=\infty$,
$\Phi^{\prime}(-\infty)=-\infty$.

To examine the convexity properties of the functional (2) under the constraint
(2.8), we need to calculate its second-order variation.

######  Proposition 4.

The second-order variation of (2) coupled to (2.8) reads

$\displaystyle\frac{\delta^{2}E[D^{0}]}{\delta
D^{2}}(D^{1},D^{1})=\bigintssss_{0}^{1}\frac{R^{2}(x)}{D^{0}}\bigg{[}$
$\displaystyle
2\Phi^{\prime\prime}(u^{0})\frac{(D^{1})^{2}}{(D^{0})^{2}}+2\Phi^{\prime\prime\prime}(u^{0})\left(\frac{D^{1}}{D^{0}}\bigintssss_{x}^{1}\frac{R(y)D^{1}}{(D^{0})^{2}}\,\mathrm{d}y+\bigintssss_{x}^{1}\frac{R(y)(D^{1})^{2}}{(D^{0})^{3}}\,\mathrm{d}y\right)$
$\displaystyle+\Phi^{\textrm{iv}}(u^{0})\left(\bigintssss_{x}^{1}\frac{R(y)D^{1}}{(D^{0})^{2}}\,\mathrm{d}y\right)^{2}\bigg{]}\,\mathrm{d}x.$
(2.9)

Moreover, if

  1. 1.

$2\Phi^{\prime\prime}(u)\Phi^{\textrm{iv}}(u)\geq\Phi^{\prime\prime\prime}(u)^{2}$
on $\mathbb{R}$,

  2. 2.

$R(x)\geq 0$ on $(0,1)$ and $\Phi^{\prime\prime\prime}(u)\geq 0$ on
$\mathbb{R}$ or $R(x)\leq 0$ on $(0,1)$ and $\Phi^{\prime\prime\prime}(u)\leq
0$ on $\mathbb{R}$,

then $E\equiv E[D]$ is convex on $L^{2}_{\gamma}(0,1):=\\{D\in
L^{2}(0,1)\,\mbox{ such that
}\,D\geq\gamma\,\mathrm{a.e.}\,\mathrm{on}\,(0,1)\\}$ for any $\gamma>0$.

###### Proof.

Fix $\varepsilon\in\mathbb{R}$. We expand $D=D^{0}+\varepsilon
D^{1}+O(\varepsilon^{2})$ with $D^{0}\geq\gamma>0$ $\mathrm{a.e.}$ on $(0,1)$.
Taylor expansion up to second order gives

$\displaystyle u[D^{0}+\varepsilon D^{1}]$
$\displaystyle=u[D^{0}]+\varepsilon\frac{\delta u[D^{0}]}{\delta
D}(D^{1})+\frac{\varepsilon^{2}}{2}\frac{\delta^{2}u[D^{0}]}{\delta
D}(D^{1},D^{1})+O(\varepsilon^{3})$ (2.10) $\displaystyle=:u^{0}+\varepsilon
u^{1}+\frac{\varepsilon^{2}}{2}u^{2}+O(\varepsilon^{3}).$ (2.11)

Plugging the above expansion into (2.1a), we can collect the $O(1)$ terms

$\displaystyle-\left(D^{0}u^{0}_{x}\right)_{x}=S,$ (2.12)

the first-order terms

$\displaystyle-\left(D^{0}u^{1}_{x}+D^{1}u^{0}_{x}\right)_{x}=0$

and the second-order ones

$\displaystyle-\left(\frac{1}{2}D^{0}u^{2}_{x}+D^{1}u^{1}_{x}\right)_{x}=0.$

Integration in $x$ of the latter one–dimensional conservation laws for
$u_{0},u_{1},u_{2}$ leads to the explicit expressions

$\displaystyle u^{0}_{x}=-\frac{R(x)}{D^{0}},$ $\displaystyle\quad
u^{0}=\bigintssss_{x}^{1}\frac{R(y)}{D^{0}}\,\mathrm{d}y$ (2.13a)
$\displaystyle u^{1}_{x}=\frac{R(x)D^{1}}{(D^{0})^{2}},$ $\displaystyle\quad
u^{1}=-\bigintssss_{x}^{1}\frac{R(y)D^{1}}{(D^{0})^{2}}\,\mathrm{d}y$ (2.13b)
$\displaystyle u^{2}_{x}=-2\frac{R(x)(D^{1})^{2}}{(D^{0})^{3}},$
$\displaystyle\quad
u^{2}=\bigintssss_{x}^{1}2\frac{R(y)(D^{1})^{2}}{(D^{0})^{3}}\,\mathrm{d}y.$
(2.13c)

Multiplication of (2.1a) by $\Phi^{\prime}(u)$ and integration by parts allow
us to recast the energy functional (2) as

$\displaystyle E[D]=\int_{0}^{1}S(x)\Phi^{\prime}(u)\,\mathrm{d}x.$

Thus, we have for the second variation of the functional

$\displaystyle\frac{\delta^{2}E[D^{0}]}{\delta D^{2}}(D^{1},D^{1})$
$\displaystyle=\frac{d^{2}}{d\varepsilon^{2}}\int_{0}^{1}S(x)\Phi^{\prime}\left(u^{0}+\varepsilon
u^{1}+\frac{\varepsilon^{2}}{2}u^{2}\right)\,\mathrm{d}x\bigg{|}_{\varepsilon=0}$
$\displaystyle=\frac{d^{2}}{d\varepsilon^{2}}\int_{0}^{1}S(x)\bigg{[}\Phi^{\prime}(u^{0})+\varepsilon\Phi^{\prime\prime}(u^{0})u^{1}+\frac{\varepsilon^{2}}{2}\left(\Phi^{\prime\prime}(u^{0})u^{2}+\Phi^{\prime\prime\prime}(u^{0})(u^{1})^{2}\right)\bigg{]}\,\mathrm{d}x\bigg{|}_{\varepsilon=0}$
$\displaystyle=\int_{0}^{1}S(x)\left(\Phi^{\prime\prime}(u^{0})u^{2}+\Phi^{\prime\prime\prime}(u^{0})(u^{1})^{2}\right)\,\mathrm{d}x.$

Using (2.12), integration by parts and substitution of (2.13) results in

$\displaystyle\frac{\delta^{2}E[D^{0}]}{\delta
D^{2}}(D^{1},D^{1})=\bigintssss_{0}^{1}\frac{R^{2}(x)}{D^{0}}\Bigg{[}$
$\displaystyle
2\Phi^{\prime\prime}(u^{0})\frac{(D^{1})^{2}}{(D^{0})^{2}}+2\Phi^{\prime\prime\prime}(u^{0})\left(\frac{D^{1}}{D^{0}}\bigintssss_{x}^{1}\frac{R(y)D^{1}}{(D^{0})^{2}}\,\mathrm{d}y+\bigintssss_{x}^{1}\frac{R(y)(D^{1})^{2}}{(D^{0})^{3}}\,\mathrm{d}y\right)$
$\displaystyle+\Phi^{\textrm{iv}}(u^{0})\left(\bigintssss_{x}^{1}\frac{R(y)D^{1}}{(D^{0})^{2}}\,\mathrm{d}y\right)^{2}\Bigg{]}\,\mathrm{d}x,$

which proves (4).

Observe that if (2) holds, we have

$\displaystyle\frac{\delta^{2}E[D^{0}]}{\delta
D^{2}}(D^{1},D^{1})\geq\bigintssss_{0}^{1}\frac{R^{2}(x)}{D^{0}}\Bigg{[}$
$\displaystyle
2\Phi^{\prime\prime}(u^{0})\frac{(D^{1})^{2}}{(D^{0})^{2}}+2\Phi^{\prime\prime\prime}(u^{0})\frac{D^{1}}{D^{0}}\bigintssss_{x}^{1}\frac{R(y)D^{1}}{(D^{0})^{2}}\,\mathrm{d}y$
$\displaystyle+\Phi^{\textrm{iv}}(u^{0})\left(\bigintssss_{x}^{1}\frac{R(y)D^{1}}{(D^{0})^{2}}\,\mathrm{d}y\right)^{2}\Bigg{]}\,\mathrm{d}x.$

Defining

$\displaystyle
v=\begin{bmatrix}\frac{D^{1}}{D^{0}}&\bigintssss_{x}^{1}\frac{R(y)D^{1}}{(D^{0})^{2}}\,\mathrm{d}y\end{bmatrix},\quad
Q=\begin{bmatrix}2\Phi^{\prime\prime}(u^{0})&\Phi^{\prime\prime\prime}(u^{0})\\\
\Phi^{\prime\prime\prime}(u^{0})&\Phi^{\textrm{iv}}(u^{0})\end{bmatrix},$

we write

$\displaystyle\frac{\delta^{2}E[D^{0}]}{\delta
D^{2}}(D^{1},D^{1})\geq\bigintssss_{0}^{1}\frac{R^{2}(x)}{D^{0}}v\cdot
Qv\,\mathrm{d}x.$

Finally, if (1) holds then $\det Q\geq 0$ and $E[D]$ is convex on
$L^{2}_{\gamma}(0,1)$ for every $\gamma>0$.

When considering the case of $\Phi(u)=u^{p}$, the energy functional is convex
if $1\leq p\leq 2$ and $R\leq 0$, or if $p\geq 4$ and $R\geq 0$. Condition (1)
is reminiscent, but not identical, of the one in [1]. We emphasize that
Condition 2 satisfies the hypothesis of Proposition 3, ensuring that $D$ never
reaches the boundary of $L^{2}_{\gamma}(0,1)$ since $D_{t}\geq 0$.

##  3 A remarkable example

Based on Proposition 3, we have identified certain cases where the diffusion
$D$ remains positive over the entire domain so that the problem is well-posed
thanks to Theorem 1. Nevertheless, it is important to note that these cases
represent only a partial enumeration of the possible scenarios, and there may
exist other cases where D remains positive throughout the domain or it becomes
zero somewhere. By systematically exploring the range of possible scenarios,
we can gain a deeper understanding of the behavior of $D$ and its implications
for the system under study. This involves analyzing the system subject to
different distributions of sources and sinks of $S(x)$. With this in mind, let
us define $S(x):=a\delta(x-x_{0})-b\delta(x-x_{1}),\quad
a,b>0,\,\,0<x_{0}<x_{1}<1,$ where $\delta(x-\bar{x})$ is the Dirac delta
distribution centered at $\bar{x}$. Then
$R(x)=a\chi_{[x_{0},x_{1}]}(x)+(a-b)\chi_{[x_{1},1]}(x)\quad\forall
x\in[0,1].$ Substituting (3), (3) into (2.8a), we can recast the evolution
equation for $D$ as

$\displaystyle D^{2}D_{t}$
$\displaystyle=\left[a\chi_{[x_{0},x_{1}]}(x)+(a-b)\chi_{[x_{1},1]}(x)\right]\int_{0}^{x}\Phi^{\prime\prime}(u)\big{(}a\delta(y-x_{0})-b\delta(y-x_{1})\big{)}\,\mathrm{d}y,$

i.e.,

$\displaystyle D^{2}D_{t}=\left\\{\begin{array}[]{ll}\displaystyle
0&x\in[0,x_{0}),\\\
a^{2}\,\Phi^{\prime\prime}(u(x_{0},t))&x\in[x_{0},x_{1}),\\\
(a-b)\big{[}a\Phi^{\prime\prime}(u(x_{0},t))-b\Phi^{\prime\prime}(u(x_{1},t))\big{]}&x\in[x_{1},1].\end{array}\right.$
(3.4)

For the sake of simplicity of presentation we choose a piecewise constant
initial datum:

$\displaystyle{D(0,x)=}\left\\{\begin{array}[]{ll}\displaystyle
D_{I,0}&x\in[0,x_{0}),\\\ D_{I,1}&x\in[x_{0},x_{1}),\\\
D_{I,2}&x\in[x_{1},1],\end{array}\right.$ (3.8)

with $D_{I,0},D_{I,1},D_{I,2}>0$. We immediately notice that in the
subintervals $[0,x_{0})$, $[x_{0},x_{1})$ and $[x_{1},1]$ the evolution of $D$
does not depend on $x$. Hence

$\displaystyle{D(t,x)=}\left\\{\begin{array}[]{ll}\displaystyle
D_{I,0}&x\in[0,x_{0}),\\\ D_{1}(t)&x\in[x_{0},x_{1}),\\\
D_{2}(t)&x\in[x_{1},1],\end{array}\right.$ (3.12)

where $D_{1}(t)\geq\gamma>0$ $\forall t\in[0,\infty)$ from (3.4).
Consequently, we rewrite (2.8b) as

$\displaystyle{u=}\left\\{\begin{array}[]{ll}\displaystyle\frac{a(x_{1}-x_{0})}{D_{1}(t)}+\frac{(a-b)(1-x_{1})}{D_{2}(t)},&x\in[0,x_{0}),\\\
\displaystyle\frac{a(x_{1}-x)}{D_{1}(t)}+\frac{(a-b)(1-x_{1})}{D_{2}(t)}&x\in[x_{0},x_{1}),\\\
\displaystyle\frac{(a-b)(1-x)}{D_{2}(t)}&x\in[x_{1},1],\end{array}\right.$
(3.16)

where we are interested in

$\displaystyle u(x_{0},t)$
$\displaystyle=\frac{a(x_{1}-x_{0})}{D_{1}(t)}+\frac{(a-b)(1-x_{1})}{D_{2}(t)},$
(3.17a) $\displaystyle u(x_{1},t)$
$\displaystyle=\frac{(a-b)(1-x_{1})}{D_{2}(t)},$ (3.17b)

and we notice that $u(x_{0},t)>u(x_{1},t)$.

As we previously noted, existence and uniqueness of solutions is guaranteed
whenever the function $D_{2}$ remains positive. Therefore, it is crucial to
investigate conditions that lead to this result. To achieve that, we will
leverage the properties of the gradient flow structure of our system.

######  Proposition 5.

Let $D$ be the solution of (3.4) with associated uniformly positive initial
condition.

In the case $a>b$ we have the following results.

  1. 1.

If
$\int_{0}^{\infty}\Phi^{\prime\prime}(u)\,\mathrm{d}u\equiv\lim_{u\rightarrow\infty}\Phi^{\prime}(u)=\infty$,
then $D(x,t)>0$ $\forall(x,t)\in[0,1]\times[0,\infty)$.

  2. 2.

If
$\int_{0}^{\infty}\Phi^{\prime\prime}(u)\,\mathrm{d}u\equiv\lim_{u\rightarrow\infty}\Phi^{\prime}(u)=:\alpha<\infty$
and $\alpha\geq\frac{E[D_{I}]}{(a-b)},$ then $D(x,t)>0$
$\forall(x,t)\in[0,1]\times[0,\infty)$.

We have the following results in the complementary case $a<b$.

  1. 3.

If
$\int_{-\infty}^{0}\Phi^{\prime\prime}(u)\,\mathrm{d}u\equiv-\lim_{u\rightarrow-\infty}\Phi^{\prime}(u)=\infty$,
then $D(x,t)>0$ $\forall(x,t)\in[0,1]\times[0,\infty)$.

  2. 4.

If
$\int_{-\infty}^{0}\Phi^{\prime\prime}(u)\,\mathrm{d}u\equiv-\lim_{u\rightarrow-\infty}\Phi^{\prime}(u)=:\beta<\infty$
and $\beta\geq\frac{E[D_{I}]}{(b-a)},$ then $D(x,t)>0$
$\forall(x,t)\in[0,1]\times[0,\infty)$.

###### Proof.

From the gradient flow theory, it follows that $E[D(t)]\leq
E[D_{I}]<\infty,\,\,\forall t\geq 0$ and the equality holds if and only if
$D=D_{I}\,\,\forall t\geq 0$. Therefore, without loss of generality, we assume
$E[D(t)]<E[D_{I}]$. Further, we can expand the energy functional (2) as

$\displaystyle E[D(t)]$
$\displaystyle=\int_{x_{0}}^{x_{1}}\frac{a^{2}}{D_{1}(t)}\Phi^{\prime\prime}(u)\,\mathrm{d}x+\int_{x_{1}}^{1}\frac{(a-b)^{2}}{D_{2}(t)}\Phi^{\prime\prime}(u)\,\mathrm{d}x$
$\displaystyle=\int_{x_{0}}^{x_{1}}\frac{a^{2}}{D_{1}(t)}\Phi^{\prime\prime}\left(\frac{a(x_{1}-x)}{D_{1}(t)}+\frac{(a-b)(1-x_{1})}{D_{2}(t)}\right)\,\mathrm{d}x+\int_{x_{1}}^{1}\frac{(a-b)^{2}}{D_{2}(t)}\Phi^{\prime\prime}\left(\frac{(a-b)(1-x)}{D_{2}(t)}\right)\,\mathrm{d}x$
$\displaystyle=:I_{1}(t)+I_{2}(t).$

Through a change of variable, it is possible to recast $I_{1}(t),I_{2}(t)$ as

$\displaystyle\displaystyle
I_{1}(t)=a\int_{\frac{(a-b)(1-x_{1})}{D_{2}(t)}}^{\frac{a(x_{1}-x_{0})}{D_{1}(t)}+\frac{(a-b)(1-x_{1})}{D_{2}(t)}}\Phi^{\prime\prime}(y)\,\mathrm{d}y,$

and

$\displaystyle\displaystyle
I_{2}(t)=(a-b)\int_{0}^{\frac{(a-b)(1-x_{1})}{D_{2}(t)}}\Phi^{\prime\prime}(y)\,\mathrm{d}y.$

Assume the proposition is false, i.e., $\exists\,0<T<\infty$ such that
$D_{2}(T)=0$.
Let $a>b$. (1) is easily proved since $E[D(t)]=\infty<E[D_{I}],$ which is a
contradiction. To prove (2), we first notice that from hypothesis
$\Phi^{\prime\prime}\in L^{1}(0,\infty)$, so that $I_{1}(T)=0$. Then

$\displaystyle(a-b)\alpha=E[D(t)]<E[D_{I}]$

leads to a contradiction of (2).
By applying analogous reasoning in the case where $a<b$, we arrive at a series
of contradictions that establish the validity of both (3) and (4).

We observe that in the case $a>b$, we have
$E[D(t)]=a\int_{0}^{\frac{a(x_{1}-x_{0})}{D_{1}(t)}+\frac{(a-b)(1-x_{1})}{D_{2}(t)}}\Phi^{\prime\prime}(y)\,\mathrm{d}y-b\int_{0}^{\frac{(a-b)(1-x_{1})}{D_{2}(t)}}\Phi^{\prime\prime}(y)\,\mathrm{d}y<a\alpha\quad\forall
t\geq 0.$ Hence, proposition 5 does not cover the cases
$(a-b)\alpha<E[D_{I}]<a\alpha,$ for which we currently do not have results. On
the other hand, whenever $a<b$, we write
$E[D(t)]=(b-a)\int_{-\frac{(b-a)(1-x_{1})}{D_{2}(t)}}^{0}\Phi^{\prime\prime}(y)\,\mathrm{d}y+a\int_{-\frac{(b-a)(1-x_{1})}{D_{2}(t)}}^{\frac{a(x_{1}-x_{0})}{D_{1}(t)}-\frac{(b-a)(1-x_{1})}{D_{2}(t)}}<b\beta+a\alpha,$
which leaves open the case $(b-a)\beta<E[D_{I}]<b\beta+a\alpha.$

######  Remark 1.

The work [15] highlights an important convex entropy generating function
$\Phi(u)=(u+1)\left(\ln(u+1)-1\right)$ for $u>-1$. This function transforms
the energy functional (1.1) into the Fisher information. When considering the
one–dimensional setting and selecting $S(x)$ as in (3), it is necessary to
ensure that $a>b$ to maintain $u>-1$ for arbitrary values of $D_{1}$, $D_{2}$.
Consequently, Condition (1) can be used to establish the global well-posedness
of the Fisher problem.

Let us now consider the evolutionary equation for $D_{2}(t)$, i.e., (3.4).
Dividing it by $\Phi^{\prime\prime}(u(x_{1},t))$, we obtain

$\displaystyle\frac{D^{2}_{2}(D_{2})_{t}}{\Phi^{\prime\prime}\left(\frac{(a-b)(1-x_{1})}{D_{2}(t)}\right)}=(a-b)\left[a\frac{\Phi^{\prime\prime}\left(\frac{a(x_{1}-x_{0})}{D_{1}(t)}+\frac{(a-b)(1-x_{1})}{D_{2}(t)}\right)}{\Phi^{\prime\prime}\left(\frac{(a-b)(1-x_{1})}{D_{2}(t)}\right)}-b\right]$

and integration between $0$ and $t$ yields

$\displaystyle\bigintsss_{D_{I,2}}^{D_{2}(t)}\frac{y^{2}\,\mathrm{d}y}{\Phi^{\prime\prime}\left(\frac{(a-b)(1-x_{1})}{y}\right)}$
$\displaystyle=(a-b)\left[a\bigintss_{0}^{t}\frac{\Phi^{\prime\prime}\left(\frac{a(x_{1}-x_{0})}{D_{1}(s)}+\frac{(a-b)(1-x_{1})}{D_{2}(s)}\right)}{\Phi^{\prime\prime}\left(\frac{(a-b)(1-x_{1})}{D_{2}(s)}\right)}\,\mathrm{d}s-bt\right]$
$\displaystyle=:(a-b)(aI(t)-bt),$ (3.18)

where

$\displaystyle
I(t):=\bigintss_{0}^{t}\frac{\Phi^{\prime\prime}\left(\frac{a(x_{1}-x_{0})}{D_{1}(s)}+\frac{(a-b)(1-x_{1})}{D_{2}(s)}\right)}{\Phi^{\prime\prime}\left(\frac{(a-b)(1-x_{1})}{D_{2}(s)}\right)}\,\mathrm{d}s.$
(3.19)

Here we have replaced $u(x_{0},t)$, $u(x_{1},t)$ by (3.17).

We have the following collection of results.

######  Proposition 6.

Let $D$ be the solution of (3.4) with associated uniformly positive initial
condition and $\Phi$ be strictly convex.

  1. 1.

If $a>b$ and

$\displaystyle\int_{1}^{\infty}\frac{\,\mathrm{d}v}{\Phi^{\prime\prime}(v)v^{4}}=\infty,$
(3.20)

then $D(x,t)>0$ $\forall(x,t)\in[0,1]\times[0,\infty)$.

  2. 2.

If $a<b$,
$\limsup_{w\rightarrow-\infty}\frac{\Phi^{\prime\prime}(z+w)}{\Phi^{\prime\prime}(w)}$
is a locally bounded function of $z\in[0,\infty)$ and

$\displaystyle\int_{1}^{\infty}\frac{\,\mathrm{d}v}{\Phi^{\prime\prime}(-v)v^{4}}=\infty,$
(3.21)

then $D(x,t)>0$ $\forall(x,t)\in[0,1]\times[0,\infty)$.

  3. 3.

If $a>b$ and

$\displaystyle\liminf_{w\rightarrow\infty}\frac{\Phi^{\prime\prime}(z+w)}{\Phi^{\prime\prime}(w)}>\frac{b}{a}$
(3.22)

uniformly for $z$ in compact subsets of $(0,\infty)$, then $D(x,t)>0$
$\forall(x,t)\in[0,1]\times[0,\infty)$.

  4. 4.

If $a<b$ and

$\displaystyle\limsup_{w\rightarrow-\infty}\frac{\Phi^{\prime\prime}(z+w)}{\Phi^{\prime\prime}(w)}<\frac{b}{a}$
(3.23)

uniformly for $z$ in compact subsets of $(0,\infty)$, then $D(x,t)>0$
$\forall(x,t)\in[0,1]\times[0,\infty)$.

###### Proof.

Assume $D_{2}(t)>0$ on $[0,T)$ and $D_{2}(T)=0$. Consequently, we rewrite
(3.18) as

$\displaystyle-\bigintsss_{0}^{D_{I,2}}\frac{y^{2}\,\mathrm{d}y}{\Phi^{\prime\prime}\left(\frac{(a-b)(1-x_{1})}{y}\right)}=(a-b)\left(aI(T)-bT\right),$
(3.24)

where $I(T)>0$ from the hypothesis of strictly convexity of $\Phi$. On the
other hand, we notice that

$\displaystyle-\bigintsss_{0}^{D_{I,2}}\frac{y^{2}\,\mathrm{d}y}{\Phi^{\prime\prime}\left(\frac{(a-b)(1-x_{1})}{y}\right)}=-\bigintss_{\frac{(a-b)(1-x_{1})}{D_{I,2}}}^{\infty}\frac{(a-b)^{3}(1-x_{1})^{3}}{\Phi^{\prime\prime}(v)v^{4}}\,\mathrm{d}v=-\infty$

from (3.20). Therefore, we should have $I(T)=-\infty$ which contradicts the
positivity of this term. Hence, $D_{2}(t)>0\,\,\forall t>0$ and (1) is proved.
To prove (2), we use a similar argument. Indeed, (3.24) holds and

$\displaystyle-\bigintsss_{0}^{D_{I,2}}\frac{y^{2}\,\mathrm{d}y}{\Phi^{\prime\prime}\left(-\frac{(b-a)(1-x_{1})}{y}\right)}=-\bigintss_{\frac{(b-a)(1-x_{1})}{D_{I,2}}}^{\infty}\frac{(b-a)^{3}(1-x_{1})^{3}}{\Phi^{\prime\prime}(-v)v^{4}}\,\mathrm{d}v=-\infty,$

which implies that $I(T)=\infty$. Next, let us define

$\displaystyle w(t)$ $\displaystyle:=-\frac{(b-a)(1-x_{1})}{D_{2}(t)},$
$\displaystyle z(t)$ $\displaystyle:=\frac{a(x_{1}-x_{0})}{D_{1}(t)}.$

Clearly, $z(t)\leq z(0)<\infty$, and from the definition (3.19) of $I(t)$
combined with the bounded limsup hypothesis we must have $I(T)<\infty$, which
leads to a contradiction.
From (3.22), (3.19) we immediately notice that the right-hand side of (3.24)
must be positive, and, conversely, the left-hand side is negative from the
hypothesis of strict convexity of $\Phi(u)$ and this gives the contradiction
to prove (3).

The statement (4) is proved using an analogous argument.

##  4 Numerical results

In this section we present several numerical simulations of the spatially
one–dimensional systems (2.8a–2.8b) and of the ODE presented in (3.4), where
the solution only exists locally in time due to $D=D(x,t)$ touching zero in
finite time. We make use of a $3^{\rm rd}$ order in time numerical scheme, and
we consider several entropy functions $\Phi$ to observe in which cases
$D=D(x,t)$ vanishes in finite time.

A semi–implicit discretization is adopted to achieve high order of accuracy in
time. In particular, we make use of implicit-explicit (IMEX) Runge–Kutta
schemes [17, 9], which are multi-step methods based on s-stages and typically
represented with the double Butcher tableau,

$\widetilde{c}$ | $\widetilde{A}$
---|---
| $\widetilde{b}^{\top}$

$c$ | $\widehat{A}$
---|---
| $\widehat{b}^{\top}$

with the matrices $(\widetilde{A},\widehat{A})\in\mathbb{R}^{s\times s}$ and
the vectors
$(\widetilde{b},\widehat{b},\widetilde{c},\widehat{c})\in\mathbb{R}^{s}$ where
$s>0$ is the number of stages of the Runge–Kutta method. The tilde symbol
refers to the explicit scheme, and $\widetilde{A}=\\{\widetilde{a}_{i,j}\\}$
is a lower triangular matrix with zero elements on the diagonal, while
$\widehat{A}=\\{\widehat{a}_{i,j}\\}$ is a triangular matrix which accounts
for the implicit scheme, thus having non-zero elements on the diagonal. Here,
we adopt IMEX schemes with $\widetilde{b}=\widehat{b}$, and the stiffly
accurate property in the implicit part.

First, we discretize in space and time the system in (2.8a–2.8b), and after we
focus on the case (3) in which the sink/source distribution $S=S(x)$ is given
by a superposition of Dirac delta distributions.

###  4.1 Discretization in space and time

In this section we discretize in space and time the spatially one–dimensional
system (2.8a–2.8b), which we recall here for the sake of the reader,

$\displaystyle D^{2}D_{t}$ $\displaystyle=R(x)V(x),\quad
x\in\Omega,\,t\in[0,t_{\rm fin}]$ (4.1) $\displaystyle u$
$\displaystyle=\int_{x}^{1}\frac{R(y)}{D}\,dy$ (4.2)

where $V(x)=\int_{0}^{x}\Phi^{\prime\prime}(u)S(y)\,dy$ and
$R(x)=\int_{0}^{x}S(y)\,dy$. The expression that we choose for the sink/source
function, $S(x)$, is the following

$\displaystyle S(x)=mx+q,\quad m,q\in\mathbb{R},$ (4.3)

where the parameters $m$ and $q$ are chosen in such a way the function
$S=S(x)$ changes sign in $[0,1]$, while its primitive $R=R(x)$ remains non-
negative for every $x$. To close the system, we prescribe the initial
condition $D(t=0)=D_{I}\in\mathbb{R}$.

For a fixed $N\in\mathbb{N}$ we discretize the spatial domain $\Omega=[0,1]$
using the equidistant grid x = { $\textbf{x}_{i}=ih;\,i=0,\cdots,N;\,h=1/N$}.
The discrete quantities shall be represented by bold letters. In particular,
${\bf S}=\\{{\bf S}_{i}\approx S({\bf x}_{i}),i=1,\cdots,N\\}$ represents the
vector that approximates the function $S(x)$ in the points ${\bf x}_{i}$.
Moreover, we define

$\displaystyle\textbf{R}=m\frac{\textbf{x}^{2}}{2}+q\textbf{x}$ $\displaystyle
0=\textbf{u}_{N}=\textbf{u}({\bf x}_{N})\approx u(x=1)$ (4.4)
$\displaystyle\textbf{u}_{N-i}=\sum_{j=1}^{i}\textbf{I}^{u}_{N-j}\,\quad
i=1,\cdots,N$ (4.5) $\displaystyle 0={\bf V}_{0}=V({\bf x}_{0})\approx V(x=0)$
(4.6) $\displaystyle\textbf{V}_{i}=\sum_{j=1}^{i}\textbf{I}^{V}_{j}\,\quad
i=1,\cdots,N$ (4.7)

where
$\textbf{I}^{u}_{i}=(\textbf{R}_{i}/\textbf{D}_{i}+\textbf{R}_{i+1}/\textbf{D}_{i+1})h/2$
and
$\textbf{I}^{V}_{i}=({\Phi^{\prime\prime}}_{i}\textbf{S}_{i}+{\Phi^{\prime\prime}}_{i-1}\textbf{S}_{i-1})h/2$.

Now we multiply and divide (4.1) by D

$\textbf{D}_{t}=M(\textbf{D})\textbf{D}$ (4.8)

where

$M(\textbf{D})=\frac{\textbf{R}(x)\textbf{V}(x)}{\textbf{D}^{3}}.$

To apply the partitioned Runge–Kutta method to (4.8), let us first set
$\textbf{D}^{1}_{E}=\textbf{D}^{n}$, then the stage fluxes for $i=1,\cdots,s$
are calculated as

$\displaystyle\textbf{D}_{E}^{i}$ $\displaystyle=\textbf{D}^{n}+\Delta
t\,\sum_{j=1}^{\rm
s-1}\widetilde{a}_{i,j}M(\textbf{D}_{E}^{j})\textbf{D}_{I}^{j},\quad
i=1,\cdots,s$ (4.9) $\displaystyle\textbf{D}_{I}^{i}$
$\displaystyle=\textbf{D}^{n}+\Delta t\,\sum_{j=1}^{\rm
s}\widehat{a}_{i,j}M(\textbf{D}_{E}^{j})\textbf{D}_{I}^{j},\quad i=1,\cdots,s$
(4.10)

and the numerical solution is finally updated with

$\displaystyle\textbf{D}^{n+1}$ $\displaystyle=\textbf{D}^{n}+\Delta
t\sum_{i=1}^{\rm s}\widehat{b}(i)M(\textbf{D}_{E}^{i})\textbf{D}_{I}^{i}$
(4.11)

where $\rm s$ is the number of stages of the scheme, and $\Delta t>0$ the time
step.

####  4.1.1 Accuracy test and Results

In this section we first check the accuracy of the discretization. We start
with a simple case, with a specific choice of the second derivative of the
entropy function, $\Phi^{\prime\prime}(u)=1$, which enables us to find an
exact solution (see Fig. 1 (a)). Moreover, we test the accuracy of the time
discretization in a general setting, with a reference solution calculated with
a fine grid (see Fig. 1 (b)).

A generic IMEX Runge–Kutta scheme is described with a triplet
$(\widehat{s},\widetilde{s},p)$ which characterizes the number $\widehat{s}$
of stages of the implicit scheme, the number $\widetilde{s}$ of stages of the
explicit scheme and the order $p$ of the resulting scheme.

A Runge–Kutta scheme is of order 3 if and only if the following conditions are
satisfied [14, Theorem 2.1]:

$\displaystyle\sum_{j}b_{j}$ $\displaystyle=1,\quad$ $\displaystyle
2\sum_{j,k}b_{j}a_{j,k}$ $\displaystyle=1,$ (4.12) $\displaystyle
3\sum_{j,k,l}b_{j}a_{j,k}a_{j,l}$ $\displaystyle=1,\quad$ $\displaystyle
6\sum_{j,k,l}b_{j}a_{j,k}a_{k,l}$ $\displaystyle=1.$ (4.13)

As third semi–implicit Runge–Kutta methods that satisfies the set of order
conditions (4.12–4.13) is given by the IMEX-SSP3(4,3,3) L-stable Scheme: SSP-
LDIRK3(4,3,3) with $\widehat{b}=\widetilde{b}$, i.e.

0 | 0 | 0 | 0 | 0
---|---|---|---|---
0 | 0 | 0 | 0 | 0
1 | 0 | 1 | 0 | 0
1/2 | 0 | 1/4 | 1/4 | 0
| 0 | 1/6 | 1/6 | 2/3

$\lambda$ | $\lambda$ | 0 | 0 | 0
---|---|---|---|---
0 | -$\lambda$ | $\lambda$ | 0 | 0
1 | 0 | $1-\lambda$ | $\lambda$ | 0
1/2 | $\mu$ | $\eta$ | $1/2-\mu-\eta-\lambda$ | $\lambda$
| 0 | 1/6 | 1/6 | 2/3

with $\lambda=0.24169426078821$, $\mu=\lambda/4$ and $\eta=0.12915286960590$
[9].

\begin{overpic}[abs,width=433.62pt,unit=1mm,scale={.25}]{Figures/accuracy_exact_solu}
\put(1.0,54.0){(a)} \end{overpic}

\begin{overpic}[abs,width=433.62pt,unit=1mm,scale={.25}]{Figures/error_plot_D}
\put(1.0,54.0){(b)} \end{overpic}

Figure 1: (a): Space and time accuracy plot of the numerical scheme defined in
(4.4–4.11) at final time $t=0.1$. The expression chosen for the second
derivative of the entropy function is ${\bf\Phi}^{\prime\prime}({\bf u})=1$,
while the Butcher Tableau are defined in SSP-LDIRK3(4,3,3) scheme. (b). Time
accuracy of the numerical scheme defined in (4.4–4.11) at final time $t=0.01$.
The expression chosen for the second derivative of the entropy function is
${\bf\Phi}^{\prime\prime}({\bf u})=(u^{2}(\sin(u)+\sigma)+1)^{-1}$, and
$\sigma$ is chosen such that the convexity of ${\bf\Phi}^{\prime\prime}({\bf
u})$ is maintained, i.e., $\sigma=2$. In both panels $m=-1.98$ and $q=1$ in
the expression for the sink/source function in (4.3).

In Fig. 1 we show the order of accuracy of the numerical scheme, defined in
(4.4–4.11). To do this, we calculate a reference solution, and then we compute
the relative error as follows

$\displaystyle{\rm error}=\frac{||{\bf D}_{h,\Delta t}-{\bf D}_{\rm
ref}||_{2}}{||{\bf D}_{\rm ref}||_{2}}$ (4.14)

where ${\bf D}_{h,\Delta t}$ is the numerical solution with space step equal
to $h$ and time step equal to $\Delta t$, while ${\bf D}_{\rm ref}$ is
calculated in two different ways. In panel (a), we consider a simple case, in
which the entropy function is defined as $\Phi^{\prime\prime}(u)=1$. In this
case we have an explicit expression for the exact solution of the system
(4.1–4.2), and in order to test the accuracy in space and time, we consider
$h=\Delta t\in\\{0.04,0.02,0.01,0.005,0.0025\\}$ for ${\bf D}_{h,\Delta t}$.
In panel (b), we fix $h=10^{5}$, and we calculate the time accuracy, fixing a
reference time step, $\Delta t_{\rm ref}=10^{-7}$, and a corresponding
reference solution ${\bf D}_{\rm ref}$. For ${\bf D}_{h,\Delta t}$ we choose
$\Delta t=0.1\cdot 2^{-k},\,k\in\\{2,3,4,5,6,7\\}$.

\begin{overpic}[abs,width=433.62pt,unit=1mm,scale={.25}]{Figures/D_phi_exponential1}
\put(1.0,54.0){(a)} \end{overpic}

\begin{overpic}[abs,width=433.62pt,unit=1mm,scale={.25}]{Figures/D_phi_exponential01_new}
\put(1.0,54.0){(b)} \end{overpic}

\begin{overpic}[abs,width=433.62pt,unit=1mm,scale={.25}]{Figures/D_phi_exponential01_zoom_new}
\put(1.0,54.0){(c)} \end{overpic}

Figure 2: Plot of the solution D (left part of each plot, in blue), together
with the function ${\bf\Phi}^{\prime\prime}({\bf u})$ (right part of each
plot, in orange), for different times. The expression of the second derivative
of the entropy function is ${\bf\Phi}^{\prime\prime}({\bf u})=\exp(-u)$, the
initial condition is D(t=0) = 1 (a), D(t=0) = 0.1 (b), and a zoom-in of panel
(b) in (c), while $m=-1.98$ and $q=1$ in the expression (3) and $\Delta
t=10^{-6},h=10^{-4}$.

Figs. 3–4 are equivalent to Fig. 2, with different expressions for the entropy
function. In Fig. 3 $\Phi^{\prime\prime}(u)=(u+1)^{-1}$, while in Fig. 4
$\Phi^{\prime\prime}(u)=(u^{2}(\sin(u)+\sigma)+1)^{-1}$. We observe that, with
the right choice of the initial condition, i.e. ${\bf D}(t=0)=0.1$, the
solution touches zero in finite time. Here, again, we choose as final time a
few time steps prior the function $\bf D$ touches zero, due to the presence of
spurious numerical oscillations arising from the division by $\bf D$.

In Fig. 5 we show that it is easier for the numerical scheme to compute the
solution $\bf D$ that approaches to zero when we refine the time step. Since
there is a division by $\bf D$ in the expression of $\bf u$ (see, e.g.,
(2.3)), spurious oscillations start to appear in the numerical simulations
when $\bf D$ is close to zero and, for this reason, we are not able to show
the exact time when $\bf D$. Fig. 5 confirms that, even if we are not able to
plot that moment, the trend of the numerical results are in agreement when
$\Delta\to 0$.

\begin{overpic}[abs,width=433.62pt,unit=1mm,scale={.25}]{Figures/D_phi_Fisher1}
\put(1.0,54.0){(a)} \end{overpic}

\begin{overpic}[abs,width=433.62pt,unit=1mm,scale={.25}]{Figures/D_phi_Fisher01_new}
\put(1.0,54.0){(b)} \end{overpic}

\begin{overpic}[abs,width=433.62pt,unit=1mm,scale={.25}]{Figures/D_phi_Fisher01_zoom_new}
\put(1.0,54.0){(c)} \end{overpic}

Figure 3: Plot of the solution D (left part of each plot, in blue), together
with the function $\bf\Phi^{\prime\prime}(u)$ (right part of each plot, in
orange), for different times. The second derivative of the entropy function is
${\bf\Phi}^{\prime\prime}({\bf u})=(u+1)^{-1}$, and the initial condition is
D(t=0) = 1 (a), D(t=0) = 0.1 (b), and zoom-in (c), while $m=-1.98$ and $q=1$
in the expression (3) and $\Delta t=h=10^{-3}$.

\begin{overpic}[abs,width=433.62pt,unit=1mm,scale={.25}]{Figures/D_phi_u2sin1}
\put(1.0,54.0){(a)} \end{overpic}

\begin{overpic}[abs,width=433.62pt,unit=1mm,scale={.25}]{Figures/D_phi_u2sin01_new}
\put(1.0,54.0){(b)} \end{overpic}

\begin{overpic}[abs,width=433.62pt,unit=1mm,scale={.25}]{Figures/D_phi_u2sin01_zoom_new}
\put(1.0,54.0){(c)} \end{overpic}

Figure 4: Plot of the solution D (left part of each plot, in blue), together
with the function $\bf\Phi^{\prime\prime}(u)$ (right part of each plot, in
orange), for different times. The second derivative of the entropy function is
$\Phi^{\prime\prime}(u)=(u^{2}(\sin(u)+\sigma)+1)^{-1}$, and the initial
condition is D(t=0) = 1 (a), D(t=0) = 0.1 (b), and zoom-in (c), while
$m=-1.98$ and $q=1$ in the expression (3) and $\Delta t=10^{-6}$ and
$h=10^{-4}$.

In Fig. 6 we consider a larger domain, i.e., $x\in[0,1.5]$, to show that the
point in which the solution touches zero does not depend neither on the
domain, nor on its boundary. On the contrary, it depends on the definition of
the function $S(x)$ and on the function $\bf\Phi^{\prime\prime}(u)$. In Fig. 6
$m=-1.98/1.5$ and $q=1$ in (4.3).

Figure 5: In this plot we show the minimum of the function $\bf D$ at final
time, changing the time step $\Delta t$, for different expressions of
${\bf\Phi}^{\prime\prime}({\bf u})$. The final time $t_{\rm fin}$ is chosen
such that at time $t_{\rm fin}+\Delta t$ spurious oscillations start to
appear.

\begin{overpic}[abs,width=433.62pt,unit=1mm,scale={.25}]{Figures/x_larger_domain}
\put(1.0,54.0){(a)} \end{overpic}

\begin{overpic}[abs,width=433.62pt,unit=1mm,scale={.25}]{Figures/x_larger_domain_zoom}
\put(1.0,54.0){(b)} \end{overpic}

Figure 6: In this specific example, we consider a larger domain, with
$x\in[0,1.5]$. Since the function $R(x)=\int S(x)dx>0$, we choose another
expression for the function $S(x)$ in (4.3), with $m=-1.98/1.5$ and $q=1$.
Here we show the solution D for different times (a) and a zoom-in (b) in
x-direction, to show in which part of the domain the solution is closer to
zero. The initial condition is D(t=0) = 0.1, $\Delta t=10^{-6}$ and
$h=10^{-4}$.

###  4.2 Time discretization for the $\delta$-system

In this section we discretize in time the expression in (3.4), where the
choice for the sink/source function is defined as a sum of delta functions, in
(3). Here we rewrite (3.4) as

$\displaystyle\displaystyle\frac{\partial D_{1}}{\partial t}=$
$\displaystyle\frac{1}{D_{1}^{2}}a^{2}\,\Phi^{\prime\prime}(u(x_{0},t))$
$\displaystyle x\in[x_{0},x_{1}),\,t\in[0,t_{\rm fin}]$ (4.15)
$\displaystyle\frac{\partial D_{2}}{\partial t}=$
$\displaystyle\frac{1}{D_{2}^{2}}(a-b)\left[a\Phi^{\prime\prime}(u(x_{0},t))-b\Phi^{\prime\prime}(u(x_{1},t))\right]$
$\displaystyle x\in[x_{1},1],\,t\in[0,t_{\rm fin}].$ (4.16)

where

$\displaystyle u(x_{0},t)$
$\displaystyle=\frac{a(x_{1}-x_{0})}{D_{1}(t)}+\frac{(a-b)(1-x_{1})}{D_{2}(t)},$
(4.17) $\displaystyle u(x_{1},t)$
$\displaystyle=\frac{(a-b)(1-x_{1})}{D_{2}(t)}.$ (4.18)

To fix ideas, we focus on the case $a>b$. The expression for the entropy we
are interested in, is the following
$\Phi^{\prime\prime}(u)=(u^{2}(\sin(u)+\sigma)+1)^{-1},$ where $\sigma$ is
selected such that

$\displaystyle\frac{\sigma-1}{\sigma+1}<\frac{b}{a},$

ensuring that we do not fall into the case (3.22) addressed by the existing
theory.

Equations (4.15–4.16) are two non-linear ODE’s, and to apply IMEX schemes as
before, we multiply and divide by the same quantities the two equations:
$D_{1}$ in the first equation and $D_{2}$ in the second one, obtaining

$\displaystyle\displaystyle\frac{\partial D_{1}}{\partial t}=$ $\displaystyle
M_{1}(D_{1},D_{2})\,D_{1}\qquad x\in[x_{0},x_{1}),$ (4.19)
$\displaystyle\frac{\partial D_{2}}{\partial t}=$ $\displaystyle
M_{2}(D_{1},D_{2})\,D_{2}\qquad x\in[x_{1},1].$ (4.20)

where

$\displaystyle M_{1}(D_{1},D_{2})$
$\displaystyle=\frac{1}{D_{1}^{3}}a^{2}\,\Phi^{\prime\prime}(u(x_{0},t))$
(4.21) $\displaystyle M_{2}(D_{1},D_{2})$
$\displaystyle=\frac{1}{D_{2}^{3}}(a-b)\left[a\Phi^{\prime\prime}(u(x_{0},t))-b\Phi^{\prime\prime}(u(x_{1},t))\right]$
(4.22)

Now we apply IMEX scheme to (4.15–4.16), and the scheme reads

$\displaystyle D_{1,E}^{i}$ $\displaystyle=D_{1,E}^{n}+\Delta
t\,\sum_{j=1}^{\rm
s-1}\widetilde{a}_{i,j}M_{1}(D_{1,E}^{j},D_{2,E}^{j})D_{1,I}^{j},\quad
i=1,\cdots,s$ (4.23) $\displaystyle D_{2,E}^{i}$
$\displaystyle=D_{2,E}^{n}+\Delta t\,\sum_{j=1}^{\rm
s-1}\widetilde{a}_{i,j}M_{2}(D_{1,E}^{j},D_{2,E}^{j})D_{2,I}^{j},\quad
i=1,\cdots,s$ (4.24) $\displaystyle D_{1,I}^{i}$
$\displaystyle=D_{1,I}^{n}+\Delta t\,\sum_{j=1}^{\rm
s}\widehat{a}_{i,j}M_{1}(D_{1,E}^{j},D_{2,E}^{j})D_{1,I}^{j},\quad
i=1,\cdots,s$ (4.25) $\displaystyle D_{2,I}^{i}$
$\displaystyle=D_{2,I}^{n}+\Delta t\,\sum_{j=1}^{\rm
s}\widehat{a}_{i,j}M_{2}(D_{1,E}^{j},D_{2,E}^{j})D_{2,I}^{j},\quad
i=1,\cdots,s.$ (4.26)

The numerical solution is finally updated with

$\displaystyle D_{1}^{n+1}$ $\displaystyle=D_{1}^{n}+\Delta t\sum_{i=1}^{\rm
s}\widehat{b}_{\rm i}M_{1}(D_{1,E}^{i},D_{2,E}^{i}){D}_{1,I}^{i}$ (4.27)
$\displaystyle D_{2}^{n+1}$ $\displaystyle=D_{1}^{n}+\Delta t\sum_{i=1}^{\rm
s}\widehat{b}_{\rm i}M_{2}(D_{1,E}^{\rm i},D_{2,E}^{\rm i}){D}_{2,I}^{\rm i}.$
(4.28)

####  4.2.1 Accuracy test and Results

In this section we first prove the accuracy of the time discretization, where
the IMEX scheme used is a third order semi–implicit Runge–Kutta methods is
given by the IMEX-SSP3(4,3,3) L-stable Scheme: SSP-LDIRK3(4,3,3).

In Fig. 7 we show the third–order accuracy in time of the numerical scheme in
(4.27), on the left of the plot (in blue) for the variable $D_{1}$, and on the
right of the plot (in orange) for the variable $D_{2}$. The initial condition
chosen is $D^{0}_{1,E}=D^{0}_{2,E}=1$, the expression for the second
derivative of the entropy function is $\Phi^{\prime\prime}(u)=u^{2}$ and the
other parameters are : $x_{0}=0.1,x_{1}=0.9,a=1,b=0.9$.

In Fig. 8 we show the time evolution of the variables $D_{1},D_{2}$ in (a), a
zoom-in in time of panel (a) in (b) and of $\Phi^{\prime\prime}(u(x_{1},t))$
in (c), where the final time $t_{\rm fin}$ is chosen such that at time $t_{\rm
fin}+\Delta t$ spurious oscillations appear. However, we are able to show the
change of sign of the variable $D_{2}$: at time $t=t_{\rm fin}-\Delta t$ its
value is positive, while at time $t=t_{\rm fin}$ it is negative. Another
agreement is showed in panel (d), where we see the minimum of the variable
$D_{2}$ at final time $t_{\rm fin}$ (chosen in such a way that at time $t_{\rm
fin}+\Delta t$ the solution $D_{2}$ becomes negative). We see that in panel
(d) $D_{2}$ goes closer to zero when we refine the time step, in alignment
with the expectations as $\Delta t\to 0$. The initial condition in these tests
is $D_{1}(t=0)=D_{2}(t=0)=0.1,$ and the other parameters in (4.15–4.18) are:
$x_{0}=0.1,x_{1}=0.9,a=1,b=0.99$ and $\Delta t=10^{-6}$. The expression for
the second derivative of the entropy function is
$\Phi^{\prime\prime}(u)=(u^{2}(\sin(u)+\sigma)+1)^{-1},$ with $\sigma=2$.

\begin{overpic}[abs,width=173.44534pt,unit=1mm,scale={.25}]{Figures/error_plot_D1D2_433_quadrat}
\end{overpic} Figure 7: Accuracy plot of the numerical scheme defined in
(4.27), where the Butcher Tableau are defined in SSP-LDIRK3(4,3,3) scheme.
Here we have in the same plot, (on the left part, in blue) the error of the
$D_{1}$ variable, and (on the right part, in orange) the error of the $D_{2}$
variable. The expression for the error is the same as the one in (4.14), with
$\Delta t_{\rm ref}=10^{-5}$. In this tests the expression for the second
derivative of the entropy function is $\Phi^{\prime\prime}(u)=u^{2}$. The
other parameters are: $x_{0}=0.1,x_{1}=0.9,a=1,b=0.9$.

\begin{overpic}[abs,width=390.25534pt,unit=1mm,scale={.25}]{Figures/D1_D2_time}
\put(1.0,49.0){(a)} \end{overpic}

\begin{overpic}[abs,width=390.25534pt,unit=1mm,scale={.25}]{Figures/D2_time_zoom}
\put(1.0,49.0){(b)} \end{overpic}

\begin{overpic}[abs,width=390.25534pt,unit=1mm,scale={.25}]{Figures/phi2_time_zoom}
\put(1.0,49.0){(c)} \end{overpic}

\begin{overpic}[abs,width=390.25534pt,unit=1mm,scale={.25}]{Figures/minD}
\put(1.0,49.0){(d)} \end{overpic}

Figure 8: Time evolution of the functions $D_{1},D_{2}$ in (a), zoom of
$D_{2}$ in (b) and of $\Phi^{\prime\prime}(u(x_{1},t))$ in (c). Initial
conditions ${D_{1}}(t=0)=D_{2}(t=0)=0.1$ and the final time $t_{\rm fin}$ is
chosen such that at time $t_{\rm fin}+\Delta t$ spurious oscillations start to
appear. In this test the expression for the second derivative of the entropy
function is $\Phi^{\prime\prime}(u)=(u^{2}(\sin(u)+2)+1)^{-1}$, with
$\sigma=2$. The other parameters are: $x_{0}=0.1,x_{1}=0.9,a=1,b=0.99$ and
$\Delta t=10^{-6}$. In panel (d), by varying the value of the time step, we
display the minimum of $D_{2}$ at the final time $t_{\rm fin}$, such that at
time $t_{\rm fin}+\Delta t$ the solution $D_{2}$ becomes negative.

##  5 Conclusions

In this paper we provided a proof of the well-posedness of self-regulating
processes that model biological transportation networks, as presented in [15].
Such models arise from a constrained minimization problem of an entropy
dissipation coupled to the conservation law for a quantity representing the
concentration of a chemical species, ions, nutrients or material pressure. Our
focus was particularly directed toward the 1D setting, where we specifically
investigate Dirichlet and Neumann boundary conditions. The analysis leads to
either a local existence and uniqueness result if $D$ reaches zero in finite
time or, alternatively, a global one. We explored several conditions under
which global existence is guaranteed.

Moreover, to gain a deeper understanding of the behavior of $D$ and its
influence on the system under investigation, we systematically studied a range
of possible scenarios. This entails analyzing the system with a signed measure
distribution of sources and sinks.

Finally, we performed several numerical tests in which the solution $D$
approaches zero. These tests validate and provide further evidence of the
local existence in specific cases, as hinted at in previous considerations.

As implied by the title, a natural extension of this work would be the
development of a theory in 2D, that is already under investigation.

## References

  * [1] A. Arnold, P. Markowich, and G. Toscani,  _On large time asymptotics for drift-diffusion-poisson systems_ ,  Transport Theory and Statistical Physics, Vol. 29, pp. 571–581, 2000.
  * [2] C. Astuto, D. Boffi, J. Haskovec, P. Markowich, and G. Russo,  _Comparison of Two Aspects of a PDE Model for Biological Network Formation_ ,  _Mathematical and Computational Applications_ , (2022), DOI: 10.3390/mca27050087
  * [3] C. Astuto, D. Boffi, J. Haskovec, P. Markowich, and G. Russo, _Asymmetry and condition number of an elliptic-parabolic system for biological network formation_ ,  Communications on Applied Mathematics and Computation, (2023), DOI: 10.1007/s42967-023-00297-3
  * [4] C. Astuto, D. Boffi, and F. Credali,  _Finite element discretization of a biological network formation system: a preliminary study_ , (submitted), DOI: arXiv:2303.10625
  * [5] S. Avgerinos, F. Bernard, A. Iollo, and G. Russo,  _Linearly implicit all Mach number shock capturing schemes for the Euler equations_ ,  Journal of Computational Physics, (2019).
  * [6] M. Barthélemy,  _Spatial networks_ ,  Physics Reports, 499 (2011), 1-101.
  * [7] M. Bernot, V. Caselles, J.-M. Morel,  _Optimal Transportation Networks: Models and Theory_ ,  LNM 1955, Springer-Verlag Berlin Heidelberg, 2009.
  * [8] S. Bohn, B. Andreotti, S. Douady, J. Munzinger, and Y. Couder,  _Constitutive property of the local organization of leaf venation networks_ ,  Physical Review E, 65 (2002).
  * [9] S. Boscarino, F. Filbet and G. Russo,  _High Order Semi–implicit Schemes for Time Dependent Partial Differential Equations_ ,  Journal of Scientific Computing, 2006
  * [10] W. Boscheri, and L. Pareschi,  _High order pressure-based semi–implicit IMEX schemes for the 3D Navier-Stokes equations at all Mach numbers_ ,  Journal of Computational Physics, (2021).
  * [11] S. Boscarino  _High-Order Semi–implicit Schemes for Evolutionary Partial Differential Equations with Higher Order Derivatives_ ,  Journal of Scientific Computing, 2023.
  * [12] F. Corson,  _Fluctuations and Redundancy in Optimal Transport Networks_ ,  _Physical Review Letters_ , 104 (2010), 048703.
  * [13] D. Gilbarg and N. Trudinger,  _Elliptic partial differential equations of second order_ ,  Springer, Vol. 224, 1977.
  * [14] E. Hairer, S.P. Norsett, and G. Wanner,  _Solving Ordinary Differential Equation. I. Non-stiff Problems_ , Springer Series in Computational Mathematics, 8 Springer, (1993).
  * [15] J. Haskovec, P. Markowich and S. Portaro,  _Emergence of biological transportation networks as self-regulated process,_ Discrete and Continuous Dynamical Systems, Vol. 43(3 & 4), pp. 1499–1515, 2023, DOI: 10.3934/dcds.2022159
  * [16] E. Katifori, G. J. Szöllosi and M. O. Magnasco,  _Damage and Fluctuations Induce Loops in Optimal Transport Networks_ ,  _Physical Review Letters_ , 104 (2010), 048704.
  * [17] L. Pareschi, and G. Russo,  _Implicit-explicit Runge-Kutta schemes for stiff systems_ ,  Recent trends in numerical analysis, 2000
  * [18] L. Pareschi, and G. Russo,  _Implicit–explicit Runge–Kutta schemes and applications to hyperbolic systems with relaxation_ ,  Journal of Scientific computing, 2005, Springer.
  * [19] A. Tero, S. Takagi, T. Saigusa, K. Ito, D. Bebber, M. Fricker, K. Yumiki, R. Kobayashi and T. Nakagaki,  _Rules for Biologically Inspired Adaptive Network Design_ ,  Science, 327 (2010), 439-442.


# Functional Equations for the Stochastic Exponential

B. Chikvinidze1), M. Mania2) and R. Tevzadze3)

( )

Abstract. We consider stochastic versions of the Cauchy exponential functional
equation and give a martingale characterization of the general solution.

1) Georgian-American University and Institute of Cybernetics of Georgian
Technical Univercity, Tbilisi, Georgia
(e-mail<EMAIL_ADDRESS>
2) Razmadze Mathematical Institute of Tbilisi State University and Georgian-
American University, Tbilisi, Georgia,
(e-mail<EMAIL_ADDRESS>
3) Georgian-American University and Institute of Cybernetics of Georgian
Technical Univercity, Tbilisi, Georgia
(e-mail<EMAIL_ADDRESS>
2010 Mathematics Subject Classification. 60G44, 60J65, 97I70

Keywords: Stochastic Exponential, Martingales, Functional Equations, Non-
anticipative functionals, Brownian Motion.

##  1 Introduction

For a continuous semimartingale $X=(X_{t},t\geq 0)$, with $X_{0}=0$, the
stochastic exponential is defined as

${\cal E}_{t}(X)=e^{X_{t}-\frac{1}{2}\langle X\rangle_{t}},\;\;\;t\geq 0,$ (1)

where $\langle X\rangle$ is the square characteristic of the martingale part
of $X$.

${\cal E}_{t}(X)$ is the unique solution of the linear stochastic differential
equation

$Z_{t}=1+\int_{0}^{t}Z_{s}dX_{s},\;\;\;t\geq 0$

and in this sense $Z_{t}={\cal E}_{t}(X)$ is the stochastic analogue of the
usual exponential function $f(x)=e^{cx}$, which is the unique solution of the
linear differential equation $f_{x}(x)=cf(x),f(0)=1.$

On the other hand, it is well known that the property

$e^{c(x+y)}=e^{cx}e^{cy},\;\;\;\text{for all}\;\;x,y\in R$ (2)

is a characterizing property of exponential functions and in a very wide class
of functions (e.g., in the class of measurable functions) $f(x)=e^{cx}$, where
$c\in R$ is some constant, is the general solution of the Cauchy exponential
functional equation

$f(x+y)=f(x)f(y)\;\;\;\text{for all}\;\;x,y\in R.$ (3)

Similar to (2) property of the stochastic exponential (4) is the equality

${\cal E}_{t}(X){\cal E}_{t}(Y)={\cal E}_{t}(X+Y+\langle
X,Y\rangle),\;\;\;t\geq 0,$ (4)

which is valid for any pair $X,Y$ of semimartingales, where the additional
term $\langle X,Y\rangle$ is the mutual characteristic of semimartingales $X$
and $Y$. For the properties of stochastic exponential and for all unexplained
notations from the martingale theory, we refer to [6] or [10] .

By definition, the stochastic exponential ${\cal E}_{t}(X)$ of a continuous
semimartingale $X$ can be expressed as a two-dimensional exponential function
of a semimartingale and its square characteristic at time $t$

${\cal E}_{t}(X)=f(\langle X\rangle_{t},X_{t}),\;\;\;t\geq 0,$ (5)

where $f(u,v)=e^{v-\frac{1}{2}u},\;\;\;u\geq 0,v\in R.$ It follows from (4),
that this function satisfies the following functional equation

$f(\langle X\rangle_{t},X_{t})f(\langle Y\rangle_{t},Y_{t})=f(\langle
X+Y\rangle_{t},X_{t}+Y_{t}+\langle X,Y\rangle_{t}),$ (6)

valid for any continuous semimartingales $X$ and $Y$ vanishing at $0$. It is
easy to verify that (6) is also true for the function

$f(u,v)=e^{cv-\frac{c}{2}u},\;\;\;u\geq 0,v\in R,$ (7)

for any constant $c\in R.$

Denote by $\cal S$ the class of continuous semimartingales and let ${\cal V}$
be a subclass of $\cal S$, for which equation (6) is satisfied. Following [9]
we call the class ${\cal V}$ the domain of validity of equation (6).

One is led to the following question: given a class of two-dimensional
functions (e.g., measurable, continuous), how small a class of semimartingales
$\cal V$ can we take, so that the solutions of (6) remain to be only functions
of the form (7).

We show in Theorem 1, that if equation (6) is satisfied for the class of
stochastic integrals $h\cdot W$ with respect to the given Brownian Motion $W$
and with deterministic integrands $0\leq h\leq 1$, then the solution of
equation (6) is of the form (7). In section 2 we prove this theorem by using
the Cauchy exponential functional equation and the corresponding ”almost”
version. In section 3 we give a martingale characterization of equation (6),
which gives also a probabilistic proof of this theorem.

If the domain of validity of equation (7) consists only of a single Brownian
Motion $W$ (${\cal V}=\\{W\\}$), which corresponds (if we take $X=Y=W$ in (6))
to the equation

$f^{2}(t,W_{t})=f(4t,2W_{t}+t),\;\;\;t\geq 0,$ (8)

then there exists a different from (7) continuous solution of equation (6),
although any analytic solution of (8) is of the form (7). See Proposition 1
and a counterexample before this proposition.

Note that the stochastic exponential ${\cal E}_{t}(X)$ transforms the class of
continuous semimartingales into itself and this mapping is non-anticipative in
the following sense: The mapping $F:R_{+}\times{\cal S}\rightarrow R$ is non-
anticipative, if for any continuous semimartingales $X$ and $Y$ and $t\geq 0$

$F(t,X)=F(t,Y),\;\;\;\text{when}\;\;\;X_{s}=Y_{s}\;\;\;\text{for
all}\;\;\;s\leq t.$

Therefore, it seems natural to consider a functional equation for stochastic
exponents in terms of non-anticipative functionals

$F(t,X)F(t,Y)=F(t,X+Y+\langle X,Y\rangle),$ (9)

where by the property (4) the stochastic exponential ${\cal E}_{t}(X)$
satisfies this equation and it represents more general form of equation (6).
But (9) will be not a characterizing property of the stochastic exponent,
since there exists a whole class of solutions of (9) which are not stochastic
exponentials. E.g., if

$F(t,X)=e^{\int_{0}^{t}{\cal K}(t,s)d(X_{s}-\frac{1}{2}\langle
X\rangle_{s})},\;\;\;X\in\cal S,$ (10)

where $(k(t,s),s\geq 0,t\geq 0)$ is bounded, measurable deterministic
function, then $F(t,X)$ defined by (10) satisfies (9), but such processes are
not always stochastic exponentials, see Theorem 3 of section 4 and Remark 2 at
the end of the paper. In Theorem 3 we prove that under some restriction on the
class of non-anticipative functionals, the general solution of (9) is of the
form (10).

##  2 Functional equation for a function of a semimartingale and its square
characteristic

Let $W=(W_{t},t\geq 0)$ be a standard Brownian Motion defined on a complete
probability space $(\Omega,\cal F,P)$. Let $F=({{\cal F}}_{t},t\geq 0)$ be a
filtration satisfying the usual conditions of right-continuity and
completeness. Let $F^{W}=({{\cal F}}^{W}_{t},t\geq 0)$ be the filtration
generated by the Brownian Motion $W$.

Let $\cal S$ (resp. $\cal M$) be the class of continuous semimartingales
(martingales) vanishing at $0$.

Denote by ${\cal M}^{W}$ the class of continuous local martingales
$M=(M_{t},t\geq 0)$ adapted to the filtration $F^{W}$ with $M_{0}=0$, i.e.,
this is the class of stochastic integrals with respect to the Brownian Motion
$W$.

Let ${\cal M}^{W}(I)$ be a sub-class of stochastic integrals $h\cdot W$ with
respect to the Brownian Motion $W$ with integrands $h$, such that
$h_{u}=I_{[s\leq u\leq t]},0\leq s\leq t$.

Let $\left(f\left(u,v\right),u\geq 0,v\in R\right)$ be a function of two
variables . We consider the functional equation

$f(\langle X\rangle_{t},X_{t})f(\langle Y\rangle_{t},Y_{t})=f(\langle
X+Y\rangle_{t},X_{t}+Y_{t}+\langle X,Y\rangle_{t}),$ (11)

for any $X,Y\in{\cal V}$ and $P$\- a.e. for each $t\geq 0$, where $\cal V$ is
some class of continuous semimartingales, the domain of validity of equation
(11).

Theorem 1. Let $\left(f\left(u,v\right),u\geq 0,v\in R\right)$ be a function
of two variables. Then the following assertions are equivalent:

$a)$ The function $f$ is a continuous strictly positive solution of the
functional equation (11) with the domain of validity $\cal V=\cal S$.

$b)$ The function $f$ is a continuous strictly positive solution of the
functional equation (11) with the domain of validity ${\cal V}={\cal
M}^{W}(I)$.

$c)$ The function $f$ is of the form

$f(u,v)=e^{cv-\frac{c}{2}u}\;\;\;\;\text{for some constant}\;\;\;c\in R.$

If we shall consider only measurable solutions, then the following two
conditions will be equivalent:

$b^{\prime})$ The function $f$ is a measurable strictly positive solution of
the functional equation (11) with the domain of validity ${\cal V}={\cal
M}^{W}(I)$.

$c^{\prime})$ The function $f=(f(u,v),u\geq 0,v\in R)$ coincides with the
function

$e^{cv-\frac{c}{2}u}\;\;\;\;\text{for some constant}\;\;\;c\in R$ (12)

almost everywhere with respect to the Lebesgue measure on $R_{+}\times R$.

###### Proof.

The implication $a)\to b)$ is evident. Let us prove the implication $b)\to
c)$.

It follows from equation (11) that for any bounded deterministic functions $h$
and $g$

$f\Big{(}\int_{0}^{t}h^{2}_{u}du,\int_{0}^{t}h_{u}dW_{u}\Big{)}f\Big{(}\int_{0}^{t}g^{2}_{u}du,\int_{0}^{t}g_{u}dW_{u}\Big{)}=$
(13)
$=f\Big{(}\int_{0}^{t}(g_{u}+h_{u})^{2}du,\int_{0}^{t}(g_{u}+h_{u})dW_{u}+\int_{0}^{t}g_{u}h_{u}du\Big{)}$

$P$\- a.e. for each $t\geq 0$.

For any fixed pair $s\leq t$ if we take $h_{u}=I_{(u<s)}$ and $g_{u}=I_{s\leq
u\leq t)}$, from (13) we obtain that $P$-a.s.

$f(s,W_{s})f(t-s,W_{t}-W_{s})=f(t,W_{t}).$ (14)

From (14) we have that

$0=EI_{(f(s,W_{s})f(t-s,W_{t}-W_{s})\neq f(t,W_{t}))}=$ (15)
$=\int_{R}\int_{R}I_{(f(s,x)f(t-s,y)\neq
f(t,x+y))}\rho(s,x)\rho(t-s,y-x)dxdy,$

where $\rho(s,x)=\frac{1}{\sqrt{2\pi s}}e^{\frac{x^{2}}{2s}}$.

Therefore, for any $s>0,t>0,s\leq t$

$f(s,x)f(t-s,y)=f(t,x+y)\;\;\;\text{for all}\;\;\;x,y\in R$ (16)

almost surely with respect to the Lebesgue measure on $R^{2}$.

If $f(t,x)$ is continuous, then (16) is fulfilled for all $s>0,t>0,s\leq
t,x\in R,y\in R$ and $f$ satisfies the two-dimensional Cauchy exponential
functional equation. The general continuous solution of this equation is of
the form (see, e.g., [1])

$f(t,x)=\exp\\{cx+bt\\}\;\;\;\text{for some}\;\;\;b,c\in R.$ (17)

On the other hand inserting $X=W$ and $Y=W$ in (11) we have

$f^{2}(t,W_{t})=f(4t,2W_{t}+t)\;\;\;t\geq 0$ (18)

$P$-a.s. and substituting $f(t,x)=\exp\\{cx+bt\\}$ in (18) we obtain that

$\exp\\{2cW_{t}+2bt\\}=\exp\\{2cW_{t}+ct+4bt\\},$

which implies that $b=c/2$, hence $f$ is of the form
$\exp\\{cx-\frac{c}{2}t\\}$.

$c)\to a)$ is evident. Indeed, if $f(u,v)=\exp\\{cv-\frac{c}{2}u\\}$ for some
$c\in R$, then for any $X,Y\in\cal S$

$f(\langle X\rangle_{t},X_{t})f(\langle
Y\rangle_{t},Y_{t})=\exp\\{cX_{t}-\frac{c}{2}\langle
X\rangle_{t}\\}\exp\\{cY_{t}-\frac{c}{2}\langle Y\rangle_{t}\\}=$
$=\exp\\{c(X_{t}+Y_{t})-\frac{c}{2}\langle X+Y\rangle_{t}+c\langle
X,Y\rangle_{t}\\}=f(\langle X+Y\rangle_{t},X_{t}+Y_{t}+\langle
X,Y\rangle_{t})$

$P$-a.s. for any $t\geq 0$.

The proof of the second part of Theorem 1 is similar if we use corresponding
results on ′′almost′′ solutions of equation (16).

$b^{\prime})\to c^{\prime})$. It follows from results of [3] and [7] that (16)
implies

$f(t,x)=\exp\\{cx+bt\\}\;\;\;\text{for some}\;\;\;b,c\in R,$ (19)

almost everywhere with respect to the Lebesgue measure on $R_{+}\times R$.

From (18) we have that

$f^{2}(t,x)=f(4t,2x+t)\;\;\;t\geq 0$ (20)

almost everywhere with respect to the Lebesgue measure on $R_{+}\times R$.
Therefore, from (19) and (20), avoiding three null sets, we obtain that

$\exp\\{2cx+2bt\\}=\exp\\{2cx+ct+4bt\\},$

almost everywhere with respect to the Lebesgue measure on $R_{+}\times R$,
which implies that $b=c/2$.

$c^{\prime})\to b^{\prime})$. Since for any $M\in{\cal M}^{W}(I)$ the random
variable $M_{t}$ admits Gaussian distribution, it follows from $c^{\prime})$
that

$f(\langle M\rangle_{t},M_{t})=\exp\\{cM_{t}-\frac{c}{2}\langle
M\rangle_{t}\\}$ (21)

$P$-a.s. for any $t\geq 0$, which will imply this implication, similarly to
the proof of corresponding implication in first part of Theorem 1.

∎

Remark. Note that, equality (21) can not be deduced from condition
$c^{\prime})$ for a general continuous semimartingale $X$.

If the domain of validity of equation (11) consists only of a single Brownian
Motion $W$ (${\cal V}=\\{W\\}$), which corresponds to the equation

$f^{2}(t,W_{t})=f(4t,2W_{t}+t),\;\;\;t\geq 0,P-a.s.$ (22)

then there exists a different from (7) continuous solution of equation (22).
Let $g(t,x)=\ln f(t,x)$. Then (22) is equivalent to equation

$g(t,W_{t})=\frac{1}{2}g(4t,2W_{t}+t),\;\;\;t\geq 0.$ (23)

Let $g(t,x)=|(2x-t)t|^{\frac{1}{3}}$. It is evident, that this function
satisfies equation (23), since

$g(4t,2x+t)=|(4x-2t)4t|^{\frac{1}{3}}=2|(2x-t)t|^{\frac{1}{3}}=2g(t,x).$

Therefore, $f(t,x)=\exp{|(2x-t)t|^{\frac{1}{3}}}$ is a continuous solution of
equation (22) different from (7).

′′ It has been shown in various works that when some additional smoothness
assumptions are imposed on set of functions of classical functional equations,
then even if the domain of validity is quite small, the set of solution does
not grow ′′ (see [12] for this comment and the references therein). In the
following proposition we show that if the domain of validity contains only a
single Brownian Motion ($\cal V=W$), then any analytic solution of (22) is of
the form (7), i.e., the set of solutions is the same as for $\cal V=\cal S$.

Proposition 1. If $f$ is analytical function satisfying (22), then $f$ is of
the form (7).

###### Proof.

It follows from (20) and continuity of $f(t,x)$, that $f$ satisfies equation

$f^{2}(t,x)=f(4t,2x+t)\;\;\;\text{for all}\;\;\;t\geq 0,x\in R.$ (24)

Taking logarithms of both sides of (24) we obtain functional equation

$h(t,x)=\frac{1}{2}h(4t,2x+t)\;\;\;\text{for all}\;\;\;t\geq 0,x\in R,$ (25)

for $h(t,x)=\ln f(t,x)$.

It is evident that $h(0,0)=0$. Differentiating equality (25) at $t$ we have

$h_{t}(t,x)=2h_{t}(4t,2x+t)+\frac{1}{2}h_{x}(4t,2x+t).$ (26)

This implies that $h_{t}(0,0)=2h_{t}(0,0)+h_{x}(0,0)$ and hence

$h_{t}(0,0)=-\frac{1}{2}h_{x}(0,0).$ (27)

It easy to see that all higher order derivatives of $h$ at $t=0,x=0$ are equal
to zero.

Indeed, it follows from (25) that for any $m\geq 1$ and $n\geq 1$

$\frac{\partial^{m+n}}{\partial t^{m}\partial
x^{n}}h(t,x)=\frac{\partial^{m+n-1}}{\partial t^{m}\partial
x^{n-1}}h_{x}(t,x)=\frac{\partial^{m+n-1}}{\partial t^{m}\partial
x^{n-1}}h_{x}(4t,2x+t)=$ $=\frac{1}{2}\frac{\partial^{m+n}}{\partial
t^{m}\partial
x^{n}}\big{(}h(4t,2x+t)\big{)}_{x}=\frac{1}{2}\frac{\partial^{m+n}}{\partial
t^{m}\partial x^{n}}h(4t,2x+t),$

which implies that

$\frac{\partial^{m+n}}{\partial t^{m}\partial
x^{n}}h(0,0)=\frac{1}{2}\frac{\partial^{m+n}}{\partial t^{m}\partial
x^{n}}h(0,0)\;\;\;\text{and}\;\;\;\frac{\partial^{m+n}}{\partial t^{m}\partial
x^{n}}h(0,0)=0.$

If $n=0$ then it follows from (25) and (26) that

$\frac{\partial^{m}}{\partial t^{m}}h(t,x)=\frac{\partial^{m-1}}{\partial
t^{m-1}}h_{t}(t,x)=\frac{\partial^{m-1}}{\partial
t^{m-1}}\big{(}2h_{t}(4t,2x+t)+\frac{1}{2}h_{x}(4t,2x+t)\big{)}=$
$=\frac{1}{2}\frac{\partial^{m}}{\partial
t^{m}}h(4t,2x+t)+\frac{1}{4}\frac{\partial^{m}}{\partial t^{m-1}\partial
x}h(4t,2x+t).$

Since $\frac{\partial^{m+n}}{\partial t^{m}\partial x^{n}}h(0,0)=0$ for any
$n\geq 1$, this implies that

$\frac{\partial^{m}}{\partial
t^{m}}h(0,0)=\frac{1}{2}\frac{\partial^{m}}{\partial t^{m}}h(0,0)$

and this derivative at point $(0,0)$ is also equal to zero.

Therefore, since $f$ is analytic, we obtain from (27) that

$h(t,x)=h_{t}(0,0)t+f_{x}(0,0)x=h_{x}(0,0)x-\frac{h_{x}(0,0)}{2}t,$

which means that $h$ is of the form $h(t,x)=cx-\frac{c}{2}t$ for some $c\in
R$. ∎

##  3 A martingale approach

In this section we give a martingale characterization of equation (11), which
gives also a probabilistic proof of Theorem 1.

Let first show that $f(t,W_{t})$ is integrable at any power.

Lemma 1. If $\left(f\left(u,v\right),u\geq 0,v\in R\right)$ is a measurable
strictly positive solution of the functional equation (11), then

$E|f(t,W_{t})|^{p}<\infty$

for any $t\geq 0$ and $p\in R$.

###### Proof.

To show that $f\left(t,W_{t}\right)$ is integrable for any $t\geq 0$ we shall
use the idea from [13] on application of the Bernstein theorem.

Substituting $s=t/2$ in (14) and taking logarithms we have

$\ln f(t/2,W_{t/2})+\ln f(t/2,W_{t}-W_{t/2})=\ln f(t,W_{t})\;\;\;\;P-a.s.$
(28)

Let

$X=\ln f(t/2,W_{t/2})\;\;\;\text{and}\;\;\;Y=\ln f(t/2,W_{t}-W_{t/2}).$

Then from (28)

$X+Y=\ln f(t,W_{t})\;\;\;\;P-a.s.$ (29)

and it is easy to see that

$X-Y=\ln f(t/2,W_{t/2})-\ln f(t/2,W_{t}-W_{t/2})=\ln
f(t,2W_{t/2}-W_{t})\;\;\;\;P-a.s.$

Since $E(2W_{t/2}-W_{t})W_{t}=0$, the random variables $2W_{t/2}-W_{t}$ and
$W_{t}$ are independent and, hence, the random variables $X+Y$ and $X-Y$ will
also be independent. Therefore, Bernstein’s theorem ([2], see also [11])
implies that $X$ and $Y$ are distributed normally. Therefore, $\ln f(t,W_{t})$
is also normally distributed. Hence $f(t,W_{t})$ is integrable at any power

$E|f(t,W_{t})|^{p}<\infty,$

as a random variable having log-normal distribution. ∎

Proposition 2. If $\left(f\left(u,v\right),u\geq 0,v\in R\right)$ is a
measurable strictly positive solution of the functional equation (11), then

$i)$ the processes

$N_{t}(1)=\frac{f(t,W_{t})}{Ef(t,W_{t})},\;\;\;\text{and}\;\;\;N_{t}(2)=\frac{f^{2}(t,W_{t})}{Ef^{2}(t,W_{t})},\;\;t\geq
0,$

are strictly positive martingales.

$ii)$ There exist constants $\lambda_{1}\in R$ and $\lambda_{2}\in R$ such
that the processes

$\frac{f(t,W_{t})}{e^{\lambda_{1}t}},\;\;\;\frac{f^{2}(t,W_{t})}{e^{\lambda_{2}t}}\;\;\;\text{and}\;\;\;\frac{f(4t,2W_{t}+t)}{e^{\lambda_{2}t}}\;\;\;t\geq
0,$

are strictly positive martingales. Besides

$\lambda_{1}\leq\frac{1}{2}\lambda_{2}.$ (30)

###### Proof.

$i)$ Taking expectations in equality (14), since the random variables
$W_{t}-W_{s}$ and $W_{s}$ are independent, we have that

$Ef(t-s,W_{t}-W_{s})=\frac{Ef(t,W_{t})}{Ef(s,W_{s})}.$ (31)

If we take the conditional expectations in the same equality (14), having in
mind that $W_{t}-W_{s}$ is independent of ${\cal F}_{s}^{W}$, we obtain

$E(f(t,W_{t})|{{\cal
F}_{s}^{W}})=f(s,W_{s})Ef(t-s,W_{t}-W_{s})\;\;\;\;P-a.s..$ (32)

Therefore, substituting the expression of $Ef(t-s,W_{t}-W_{s})$ from (31) into
(32) we get the martingale equality

$E\Big{(}\frac{f(t,W_{t})}{Ef(t,W_{t})}\Big{|}{\cal
F}_{s}\Big{)}=\frac{f(s,W_{s})}{Ef(s,W_{s})},\;\;\;\;P-\text{a.s.}.$

The martingale property of the process $N_{t}(2)$ is proved similarly, if we
use equality

$f^{2}(s,W_{s})f^{2}(t-s,W_{t}-W_{s})=f^{2}(t,W_{t})\;\;\;\;P-a.s.$ (33)

instead of equality (14).

$ii)$ Let $g(t)\equiv Ef(t,W_{t})$. Then $Ef(t-s,W_{t}-W_{s})=g(t-s)$ and it
follows from (31) that $(g(t),t\geq 0)$ satisfies the Cauchy exponential
functional equation

$g(s)g(t-s)=g(t),\;\;\;s\geq 0,t\geq 0,s\leq t$

a general measurable solution of which is of the form $e^{\lambda t}$ for some
$\lambda\in R$ (see, e.g., [1]). Therefore the proof of $ii)$ follows from
assertion $i)$.

Since $f^{2}(t,W_{t})e^{-\lambda_{2}t}$ is a martingale, the equality (18)
implies that the process $f(4t,2W_{t}+t)e^{-\lambda_{2}t}$ is also a
martingale.

The inequality (30) follows from the martingale property of $N(1)$ and $N(2)$
and from the Hôlder inequality, since

$1=E\frac{f(t,W_{t})}{e^{\lambda_{1}t}}\leq
E^{1/2}\frac{f^{2}(t,W_{t})}{e^{2\lambda_{1}t}}=$
$=E^{1/2}\frac{f^{2}(t,W_{t})}{e^{\lambda_{2}t}}\cdot\frac{e^{\lambda_{2}t}}{e^{2\lambda_{1}t}}=e^{\big{(}\frac{\lambda_{2}}{2}-\lambda_{1}\big{)}t}.$

∎

Theorem 2. Each assertion $a),b),c)$ of Theorem 1 is equivalent to

$d)$ $\left(f\left(u,v\right),u\geq 0,v\in R\right)$ is a continuous strictly
positive function such that $f(0,0)=1$ and the processes

$\frac{f(t,W_{t})}{e^{\lambda_{1}t}},\;\;\;\frac{f^{2}(t,W_{t})}{e^{\lambda_{2}t}},\;\;\;\text{and}\;\;\;\frac{f(4t,2W_{t}+t)}{e^{\lambda_{2}t}},\;\;\;t\geq
0,$

are martingales for some constants $\lambda_{1}\in R$ and $\lambda_{2}\in R$.

###### Proof.

The implication $b)\to d)$ follows from Proposition 2. Let us show that if
$d)$ is satisfied then the function $f(t,x)$ is of the form (12).

Let us define two functions:

$U(t,x)=E\big{[}f(T,W_{T})e^{-\lambda_{1}T}\big{|}W_{t}=x\big{]}\;\;\text{and}\;\;V(t,x)=E\big{[}f^{2}(T,W_{T})e^{-\lambda_{2}T}\big{|}W_{t}=x\big{]}.$

Since $U(t,x)$ and $V(t,x)$ are positive, they will be of the class $C^{1.2}$
on $(0,T)\times R$ and satisfy the ”backward” heat equation (see, e.g. [5]
page 257)

$\frac{\partial Y}{\partial t}+\frac{1}{2}\frac{\partial^{2}Y}{\partial
x^{2}}=0,\;\;\;\;0<t<T,x\in R.$ (34)

Since the processes $f(t,W_{t})e^{-\lambda_{1}t}$ and
$f^{2}(t,W_{t})e^{-\lambda_{2}t}$ are martingales, by the Markov property of
$W$ we shall have that

$U(t,W_{t})={f(t,W_{t})}{e^{-\lambda_{1}t}}\;\;\;\text{and}\;\;\;V(t,W_{t})=f^{2}(t,W_{t})e^{-\lambda_{2}t}$

and hence, by continuity of the function $f$

$U(t,x)={f(t,x)}{e^{-\lambda_{1}t}}\;\;\;\text{and}\;\;\;V(t,x)={f^{2}(t,x)}{e^{-\lambda_{2}t}},\;\;\;\text{for
all}\;\;\;t\geq 0,x\in R.$

Since

$U_{t}(t,x)=f_{t}(t,x)e^{-\lambda_{1}t}-\lambda_{1}f(t,x)e^{-\lambda_{1}t};\;\;\;U_{xx}(t,x)=\frac{1}{2}f_{xx}(t,x)e^{-\lambda_{1}t}$

and $U(t,x)$ satisfies equation (34), we shall have that $f$ also belongs to
the class $C^{1.2}$ on $(0,T)\times R$ and satisfies equation

$f_{t}(t,x)+\frac{1}{2}f_{xx}(t,x)-\lambda_{1}f(t,x)=0.$ (35)

Similarly, since

$V_{t}(t,x)=2f(t,x)f_{t}(t,x)e^{-\lambda_{2}t}-\lambda_{2}f^{2}(t,x)e^{-\lambda_{2}t};$
$V_{xx}(t,x)=2(f_{x}(t,x))^{2}e^{-\lambda_{2}t}+2f(t,x)f_{xx}(t,x)e^{-\lambda_{2}t}$

and $V(t,x)$ satisfies equation (34), we obtain that the function $f$
satisfies also the following PDE:

$2f(t,x)f_{t}(t,x)e^{-\lambda_{2}t}-\lambda_{2}f^{2}(t,x)e^{-\lambda_{2}t}+$
$\frac{1}{2}\Big{[}2(f_{x}(t,x))^{2}e^{-\lambda_{2}t}+2f(t,x)f_{xx}(t,x)e^{-\lambda_{2}t}\Big{]}=0,$

which after simplifying takes the form:

$2f(t,x)f_{t}(t,x)-\lambda_{2}f^{2}(t,x)+(f_{x}(t,x))^{2}+f(t,x)f_{xx}(t,x)=0.$
(36)

From (35) we have that $f_{xx}(t,x)=2\lambda_{1}f(t,x)-2f_{t}(t,x)$ and
substituting the right side of this equality instead of $f_{xx}$ in (36) we
get that

$f^{2}_{x}(t,x)=(\lambda_{2}-2\lambda_{1})f^{2}(t,x),$ (37)

where $\lambda_{2}-2\lambda_{1}\geq 0$ according to Proposition 2.

To find the general solution of (37) we consider two cases:
I: $\;f_{x}(t,x)=\sqrt{\lambda_{2}-2\lambda_{1}}f(t,x)$. It is evident that
$f(t,x)=\phi(t)e^{\sqrt{\lambda_{2}-2\lambda_{1}}x}$, where $\phi$ is some
function of $t$. It follows from Proposition 2 that
$Ef(t,W_{t})=e^{\lambda_{1}t}$. So
$e^{\lambda_{1}t}=\phi(t)Ee^{\sqrt{\lambda_{2}-2\lambda_{1}}W_{t}}=\phi(t)e^{\frac{\lambda_{2}-2\lambda_{1}}{2}t}$,
which implies that $\phi(t)=e^{\frac{4\lambda_{1}-\lambda_{2}}{2}t}$ and
finally $f(t,x)=e^{cx+bt}$, where $c=\sqrt{\lambda_{2}-2\lambda_{1}}$ and
$b=\frac{4\lambda_{1}-\lambda_{2}}{2}$.
II: $\;f_{x}(t,x)=-\sqrt{\lambda_{2}-2\lambda_{1}}f(t,x)$. By the same manner
as in the case I, we obtain that $f(t,x)=e^{cx+bt}$, where
$c=-\sqrt{\lambda_{2}-2\lambda_{1}}$ and
$b=\frac{4\lambda_{1}-\lambda_{2}}{2}$.

After that we can use equality $Ef^{2}(t,W_{t})=Ef(4t,2W_{t}+t)$ to get
equality $b=-c/2$ and the representation $f(t,x)=e^{cx-\frac{c}{2}t}$. ∎

##  4 A functional equation for non-anticipative functionals

The mapping $h:[0,T]\times{C[0,T]}\rightarrow R$ is non-anticipative, if for
any $\omega,\omega^{\prime}\in C[0,T]$ and $t\in[0,T]$

$h(t,\omega)=h(t,\omega^{\prime}),\;\;\;\text{when}\;\;\;\omega_{s}=\omega^{\prime}_{s}\;\;\;\text{for
all}\;\;\;s\leq t.$

Consider the class of functions $F:[0,T]\times{\cal{S}}\rightarrow R$ defined
by

$\displaystyle{\cal C}=\\{F;\;F(t,X)=e^{h(t,X-\frac{1}{2}\langle
X\rangle)},X\in{\cal S},$ $\displaystyle\text{for some continuous, non-
anticipative functional}\;h(t,\omega)\\}.$

Theorem 3. The general solution of the functional equation

$F(t,X)F(t,Y)=F(t,X+Y+\langle X,Y\rangle),\;\;\;\text{for all}\;\;\;X,Y\in\cal
S,$ (38)

in the class ${\cal C}$, is of the form

$F(t,X)=e^{\int_{0}^{t}{\cal K}(t,s)d(X_{s}-\frac{1}{2}\langle
X\rangle_{s})},\;\;\;X\in\cal S,$ (39)

where $({\cal K}(t,s),s\geq 0,t\geq 0)$ is a deterministic function with
${\cal K}(t,s)=0,t\leq s$, such that

i)

${\cal K}(t,\cdot)\;\;\text{ is cadlag and has a finite variation, for
each}\;\;t\in[0,T],$

ii)

$\int_{0}^{T}\omega_{s}{\cal K}(\cdot,ds)\;\;\;\text{ is continuous, for
each}\;\;\omega\in C[0,T].$

###### Proof.

Let $F$ be a solution of equation (38) from the class $\cal C$. Then it
follows from (38) and from the definition of the class $\cal C$ that

$\displaystyle h(t,X)+h(t,Y)=\ln F(t,X+\frac{1}{2}\langle X\rangle)+\ln
F(t,Y+\frac{1}{2}\langle Y\rangle)$ $\displaystyle=\ln
F(t,X+Y+\frac{1}{2}(\langle X\rangle+\langle Y\rangle)+\langle X,Y\rangle)$
$\displaystyle=\ln F(t,X+Y+\frac{1}{2}(\langle X+Y\rangle))=h(t,X+Y)$

for any $X,Y\in{\cal S}$.

Since a deterministic function is a semimartingale if and only if it is of
finite variation and the functions of finite variations are dense in $C$, it
follows from the continuity of $h$ that

$h(t,\omega+\omega^{\prime})={h}(t,\omega)+{h}(t,\omega^{\prime})\;\;\text{for
all}\;\;\omega,\omega^{\prime}\in{C}.$ (40)

By the Riesz theorem (see, e.g. [4]) for each $t$ there exists a cadlag
function $G(t,\cdot)$ of finite variation, such that $G(t,s)=G(t,T),s\geq t$
and

$h(t,\omega)=\int_{0}^{T}\omega_{s}G(t,ds).$ (41)

Using integration by part formula we get from (41) that

$\displaystyle
h(t,X)=\int_{0}^{t}X_{s}G(t,ds)=X_{t}G(t,t)-\int_{0}^{t}G(t,s)dX_{s}$
$\displaystyle=\int_{0}^{t}(G(t,t)-G(t,s))dX_{s}=\int_{0}^{t}{\cal
K}(t,s)dX_{s},$

where ${\cal K}(s,t)=G(t,t)-G(t,s)$ is of finite variation for each $t$ and by
continuity of $h$ the condition $ii)$ is also satisfied.

Now let us show that the non-anticipative functional
$F(t,X)=e^{\int_{0}^{t}{\cal K}(t,s)d(X_{s}-\frac{1}{2}\langle X\rangle_{s})}$
satisfies equation (38). Let first show that $F$ belongs to the class $\cal
C$.

Consider the non-anticipative functional
$h(t,\omega)=\int_{0}^{T}\omega_{s}{\cal K}(\cdot,ds)$, where $\cal K$
satisfies conditions $i)-ii)$ of the theorem. By $ii)$ this functional is
continuous and the integration by part formula gives the equality

$h(t,X)=-\int_{0}^{t}X_{s}{\cal K}(t,ds)=\int_{0}^{t}{\cal
K}(t,s)dX_{s},\;X\in{\cal S}.$ (42)

Therefore, (42) implies that

$F(t,X)=e^{\int_{0}^{t}{\cal K}(t,s)d(X_{s}-\frac{1}{2}\langle
X\rangle_{s})}=e^{h(t,X-\frac{1}{2}\langle X\rangle)},\;\;\;X\in\cal S,$ (43)

which means that $F\in\cal C$.

It is evident that $F(t,X)$ defined by (39) satisfies (38), since for any
continuous semimartingales $X$ and $Y$

$F(t,X)F(t,Y)=e^{\int_{0}^{t}k(t,s)d(X_{s}+Y_{s}-\frac{1}{2}\langle
X\rangle_{s}-\frac{1}{2}\langle Y\rangle_{s})}=$
$=e^{\int_{0}^{t}k(t,s)d(X_{s}+Y_{s}+\langle X,Y\rangle_{s}-\frac{1}{2}\langle
X+Y\rangle_{s})}=F(t,X+Y+\langle X,Y\rangle).$

∎

Remark 1. If ${\cal K}(t,s)$, ${\cal K}(t,s)=0,t\leq s$ is a measurable
function, such that $\int_{0}^{T}|{\cal K}(t,s)-{\cal
K}(t^{\prime},s)|^{2}ds\leq const|t-t^{\prime}|^{2}$, then by the Kolmogorov
theorem $\int_{0}^{t}{\cal K}(t,s)dX_{s}$ has a continuous modification, for
each $X_{t}=X_{0}+\int_{0}^{t}\sigma_{s}dW_{s}+\int_{0}^{t}b_{s}ds$ with
bounded, predictable $b,\sigma$.

Remark 2. Note that $F(t,X)$ is not always a stochastic exponential. For
example, if we take ${\cal K}(t,s)=v(t)$, where $v$ is not of finite
variation, then $F(t,X)=e^{v(t)(X_{t}-\frac{1}{2}\langle X\rangle_{t})}$ will
be not a stochastic exponential for any $X\in{\cal S}$. E.g., for $X=W$ the
process $v(t)(W_{t}-\frac{1}{2}t)$ will be not a semimartingale and, hence
$F(t,W)$ will not be a stochastic exponential.

## References

  * [1] J.Aczel, Lectures on Functional Equations and their Applications, Academic Press 1966.
  * [2] S. N. Bernstein, Ob odnom svoystve, harakterizuyushhem zakon Gaussa (in Russian) [On a characteristic property of the normal law]. Tr. leningr.Polytech. Inst.3, (1941), pp. 21-22
  * [3] N. G. De Bruijn, On Almost Additive Functions, Colloquium Mathematicum , Vol.XV, 1966, 59-63.
  * [4] N. Dunford and J. T. Schwartz, Linear Operators, part 1, Wiley, 1988.
  * [5] I. Karatzas and S. E. Shreve, Brownian Motion and Stochastic Calculus, Springer, 1991.
  * [6] J. Jacod, Calcul Stochastique et Problemes de Martingales, Lecture Notes in Mathematics, Vol 714, Springer, 1979.
  * [7] W. Jurkat, On Cauchy’s functional equation, Proceedings of the American Mathematical Society, Vol. 16, No. 4, 1965, pp. 683-686.
  * [8] K. Kuczma, A Survey of the Theory of Functional Equations, Publikacije Elektrotehničkog fakulteta. Serija Matematika i fizika No. 130 (1964), pp. 1-64.
  * [9] K. Kuczma, Functional equations on restricted domains Survey of the Theory of Functional Equations, Aequat. Math.Vol. 18, (1978), pp. 1-34.
  * [10] R. Sh. Liptser and A. N. Shiryaev, Theory of Martingales, Moskva ”Nauka”, 1986.
  * [11] M. P. Quine, On Three Characterizations of the Normal Distribution, Probability Theory and Mathematical Statistics, Vol 14, Fasc 2, (1993), pp. 257-263.
  * [12] O.M. Schalit, The overdeterminedness of a class of functional equations Aequationes mathematicae, volume 74, (2007), pp. 242–248.
  * [13] S. N. Smirnov, A probabilistic note on the Cauchy functional equation, Aequat. Math.Vol. 93, Issue 2, (2019), pp. 445-449.