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

Clarissa Astuto^*^**Mathematical and Computer Sciences and Engineering Division, King Abdullah University of Science and Technology, Thuwal 23955-6900, Kingdom of Saudi Arabia; [email protected] Jan Haskovec^†^††Mathematical and Computer Sciences and Engineering Division, King Abdullah University of Science and Technology, Thuwal 23955-6900, Kingdom of Saudi Arabia; [email protected] Peter Markowich^‡^‡‡Mathematical and Computer Sciences and Engineering Division, King Abdullah University of Science and Technology, Thuwal 23955-6900, Kingdom of Saudi Arabia; [email protected], and Faculty of Mathematics, University of Vienna, Oskar-Morgenstern-Platz 1, 1090 Vienna; [email protected] Simone Portaro^§^§§Mathematical and Computer Sciences and Engineering Division, King Abdullah University of Science and Technology, Thuwal 23955-6900, Kingdom of Saudi Arabia; [email protected]

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

\displaystyle E[D]=\int_{\Omega}\Phi^{\prime\prime}(u)\nabla u\cdot D\nabla u\mathrm{d}x.

(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

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

(1.2)

where $S=S(x)$ represents the prescribed, time-independent distribution of sources and sinks of the system. We equip (1.2) with the Dirichlet boundary condition

\displaystyle 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.,

\displaystyle\Phi^{\prime}(c)=0.

(1.3)

The formal $L^{2}$ -gradient flow of the energy functional (1.1) constrained by (1.2) is of the form

\displaystyle\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,

(1.4)

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

\displaystyle-\mathrm{div}(D\nabla\sigma)=\Phi^{\prime\prime\prime}(u)\nabla u\cdot D\nabla u\qquad\mathrm{in}\,\,(0,\infty)\times\Omega,

(1.5)

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.4–1.5). 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.4–1.5) 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.4–1.5).

2. The spatially 1D setting

Without loss of generality, let $\Omega=(0,1)$ and the equilibrium $c=0$ . The one–dimensional version of (1.2), (1.5) and (1.4) 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

\displaystyle D(0,x)=D_{I}(x)\geq\gamma>0\quad\forall x\in[0,1],\quad D_{I}\in L^{\infty}(0,1).

(2.2)

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.3)

where $\alpha=\alpha(t)$ is constant in space and

\displaystyle R(x):=\int_{0}^{x}S(y)\,\mathrm{d}y.

Another subsequent integration of (2.3) 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.4)

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

\displaystyle\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,

(2.5)

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.5) and employing (1.3). Finally, substitution of (2.3) and (2.5) 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.6)

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.4). 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.7)

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.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

\displaystyle\|\Phi^{\prime}(u(\cdot,t))\|_{L^{\infty}(0,1)}\leq C(T)\quad\forall t\in[0,T].

(2.8)

Proof.

Fix $\varepsilon>0$ . Using (2.3) 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.9)

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

\displaystyle\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.3) 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.10)

A combination of (2.9) and (2.10) 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.4–2).

Theorem 1.

The system (2.4–2) complemented with the initial condition (2.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.4):


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

complemented with the initial condition (2.2), and the energy functional (1.1) can be written as

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

(2.12)

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.11a), (2.11b) with the initial condition (2.2). If any of the following conditions hold:

(1)

$S(x)\geq 0$ or $S(x)\leq 0$ almost everywhere on $(0,1)$ ,
(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)

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

\displaystyle\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)

(2.13)

\displaystyle\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),

(2.14)

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.11a) 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 (2.13) 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 (2.14) 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.12) under the constraint (2.11), we need to calculate its second-order variation.

Proposition 4.

The second-order variation of (2.12) coupled to (2.11) 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.15)

Moreover, if

(1)

$2\Phi^{\prime\prime}(u)\Phi^{\textrm{iv}}(u)\geq\Phi^{\prime\prime\prime}(u)^{2}$ on $\mathbb{R}$ ,
(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.16)
		$\displaystyle=:u^{0}+\varepsilon u^{1}+\frac{\varepsilon^{2}}{2}u^{2}+O(\varepsilon^{3}).$		(2.17)

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.18)

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.19a)
$\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.19b)
$\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.19c)

Multiplication of (2.1a) by $\Phi^{\prime}(u)$ and integration by parts allow us to recast the energy functional (2.12) 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.18), integration by parts and substitution of (2.19) 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

\displaystyle S(x):=a\delta(x-x_{0})-b\delta(x-x_{1}),\quad a,b>0,\,\,0<x_{0}<x_{1}<1,

(3.1)

where $\delta(x-\overline{x})$ is the Dirac delta distribution centered at $\overline{x}$ . Then

\displaystyle R(x)=a\chi_{[x_{0},x_{1}]}(x)+(a-b)\chi_{[x_{1},1]}(x)\quad\forall x\in[0,1].

(3.2)

Substituting (3.1), (3.2) into (2.11a), 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.6)

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.10)

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.14)

where $D_{1}(t)\geq\gamma>0$ $\forall t\in[0,\infty)$ from (3.6). Consequently, we rewrite (2.11b) 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.18)

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.19a)
$\displaystyle u(x_{1},t)$	$\displaystyle=\frac{(a-b)(1-x_{1})}{D_{2}(t)},$	(3.19b)

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.6) with associated uniformly positive initial condition.

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

(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)

If $\int_{0}^{\infty}\Phi^{\prime\prime}(u)\,\mathrm{d}u\equiv\lim_{u\rightarrow\infty}\Phi^{\prime}(u)=:\alpha<\infty$ and

$\displaystyle\alpha\geq\frac{E[D_{I}]}{(a-b)},$ (3.20)

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$ .

(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)$ .
(4)

If $\int_{-\infty}^{0}\Phi^{\prime\prime}(u)\,\mathrm{d}u\equiv-\lim_{u\rightarrow-\infty}\Phi^{\prime}(u)=:\beta<\infty$ and

$\displaystyle\beta\geq\frac{E[D_{I}]}{(b-a)},$ (3.21)

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.12) 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

\displaystyle 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 (3.20).
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

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

\displaystyle(a-b)\alpha<E[D_{I}]<a\alpha,

(3.22)

for which we currently do not have results. On the other hand, whenever $a<b$ , we write

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

\displaystyle(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.1), 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.6). 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.23)

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.24)

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

We have the following collection of results.

Proposition 6.

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

(1)

If $a>b$ and

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

then $D(x,t)>0$ $\forall(x,t)\in[0,1]\times[0,\infty)$ .
(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.26)

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

If $a>b$ and

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

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

If $a<b$ and

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

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.23) 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.29)

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.25). 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.29) 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.24) of $I(t)$ combined with the bounded limsup hypothesis we must have $I(T)<\infty$ , which leads to a contradiction.
From (3.27), (3.24) we immediately notice that the right-hand side of (3.29) 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.11a–2.11b) and of the ODE presented in (3.6), 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.11a–2.11b), and after we focus on the case (3.1) 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.11a–2.11b), 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

\displaystyle 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
1	1	0	0
1/2	1/4	1/4	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.1) and

\Delta t=10^{-6},h=10^{-4}

In Fig. 2 we show the solution ${\bf D}$ (left part of each plot, in blue), together with ${\bf\Phi}^{\prime\prime}(u)$ (right part of each plot, in orange) with two different initial conditions: ${\bf D}(t=0)=1$ in (a) and ${\bf D}(t=0)=0.1$ in (b) and (c) panels. The expression for the second derivative of the entropy function is $\Phi^{\prime\prime}(u)=\exp(-u)$ , while $m=-1.98$ and $q=1$ in the expression (4.3). In these plots we show ${\bf D}$ and ${\bf\Phi}^{\prime\prime}(u)$ in different moments, to show how they evolve in time. It is interesting to see that after finite time the solution ${\bf D}$ touches zero if the initial condition is small enough. We show ${\bf\Phi}^{\prime\prime}(u)$ in the same plot because, analyzing equation (2.11a), it is evident that ${\bf\Phi}^{\prime\prime}(u)$ acts as a weight, dependent on the solution itself ${\bf D}$ , for the integral involving $S(x)$ . To induce a change in sign in ${\bf D}$ , it is necessary for the integral to become negative at a rate that is sufficient to drive the solution towards zero. Consequently, we select $S(x)$ to ensure that ${\bf\Phi}^{\prime\prime}(u)$ becomes significantly big when $S(x)$ is negative.

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.4)), 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.1) 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.1) 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).

Refer to caption — 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.6), where the choice for the sink/source function is defined as a sum of delta functions, in (3.1). Here we rewrite (3.6) 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

\displaystyle\Phi^{\prime\prime}(u)=(u^{2}(\sin(u)+\sigma)+1)^{-1},

(4.19)

where $\sigma$ is selected such that

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

ensuring that we do not fall into the case (3.27) 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.20)
	$\displaystyle\frac{\partial D_{2}}{\partial t}=$	$\displaystyle M_{2}(D_{1},D_{2})\,D_{2}\qquad x\in[x_{1},1].$		(4.21)

where

	$\displaystyle M_{1}(D_{1},D_{2})$	$\displaystyle=\frac{1}{D_{1}^{3}}a^{2}\,\Phi^{\prime\prime}(u(x_{0},t))$		(4.22)
	$\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.23)

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.24)
$\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.25)
$\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.26)
$\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.27)

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.28)
	$\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.29)

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.28), 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=260.17464pt,unit=1mm,scale={.25}]{Figures/error_plot_D1D2_433_quadrat} \end{overpic}

Figure 7. Accuracy plot of the numerical scheme defined in (4.28), 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=433.62pt,unit=1mm,scale={.25}]{Figures/D1_D2_time} \put(1.0,54.0){(a)} \end{overpic}

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

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

\begin{overpic}[abs,width=433.62pt,unit=1mm,scale={.25}]{Figures/minD} \put(1.0,54.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).
[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).
[4] C. Astuto, D. Boffi, and F. Credali, Finite element discretization of a biological network formation system: a preliminary study, (submitted).
[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.
[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.

	$\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,$

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

abstract

1. Introduction

2. The spatially 1D setting

Proposition 1.

Proof.

Proposition 2.

Proof.

Theorem 1.

Proposition 3.

Proof.

Proposition 4.

Proof.

3. A remarkable example

Proposition 5.

Proof.

Remark 1.

Proposition 6.

Proof.

4. Numerical results

4.1. Discretization in space and time

4.1.1. Accuracy test and Results

4.2. Time discretization for the δ\delta-system

4.2.1. Accuracy test and Results

5. Conclusions

References

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