Semilinear multi-term fractional in time diffusion with memory

In this study, the initial-boundary value problems to semilinear integro-differential equations with multi-term fractional Caputo derivatives are analyzed. A particular case of these equations models oxygen diffusion through capillaries. Under proper requirements on the given data in the model, the classical and strong solvability of these problems for any finite time interval [0, T] are proved via so-called continuation method. The key point in this approach is finding suitable a priori estimates of a solution in the fractional Hölder and Sobolev spaces.

1 Introduction

Complex phenomena in the engineering and scientific fields are modeled utilizing the fractional differential equations (FDEs). Nowadays, the fractional calculus is an efficient tool for describing dynamic behavior of living systems and hereditary properties of various materials: the relaxation process in polymers [1], chaotic neuron model [2], longtime memory in financial time series via fractional Langevin equations [3], and tumor growth models [4] (see also references therein). We also refer to [5, 6], where the authors propose the advanced mathematical model for oxygen delivery to tissue through a capillary in both (transverse and longitudinal) directions. In these studies, conveying oxygen from a capillary to the surrounding tissue is described by means of a subdiffusion equation having two fractional derivatives in time, that is

$\begin{array}{l}{D}_t^{\nu }\text{C}-\tau D_t^{\mu }\text{C}=\text{div}(a_0\nabla \text{C})-\text{k}-{\displaystyle {\int }_0^t\frac{{(t-s)}^{\nu -1}}{\Gamma (\nu )}}(b_1(x,s)\nabla \text{C}(x,s)\\ +b_0(x,s)\text{C}(x,s))ds\end{array}$

with 0 < μ < ν < 1. Here, C represents the concentration of oxygen, τ is the time lag in concentration of oxygen along the capillaries (in the present model, this parameter is a positive constant), k is the rate of consumption per volume of tissue, and a₀ and b_i are the diffusion coefficients of oxygen. In addition, the term $D_t^{\nu }\text{C}-\tau D_t^{\mu }\text{C}$ details the net diffusion of oxygen to all tissues.

In this equation, the symbol $D_t^{\theta }$ stands for the Caputo fractional derivative with respect to time of order θ ∈ (0, 1),

$D_t^{\theta }\text{C}(x,t)=\{\begin{array}{l}\frac{1}{\Gamma (1-\theta )}\frac{\partial }{\partial t}{\displaystyle \underset{0}{\overset{t}{\int }}\frac{\text{C}(x,s)-\text{C}(x,0)}{{(t-s)}^{\theta }}}ds\text{\hspace{1em}if\hspace{1em}}\theta \in (0,1),\hfill \\ \frac{\partial \text{C}}{\partial t}(x,t)\text{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} if\hspace{1em}}\theta =1,\hfill \end{array}$

where Γ is the Euler's Gamma function. An equivalent definition of this derivative in the case of absolutely continuous functions reads

$D_t^{\theta }\text{C}(x,t)=\{\begin{array}{l}\frac{1}{\Gamma (1-\theta )}{\displaystyle \underset{0}{\overset{t}{\int }}{(t-s)}^{-\theta }}\frac{\partial \text{C}}{\partial s}(x,s)ds\text{\hspace{1em}if\hspace{1em}}\theta \in (0,1),\hfill \\ \frac{\partial \text{C}}{\partial t}(x,t)\text{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} if\hspace{1em}}\theta =1.\hfill \end{array}$

In this art, motivated by the discussion above, we focus on the analytical study of the semilinear integro-differential equation with memory terms. Let Ω ⊂ ℝⁿ, n ≥ 2, be a bounded domain with a smooth boundary ∂Ω, and for any T > 0, we set

$\begin{array}{l}{\Omega }_T=\Omega \times (0,T)\text{and\hspace{1em} }\partial {\Omega }_T=\partial \Omega \times \left[0,T\right].\end{array}$

We consider the initial-value problems to the multi-term time-fractional semilinear diffusion equation in the unknown function u = u(x, t):Ω_T → ℝ,

$\begin{array}{ll}{D}_{t}u-{ℒ}_{1}u-\mathcal{K}*{ℒ}_{2}u+f(u)=g(x,t)\text{\hspace{1em}in\hspace{1em}}{\Omega }_T,& (1.1)\end{array}$

subject to the following initial and boundary conditions:

$\begin{array}{ll}\{\begin{array}{l}u(x,0)=u_0(x)\text{\hspace{1em}in\hspace{1em}}\overline{\Omega },\hfill \\ u={\phi }_1(x,t)\text{\hspace{1em}on\hspace{1em}}\partial {\Omega }_T\text{\hspace{1em}in }DBC\text{ case},\hfill \\ \text{or}\hfill \\ ℳu+{\mathcal{K}}_1*ℳu-c_{0}u={\phi }_2(x,t)\text{\hspace{1em}on\hspace{1em}}\partial {\Omega }_T\text{\hspace{1em}in }3BC\text{ case},\hfill \end{array}& (1.2)\end{array}$

where the abbreviations DBC and 3BC mean the Dirichlet boundary condition and the boundary condition of the third kind, respectively.

Here, c₀ is given positive number, g, u₀, φ_i are given functions, and ${\mathcal{K}}_1$ and $\mathcal{K}$ are prescribed memory kernels.

Here, the symbol * stands for the usual time-convolution product on (0, t),

$({\mathfrak{h}}_1*{\mathfrak{h}}_2)(t)={\displaystyle \underset{0}{\overset{t}{\int }}{\mathfrak{h}}_1}(t-s){\mathfrak{h}}_2(s)ds.$

The operator D_t is the linear combination of Caputo fractional derivatives with respect to time, namely

$\begin{array}{ll}{D}_{t}u=D_t^{\nu }(\rho u)+{\displaystyle \sum _{i=1}^M{D}_t^{{\nu }_i}}({\rho }_{i}u)-{\displaystyle \sum _{j=1}^N{D}_t^{{\mu }_j}}({\gamma }_{j}u),& (1.3)\end{array}$

where ν ∈ (0, 1) and ν_i, μ_j ∈ (0, ν) are arbitrary but fixed, and ρ = ρ(x, t), ρ_i = ρ_i(x, t) and γ_j = γ_j(x, t) are given positive functions.

Coming to the remaining operators, ${\mathcal{L}}_i$, i = 1, 2, are linear elliptic operators of the second order with time-dependent coefficients, while $\mathcal{M}$ is a first-order differential operator. Their precise forms will be given in Sections 3, where we detail the main assumptions in the model.

Published works concerning the multi-term fractional diffusion/wave equations, i.e., the equation with the operator

$\begin{array}{ll}{D}_{t}u={\displaystyle \sum _{i=1}^N{q}_i}{D}_t^{{\nu }_i}u,& (1.4)\end{array}$

with q_i being positive, and 0 ≤ ν₁ < ν₂ < ... < ν_M, are quite limited in spite of rich literature on their single-term version. Exact solutions of linear multi-term fractional diffusion equations with q_i being positive constants on bounded domains are searched employing eigenfunction expansions in Daftardar-Gejji and Bhalekar [7] and Morales-Delgado et al. [5]. We quote Srivastava and Rai [6] and Morales-Delgado et al. [5], where certain numerical solutions are constructed to the corresponding initial-boundary value problems to evolution equations with D_t given via Equation 1.4. Finally, we mention [8], where existence and non-existence of the mild solutions to the Cauchy problem for semilinear subdiffusion equation with the operator Equation 1.4 are discussed. In particular, the authors obtain the Fujita-type and Escobedo-Herrero-type critical exponents for this equation and the system. It is worth noting that, all these works concern to evolution equations with the operator Equation 1.4 which can be rewritten in the form of a generalized fractional derivative with a non-negative locally integrable kernel $\mathcal{N}(t)$, that is

$\begin{array}{ll}{\text{D}}_{t}u(x,t)=\frac{\partial }{\partial t}{\displaystyle {\int }_0^t}\mathcal{N}(t-\tau )u(x,\tau )d\tau -\mathcal{N}(t)u(x,0),t>0.& (1.5)\end{array}$

Coming to the initial-boundary value problems associated with Equation 1.1 with the operator D_t given by Equation 1.3, we point out two principal differences with respect to the aforementioned articles. The first deals with the presence of Caputo fractional derivatives of the product of two functions: the desired solution u and the prescribed coefficients ρ, ρ_i, γ_j. Incidentally, we recall that the well-known Leibniz rule does not work in the case of fractional derivatives. The second distinction is that the operator D_t given by Equation 1.3 (under certain assumptions on the coefficients) can be represented in the form Equation 1.4 but with a negative kernel. Indeed, setting in Equation 1.3

$\begin{array}{l}M=0,N=1,\rho =C_{\rho },{\gamma }_1=1+C_{\rho },\end{array}$

where C_ρ is a positive constant, we have the representation

$\begin{array}{l}{\text{D}}_{t}u=\frac{\partial }{\partial t}(\mathcal{N}*(u-u_0))\end{array}$

with

$\mathcal{N}=C_{\rho }[\frac{t^{-\nu }}{\Gamma (1-\nu )}-\frac{t^{-{\mu }_1}}{\Gamma (1-{\mu }_1)}]-\frac{t^{-{\mu }_1}}{\Gamma (1-{\mu }_1)}$

being negative for $t>e^{-C_{\gamma }}$ [see Janno and Kinash [9], Lemma 4]; C_γ is the Euler–Mascheroni constant. In fact, the non-negativity of the kernel $\mathcal{N}$ is a principal assumption in the aforementioned studies.

The linear version of Equations 1.1, 1.3 subject to various type boundary conditions with the coefficients in D_t being alternating sign is discussed in Pata et al. [10] and Vasylyeva [11]. For any fixed time T, the existence and uniqueness of a solution to semilinear problem (Equations 1.1, 1.3) with the Dirichlet or the Neumann boundary conditions are analyzed in Siryk and Vasylyeva [12] and Vasylyeva [11]. Namely, if the coefficients of the operator D_t are only time-dependent and non-decreasing functions, then the well-posedness of these problems in the fractional Hölder and Sobolev spaces is established in the one-dimensional case in Siryk and Vasylyeva [12]. As for the multidimensional case, the classical solvability of the Cauchy-Dirichlet problem to semilinear (Equation 1.1) in the case of two-term fractional derivatives in D_t (i.e., either M = 1, N = 0 or M = 0, N = 1) is proved in Vasylyeva [11]. In this study, the coefficients in Equation 1.1 are time and space dependent but instead of their non-decreasing in time, they have to satisfy more complex assumption. Indeed, if M = 1, N = 0, then the function $\frac{\rho }{{\rho }_1}$ should be decreasing. Finally, we remark that in Vasylyeva [11], the non-linear term is the local Lipschitz.

The goal of this study is finding sufficient conditions on the coefficients of the operator D_t, the order fractional derivatives ν, μ_j, and ν_i, j, i ≥ 1, and the non-linearity f which provide one-to-one classical and strong solvability (for any fixed T) in the case of the DBC or the 3BC. Actually, we consider two types of the non-linearity f(u). The first is f satisfying the local Lipschitz condition and having the linear growth. As for the second, f is a continuous differentiable on ℝ with a super-linear growth. For example, f is a polynomial of odd degree with the positive leading coefficient (see Giorgi et al. [13]). Coming to the coefficients in the fractional operator D_t, we discuss both the non-decreasing coefficients and the coefficients satisfying the properties of Theorem 2 [11].

We notice that the key ingredient in the proof of the classical solvability is the continuation approach, based on the introduction of a family of auxiliary problems depending on a parameter λ ∈ [0, 1]. Then, one has to produce a priori estimates in the fractional Hölder spaces for the solution which are independent of λ. One of the crucial points in the arguments is concerned to the estimates of $\left|\right|u|{|}_{\mathcal{C}({\overline{\Omega }}_T)}$, obtained via integral iteration technique adopted to the multi-term fractional case. As for the strong solvability, it is proved via the construction of this solution as a limit of approximate smooths solutions and exploiting a priori estimates in the Sobolev spaces.

Finally, we notice that assumptions on the coefficients and the memory kernels in the one-dimensional and multidimensional cases are different. It is related with using various approaches to get a priori estimates of the solutions if n = 1 and n ≥ 2. Namely, if n ≥ 2, we relax assumptions on the coefficients of D_t, in particular, we allow coefficients depending on time and space in Equation 1.3. However, we require more regular memory kernel in Equation 1.1, $\mathcal{K}\in {\mathcal{C}}^1(\left[0,T\right])$.

Outline of the study

This article is organized as follows: in Section 2, we introduce the notations and the functional spaces. The general assumptions and main results (Theorems 3.1, 3.2) are stated in Section 3. Theorem 3.1 is devoted to the one-valued classical solvability to Equation 1.1 with the DBC or the 3BC in the multidimensional case, while the strong solvability is established in Theorem 3.2. Section 4 is auxiliary and contains some technical and preliminary results from fractional calculus, playing a key role in the course of the investigation. Section 5 concerns to the obtaining a priori estimates in the fractional Hölder and Sobolev spaces, which will be a crucial point in the proof of the main results. Here, the key bound is the estimate of $\Vert u{\Vert }_{{\mathcal{C}}^{\alpha ,\alpha \nu /2}({\overline{\Omega }}_T)}$, produced via integral iteration techniques adapted to the case of multi-term fractional derivatives. The proof of Theorems 3.1 and 3.2 is carried out in Section 6.

2 Functional spaces and notation

Throughout this study, the symbol C will denote a generic positive constant, depending only on the structural quantities of the problem.

In the course of our study, we will exploit the fractional Hölder and Sobolev spaces. To this end, in what follows, we take two arbitrary (but fixed) parameters

$\begin{array}{l}\alpha \in (0,1)\text{and\hspace{1em}}\nu \in (0,1).\end{array}$

For any non-negative integer l, any p ≥ 1, s ≥ 0, and any Banach space (X, ||·||_X), we consider the usual spaces

${\mathcal{C}}^{l+\alpha }(\overline{\Omega }),W^{s,p}(\Omega ),L_p(\Omega ),{\mathcal{C}}^s([0,T],X),W^{s,p}((0,T),X).$

Recall that for non-integer s, the space W^{s, p} is called Sobolev-Slobodeckii space [for its definition and properties see, e.g., Adams and Fournier [14], Chapter 1].

Denoting for β ∈ (0, 1)

$\begin{array}{c}{\langle v\rangle }_{x,{\Omega }_T}^{(\beta )}\\ =\sup \{\frac{|v(x_1,t)-v(x_2,t)|}{|x_1-x_2{|}^{\beta }}:x_2\ne x_1,x_1,x_2\in \overline{\Omega },\text{\hspace{1em}}t\in [0,T]\},\\ {\langle v\rangle }_{t,{\Omega }_T}^{(\beta )}\\ =\sup \{\frac{|v(x,t_1)-v(x,t_2)|}{|t_1-t_2{|}^{\beta }}:t_2\ne t_1,x\in \overline{\Omega },t_1,t_2\in [0,T]\}.\end{array}$

Then, we assert the following definition.

Definition 2.1. A function v = v(x, t) belongs to the class ${\mathcal{C}}^{l+\alpha ,\frac{l+\alpha }{2}\nu }({\bar{\Omega }}_T)$, for l = 0, 1, 2, if the function v and its corresponding derivatives are continuous and the norms here below are finite:

$\begin{array}{c}\Vert v{\Vert }_{{\mathcal{C}}^{l+\alpha ,\frac{l+\alpha }{2}\nu }({\overline{\Omega }}_T)}\\ =\{\begin{array}{l}\Vert v{\Vert }_{\mathcal{C}([0,T],{\mathcal{C}}^{l+\alpha }(\overline{\Omega }))}+{\displaystyle {\sum }_{\left|j\right|=0}^l{\langle D_x^{j}v\rangle }_{t,{\Omega }_T}^{(\frac{l+\alpha -\left|j\right|}{2}\nu )}},l=0,1,\hfill \\ \Vert v{\Vert }_{\mathcal{C}([0,T],{\mathcal{C}}^{2+\alpha }(\overline{\Omega }))}+\Vert D_t^{\nu }v{\Vert }_{{\mathcal{C}}^{\alpha ,\frac{\alpha }{2}\nu }({\overline{\Omega }}_T)}\hfill \\ +{\displaystyle {\sum }_{\left|j\right|=1}^2{\langle D_x^{j}v\rangle }_{t,{\Omega }_T}^{(\frac{2+\alpha -\left|j\right|}{2}\nu )}},l=2.\hfill \end{array}\end{array}$

In a similar way, for l = 0, 1, 2, we introduce the space ${\mathcal{C}}^{l+\alpha ,\frac{l+\alpha }{2}\nu }(\partial {\Omega }_T)$.

The properties of these spaces have been discussed in Krasnoschok et al. [15] (Section 2). As for the limiting case ν = 1, these classes boil down to the usual parabolic Hölder spaces.

Finally, we will say that a function v defined in Ω_T belongs to ${ℌ}_p^{s_1,s_2}({\Omega }_T)$ with p > 1 and s₁, s₂ ≥ 0, if $v\in W^{s_1,p}((0,T),L_p(\Omega ))\cap L_p((0,T),W^{s_2,p}(\Omega ))$, and the norm here below is finite

$\begin{array}{l}\left|\right|v|{|}_{{ℌ}_p^{s_1,s_2}({\Omega }_T)}=|\left|v\right|{|}_{W^{s_1,p}((0,T),L_p(\Omega ))}+\left|\right|v|{|}_{L_p((0,T),W^{s_2,p}(\Omega ))}.\end{array}$

The space ${ℌ}_p^{s_1,s_2}(\partial {\Omega }_T)$ is defined in the similar manner.

3 Main results

First, we state additional requirements on the given data in Equations 1.1, 1.2.

• h1 (Conditions on the fractional order of the derivatives in Equation 1.3): We assume that

$\begin{array}{c}\nu \in (0,1),{\nu }_i,{\mu }_j\in (0,\frac{\nu (2-\alpha )}{2}),{\nu }_i\ne {\mu }_j,\\ i=1,2,\mathrm{...},M,j=1,2,\mathrm{...},N,\\ 0<{\nu }_1<{\nu }_2<\mathrm{...}<{\nu }_M<\nu ,0<{\mu }_1<{\mu }_2<\mathrm{...}<{\mu }_N<\nu .\end{array}$

• h2 (Conditions on the operators): The operators appearing in Equations 1.1, 1.2 read

$\begin{array}{ll}\{\begin{array}{l}{ℒ}_1=\sum _{ij=1}^n\frac{\partial }{\partial x_i}(a_{ij}(x,t)\frac{\partial }{\partial x_j})+\sum _{i=1}^n{a}_i(x,t)\frac{\partial }{\partial x_i}+a_0(x,t),\hfill \\ {ℒ}_2=\sum _{ij=1}^n\frac{\partial }{\partial x_i}(a_{ij}(x,t)\frac{\partial }{\partial x_j})+\sum _{i=1}^n{b}_i(x,t)\frac{\partial }{\partial x_i}+b_0(x,t),\hfill \\ ℳ=-\sum _{ij=1}^n{a}_{ij}(x,t)N_i\frac{\partial }{\partial x_j},\hfill \end{array}& (3.1)\end{array}$

where N_i is a component of the outward normal N = {N₁, ..., N_n} to Ω; the fractional operator D_t in Equation 1.1 is given by Equation 1.3.

There are positive constants 0 < δ₁ < δ₂, such that

$\begin{array}{l}{\delta }_1|\xi {|}^2\le {\displaystyle \sum _{ij=1}^n}{a}_{ij}(x,t){\xi }_i{\xi }_j\le {\delta }_2|\xi {|}^2\end{array}$

for any $(x,t,\xi )\in {\overline{\Omega }}_T\times {ℝ}^n$.

Moreover, we require that

$\begin{array}{c}{a}_0,b_0\in {\mathcal{C}}^{\alpha ,\alpha \nu /2}({\overline{\Omega }}_T),a_{ij},a_i,b_j\in {\mathcal{C}}^{1+\alpha ,(1+\alpha )\nu /2}({\overline{\Omega }}_T),\\ i,j=1,\mathrm{...},n.\end{array}$

• h3 (Conditions on the coefficients of D_t): We require that for

$\begin{array}{l}{\nu }_0\ge max\left\{1,\nu (2+\alpha )/2\right\}\end{array}$

the relations hold

$\begin{array}{l}\rho (x,t),{\rho }_i(x,t),{\gamma }_j(x,t)\in {\mathcal{C}}^{{\nu }_0}(\left[0,T\right],{\mathcal{C}}^1(\overline{\Omega }));\end{array}$

and there are positive constants δ, δ₃, δ₄, such that

$\begin{array}{l}\rho \ge \delta >0,{\rho }_i\ge {\delta }_3>0,{\gamma }_j\ge {\delta }_4>0\end{array}$

for each $(x,t)\in {\overline{\Omega }}_T$ and for all i = 1, 2, ..., M, j = 1, 2, ..., N.

In addition, if N ≥ 1, then

$\begin{array}{l}\rho (x,t)={\rho }_0(x,t)+{\displaystyle \sum _{j=1}^N}{\gamma }_j(x,t),\end{array}$

where the function ${\rho }_0\in {\mathcal{C}}^{{\nu }_0}(\left[0,T\right],{\mathcal{C}}^1(\overline{\Omega }))$ is positive for all t ∈ [0, T] and $x\in \overline{\Omega }$.

Moreover, we require that the one of the following conditions holds:

(i) either $\frac{\partial \rho }{\partial t},$ $\frac{\partial {\rho }_0}{\partial t},$ $\frac{\partial {\rho }_i}{\partial t},$ $\frac{\partial {\gamma }_j}{\partial t}$ are non-negative for all $(x,t)\in {\overline{\Omega }}_T$;

(ii) or

$\{\begin{array}{l}\frac{\partial }{\partial t}\left(\frac{{\rho }_0}{{\rho }_i}\right),\frac{\partial }{\partial t}\left(\frac{{\rho }_0}{{\gamma }_j}\right)\le 0\text{\hspace{1em}if\hspace{1em}\hspace{1em}}N\ge 1,\hfill \\ \frac{\partial }{\partial t}\left(\frac{\rho }{{\rho }_i}\right)\le 0,\text{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} if\hspace{1em}}N=0,\hfill \end{array}$

for all i = 1, ..., M, j = 1, ..., N, and any $(x,t)\in {\overline{\Omega }}_T$.

• h4 (Conditions on the right-hand sides): The given functions have the following regularity:

$\begin{array}{r}g\in {\mathcal{C}}^{\alpha ,\frac{\nu \alpha }{2}}({\overline{\Omega }}_T),u_0\in {\mathcal{C}}^{2+\alpha }(\overline{\Omega }),\\ {\phi }_1\in {\mathcal{C}}^{2+\alpha ,\frac{2+\alpha }{2}\nu }(\partial {\Omega }_T),{\phi }_2\in {\mathcal{C}}^{1+\alpha ,\frac{1+\alpha }{2}\nu }(\partial {\Omega }_T),\end{array}$

• h5 (Conditions on the memory kernels):

$\begin{array}{l}\mathcal{K}(t)\in {\mathcal{C}}^1(\left[0,T\right]),{\mathcal{K}}_1\in L_1(0,T).\end{array}$

• h6 (Compatibility conditions): The following compatibility conditions hold for every x ∈ ∂Ω at the initial time t = 0,

$\begin{array}{l}\begin{array}{c}{\phi }_1(x,0)=u_0(x)\text{and}\\ {\text{D}}_t{\phi }_1{|}_{t=0}={\mathcal{L}}_1{u}_0(x){|}_{t=0}-f(u_0)+g(x,0),\end{array}\end{array}$

if the DBC holds, and there is

$\begin{array}{l}\mathcal{M}{u}_0(x){|}_{t=0}-{\mathfrak{c}}_0{u}_0(x)={\phi }_2(x,0)\end{array}$

in the 3BC case.

• h7 (Conditions on the nonlinearity): We assume that the one of the following requirements holds:

• h7.I: either f(u) is the local Lipschitz and has a linear growth, i.e., for every ϱ > 0, there exists a positive constant C_ϱ, such that

$\begin{array}{l}|f(u_1)-f(u_2)|\le C_{\varrho }|u_1-u_2|\end{array}$

for any u₁, u₂ ∈ [−ϱ, ϱ]; and

there is a positive constant L, such that

$\begin{array}{l}\left|f(u)\right|\le L(1+\left|u\right|)\text{for any\hspace{1em}}u\in ℝ;\end{array}$

• h7.II: or $f\in {\mathcal{C}}^1(ℝ)$, and for some non-negative constants L_i, i = 1,2,3,4, and q ≥ 0, the inequalities hold

$\{\begin{array}{l}|f(u)|\le L_1(1+|u{|}^q),\hfill \\ uf(u)\ge -L_2+L_3|u{|}^{q+1},\hfill \\ f^{\prime }(u)\ge -L_4.\hfill \end{array}$

Remark 3.1. It is apparent that if the positive functions ρ, ρ_i, γ_j are time-independent, then condition h3(i) boils down to h3(ii).

Example 3.1. The simplest example of the functions satisfying h3 is

$\begin{array}{l}\rho ={\mathcal{C}}_0,{\gamma }_j={\mathcal{C}}_j,{\rho }_i={\overline{\mathcal{C}}}_i,i=1,...,M,j=1,...,N,\end{array}$

where ${\mathcal{C}}_0$, ${\mathcal{C}}_j,$ ${\overline{\mathcal{C}}}_i$ are positive constants, such that

$\begin{array}{l}{\mathcal{C}}_0-{\displaystyle \sum _{j=1}^N}{\mathcal{C}}_j>0.\end{array}$

Now, we are in the position to state the one-valued classical solvability of Equations 1.1, 1.2.

Theorem 3.1. Let T > 0 be arbitrarily given, $\partial \Omega \in {\mathcal{C}}^{2+\alpha },$ n ≥ 2, and let assumptions h1–h6 hold. We assume that f(u) meets the requirement h7.I if N ≥ 1, while in the case of N = 0, f(u) satisfies h7. Then, initial-boundary value problem Equations 1.1, 1.2 admits a unique classical solution u = u(x, t) satisfying the regularity:

$\begin{array}{r}u\in {\mathcal{C}}^{2+\alpha ,\frac{2+\alpha }{2}\nu }({\overline{\Omega }}_T),{\text{D}}_t^{{\nu }_i}u,{\text{D}}_t^{{\mu }_j}u\in {\mathcal{C}}^{\alpha ,\frac{\alpha \nu }{2}}({\overline{\Omega }}_T),\\ i=1,...,M,j=1,...N.\end{array}$

The next assertion is related to the strong solvability of Equations 1.1, 1.2.

Theorem 3.2. Let N = 0, n ≥ 2, $\partial \Omega \in {\mathcal{C}}^{2+\alpha },$ and let T > 0 be arbitrarily given. We assume that h1–h5 and h7 hold and

$\begin{array}{l}{\psi }_1,{\psi }_2,u_0\equiv 0,f\in L_p({\Omega }_T)\cap W^{s_1,r}((0,T),W^{s_2,r}(\Omega )),\end{array}$

where $p>max\left\{n+\frac{2}{\nu };\frac{1}{\nu -{\nu }_M}\right\}$, r ≥ n + 1, $s_1\in (r^{-1},1)$, and $s_2\in ((n+1)r^{-1},1)$. Moreover, in the DBC case, we require

$\begin{array}{l}f(0){|}_{\partial \Omega }=g(x,0){|}_{\partial \Omega }.\end{array}$

Then, the initial-boundary value problem Equations 1.1, 1.2 admits a unique strong solution in the class ${ℌ}_p^{\nu ,2}({\Omega }_T)$.

Remark 3.2. Theorems 3.1 and 3.2 hold if D_tu in Equation 1.3 is changed by

$\begin{array}{l}{\text{D}}_{t}u=\rho (x,t){\text{D}}_t^{\nu }u+{\displaystyle \sum _{i=1}^M}{\rho }_i(x,t){\text{D}}_t^{{\nu }_i}u-{\displaystyle \sum _{j=1}^N}{\gamma }_j(x,t){\text{D}}_t^{{\mu }_j}u,\end{array}$

where ρ, ρ_i, γ_j satisfies h3, but the requirement on the regularity of these functions can be relaxed. Namely, $\rho ,{\rho }_i,{\gamma }_j\in {\mathcal{C}}^{\alpha ,\alpha \nu /2}({\overline{\Omega }}_T)$.

The remaining part of this study is devoted to the verification of Theorems 3.1, 3.2. Here, we proceed with a detailed proof of Theorem 3.1 in the most difficult case, i.e., if N ≥ 1, M ≥ 1 in Equation 1.3. This means that the non-linear term f(u) satisfies h7.I. The verification of the remaining cases is simpler and repeats the main steps (with minor changes) in the arguments related with the cases N, M ≥ 1.

4 Technical results

In this section, we collect some useful properties of fractional derivatives and integrals, as well as several preliminaries results that will be significant in our investigation. Throughout this art, for any θ > 0, we use the notation

$\begin{array}{l}{\omega }_{\theta }=\frac{t^{\theta -1}}{\Gamma (\theta )}\end{array}$

and define the fractional Riemann-Liouville integral and the derivative of order θ, respectively, of a function v = v(x, t) with respect to time t as

$\begin{array}{l}{I}_t^{\theta }v(x,t)=({\omega }_{\theta }*v)(x,t){\partial }_t^{\theta }v(x,t)=\frac{{\partial }^{\lceil \theta \rceil }}{\partial t^{\lceil \theta \rceil }}({\omega }_{\lceil \theta \rceil -\theta }*v)(x,t),\end{array}$

where ⌈θ⌉ is the ceiling function of θ (i.e., the smallest integer greater than or equal to θ).

It is apparent that, for θ ∈ (0, 1), there holds

$\begin{array}{l}{\partial }_t^{\theta }v(x,t)=\frac{\partial }{\partial t}({\omega }_{1-\theta }*v)(x,t).\end{array}$

Accordingly, the Caputo fractional derivative of the order θ ∈ (0, 1) to the function v(x, t) can be represented as

$\begin{array}{ll}\begin{array}{l}{D}_t^{\theta }v(·,t)=\frac{\partial }{\partial t}({\omega }_{1-\theta }*v)(·,t)-{\omega }_{1-\theta }(t)v(·,0)\hfill \\ ={\partial }_t^{\theta }v(·,t)-{\omega }_{1-\theta }(t)v(·,0)\hfill \end{array}& (4.1)\end{array}$

provided that both derivatives exist.

At this point, we subsume [16, Proposition 4.1], [11, Proposition 1] as the following claim.

Proposition 4.1. The following hold.

(i) For any given positive numbers θ₁ and θ₂ and a summable kernel k = k(t), there are relations

$\begin{array}{l}({\omega }_{{\theta }_1}*{\omega }_{{\theta }_2})(t)={\omega }_{{\theta }_1+{\theta }_2}(t),(1*{\omega }_{{\theta }_1})(t)={\omega }_{1+{\theta }_1}(t),\\ {\omega }_{{\theta }_1}(t)\ge C{T}^{{\theta }_1-1},({\omega }_{{\theta }_1}*k)(t)\le C{\omega }_{{\theta }_1}(t).\end{array}$

Here, the positive constant C depends only on T, θ₁, and ||k||_{L1(0, T)}.

(ii) Let $k(t)\in {\mathcal{C}}^1(\left[0,T\right])$, θ ∈ (0, 1), θ₁ ≥ 1, $v=v(t)\in {\mathcal{C}}^{\theta }(\left[0,T\right]),$ ${\text{D}}_t^{\theta }v(t)\in \mathcal{C}(\left[0,T\right])$, $w=w(t)\in {\mathcal{C}}^{{\theta }_1}(\left[0,T\right])$. Then, the equality holds

$\begin{array}{l}(k*w{D}_t^{\theta }v)(t)=k(0)w(t)({\omega }_{1-\theta }*[v-v(0)])(t)\\ +(k^{\prime }*w({\omega }_{1-\theta }*[v-v(0)]))(t)\\ +(k*w^{\prime }({\omega }_{1-\theta }*[v-v(0)]))(t),t\in [0,T].\end{array}$

The next result is key inequalities in the fractional calculus and includes [12, Proposition 5.1, Corollaries 5.2-5.3].

Proposition 4.2. The following holds.

(i) Let θ, θ₁ ∈ (0, 1) and θ₁ > θ/2, $v\in {\mathcal{C}}^{{\theta }_1}(\left[0,T\right])$. For any even integer p ≥ 2, the inequalities are true

${\partial }_t^{\theta }{v}^p(t)\le {\partial }_t^{\theta }{v}^p(t)+(p-1)v^p(t){\omega }_{1-\theta }(t)\le p{v}^{p-1}(t){\partial }_t^{\theta }v(t).$

If v is non-negative, then these bounds hold for any integer odd p.

(ii) Let 0 < θ₁ < θ ≤ 1, θ₂ ∈ (θ₁, 1), and $v\in {\mathcal{C}}^{{\theta }_2}(\left[0,T\right]).$ Then, there is positive value T₁ = T₁(θ), such that the following inequalities hold:

$\begin{array}{l}\begin{array}{c}{\mathcal{N}}_1=\mathcal{N}(t;\theta ,{\theta }_1)={\omega }_{1-\theta }(t)-{\omega }_{1-{\theta }_1}(t)\ge 0\text{for all}\\ t\in \left[0,T_1\right];\\ \frac{d}{dt}({\mathcal{N}}_1*v^p)(t)\le \frac{d}{dt}({\mathcal{N}}_1*v^p)(t)+(p-1)v^p(t){\mathcal{N}}_1(t)\\ \le p{v}^{p-1}(t)\frac{d}{dt}({\mathcal{N}}_1*v)(t)\text{for all}\\ t\in \left[0,min\left\{T,T_1\right\}\right],\end{array}\end{array}$

where p meets requirements of (i).

At this point, for given functions w₁ and w₂, we define

$\begin{array}{l}{\mathcal{J}}_{\theta }(t)={\mathcal{J}}_{\theta }(t;w_1,w_2)={\displaystyle {\int }_0^t\frac{[w_1(t)-w_1(s)]}{{(t-s)}^{1+\theta }}}[w_2(s)-w_2(0)]ds,\\ \mathcal{W}(w_1)=\mathcal{W}(w_1;t,\tau )={\displaystyle {\int }_0^1\frac{\partial w_1}{\partial z}}(z)ds,\text{\hspace{1em}where}\\ z=st+(1-s)\tau ,0<\tau <t<T,\end{array}$

and assert the results obtained in ([12], Proposition 5.5) and related to the fractional differentiation of the product.

Proposition 4.3. Let θ ∈ (0, 1), $w_1\in {\mathcal{C}}^1(\left[0,T\right]),$ $w_2\in \mathcal{C}(\left[0,T\right])$.

(i) If ${\text{D}}_t^{\theta }{w}_2$ belongs either to $\mathcal{C}(\left[0,T\right])$ or to L_p(0, T), p ≥ 2, then, there are equalities:

$\begin{array}{c}{D}_t^{\theta }(w_1{w}_2)=w_1(t)D_t^{\theta }{w}_2(t)+w_2(0)D_t^{\theta }{w}_1(t)\\ +\frac{\theta }{\Gamma (1-\theta )}{\mathcal{J}}_{\theta }(t;w_1,w_2),\\ {\partial }_t^{\theta }(w_1{w}_2)=w_1(t)D_t^{\theta }{w}_2(t)+w_2(0){\partial }_t^{\theta }{w}_1(t)\\ +\frac{\theta }{\Gamma (1-\theta )}{\mathcal{J}}_{\theta }(t;w_1,w_2),\end{array}$

and ${\text{D}}_t^{\theta }(w_1{w}_2),$ ${\partial }_t^{\theta }(w_1{w}_2)$ have the regularity:

$D_t^{\theta }(w_1{w}_2),{\partial }_t^{\theta }(w_1{w}_2)\in \{\begin{array}{l}\mathcal{C}([0,T]),\text{\hspace{1em}if\hspace{1em}}{D}_t^{\theta }{w}_2\in \mathcal{C}([0,T]),\hfill \\ L_p(0,T),\text{\hspace{1em}if\hspace{1em}}{D}_t^{\theta }{w}_2\in L_p(0,T).\hfill \end{array}$

(ii) For any θ₁ ≥ θ > 0 and each t ∈ [0, T], there hold

$\begin{array}{l}{I}_t^{{\theta }_1}(w_1{\partial }_t^{\theta }{w}_2)(t)=I_t^{{\theta }_1-\theta }(w_1{w}_2)(t)-w_2(0)\\ \times [I_t^{{\theta }_1-\theta }{w}_1-I_t^{{\theta }_1}(w_1{\omega }_{1-\theta })(t)]\\ -\theta I_t^{1+{\theta }_1-\theta }(\mathcal{W}(w_1)w_2)(t),\\ I_t^{{\theta }_1}(w_1{D}_t^{\theta }{w}_2)(t)=I_t^{{\theta }_1-\theta }(w_1{w}_2)(t)-w_2(0)I_t^{{\theta }_1-\theta }{w}_1\\ -\theta I_t^{1+{\theta }_1-\theta }(\mathcal{W}(w_1)w_2)(t).\end{array}$

5 A priori estimates

First, recasting step-by-step the proof of ([11], Theorem 1) and additionally exploiting [17, Theorem 3.4] and arguments leading to ([18], Theorem 4.1) in the 3BC case, we claim the following result.

Lemma 5.1. Let f(u) ≡ 0, n ≥ 2, ν, μ_j, ν_i satisfy h1, and

$p>\{\begin{array}{l}\max \{n+\frac{2}{\nu };\frac{1}{\nu -{\nu }_M};\frac{1}{\nu -{\mu }_N}\},\text{\hspace{1em} if\hspace{1em}}N\ge 1,M\ge 1,\hfill \\ \max \{n+\frac{2}{\nu };\frac{1}{\nu -{\nu }_M}\},\text{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} if\hspace{1em}}N=0,M\ge 1;\hfill \\ \max \{n+\frac{2}{\nu };\frac{1}{\nu -{\mu }_N}\},\text{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} if\hspace{1em}}N\ge 1,M=0.\hfill \end{array}$

We require that

$\begin{array}{l}g\in L_p({\Omega }_T),u_0\in W^{2-\frac{2}{p\nu },p}(\Omega ),{\phi }_1\in {ℌ}_p^{\nu (1-\frac{1}{2p}),2-\frac{1}{p}}(\partial {\Omega }_T),\\ {\phi }_2\in {ℌ}_p^{\nu (\frac{1}{2}-\frac{1}{2p}),1-\frac{1}{p}}(\partial {\Omega }_T).\end{array}$

Moreover, in the DBC case, we additionally assume

$\begin{array}{l}{u}_0(x){|}_{\partial \Omega }={\phi }_1(x,0).\end{array}$

Under assumptions h2–h5, the classical solution $u\in {\mathcal{C}}^{2+\alpha ,\frac{2+\alpha }{2}\nu }({\overline{\Omega }}_T)$ of Equations 1.1, 1.2 satisfies the estimate

$\begin{array}{l}\Vert u{\Vert }_{ℌ{}_p^{\nu ,2}({\Omega }_T)}+\Vert u{\Vert }_{{\mathcal{C}}^{\alpha ,\frac{\alpha \nu }{2}}({\overline{\Omega }}_T)}+{\displaystyle \sum _{i=1}^M\Vert D_t^{{\nu }_i}u{\Vert }_{L_p({\Omega }_T)}}\\ +{\displaystyle \sum _{j=1}^N\Vert D_t^{{\mu }_j}u{\Vert }_{L_p({\Omega }_T)}}\le \mathcal{C}\{\Vert g{\Vert }_{L_p({\Omega }_T)}+\Vert u_0{\Vert }_{W^{2-\frac{2}{p\nu },p}(\Omega )}+|\phi \left|\right\},\end{array}$

where

$|\phi |=\{\begin{array}{l}\Vert {\phi }_1{\Vert }_{ℌ{}_p^{\nu (1-\frac{1}{2p}),2-\frac{1}{p}}(\partial {\Omega }_T)}\text{\hspace{1em}\hspace{1em}in the }DBC\text{ case},\hfill \\ \Vert {\phi }_2{\Vert }_{ℌ{}_p^{\nu (\frac{1}{2}-\frac{1}{2p}),1-\frac{1}{p}}(\partial {\Omega }_T)}\text{\hspace{1em}\hspace{1em}in the }3BC\text{ case}.\hfill \end{array}$

Here, the generic constant C is independent of the right-hand sides in Equations 1.1, 1.2.

Our next result connects with a priori estimates in the fractional Hölder space to the function u satisfying the family of equations for each λ ∈ [0, 1]:

$\begin{array}{ll}{\text{D}}_{t}u-{\mathcal{L}}_{1}u-\mathcal{K}*{\mathcal{L}}_{2}u+\lambda f(u)=g(x,t)\text{\hspace{1em}in\hspace{1em}}{\Omega }_T& (5.1)\end{array}$

and homogeneous conditions Equation 1.2.

Lemma 5.2. Let assumptions of Theorem 3.1 hold, and

$\begin{array}{l}{\phi }_1,{\phi }_2,u_0\equiv 0.\end{array}$

We assume also $u\in {\mathcal{C}}^{2+\alpha ,\frac{2+\alpha }{2}\nu }({\overline{\Omega }}_T)$ be solution to Equations 5.1, 5.2. Then, for any λ ∈ [0, 1], there are the following estimates:

$\begin{array}{ll}\left|\right|u|{|}_{\mathcal{C}({\overline{\Omega }}_T)}\le \mathcal{C}\left[1+\left|\right|g|{|}_{\mathcal{C}({\overline{\Omega }}_T)}\right],& (5.2)\end{array}$

$\begin{array}{ll}\begin{array}{c}\Vert u{\Vert }_{{\mathcal{C}}^{2+\alpha ,\frac{2+\alpha }{2}{\nu }_0}({\overline{\Omega }}_T)}+{\displaystyle {\sum }_{i=1}^M\Vert D_t^{{\nu }_i}u{\Vert }_{{\mathcal{C}}^{\alpha ,\frac{\nu \alpha }{2}}({\overline{\Omega }}_T)}}\\ +{\displaystyle {\sum }_{j=1}^N\Vert D_t^{{\mu }_j}u{\Vert }_{{\mathcal{C}}^{\alpha ,\frac{\nu \alpha }{2}}({\overline{\Omega }}_T)}}\le C[1+\Vert g{\Vert }_{{\mathcal{C}}^{\alpha ,\frac{\nu \alpha }{2}}({\overline{\Omega }}_T)}].\end{array}& (5.3)\end{array}$

The positive constant C is independent of λ and the right-hand sides of Equations 5.1, 5.2 and depends only on T and the structural parameters in the model.

First of all, we notice that estimate Equation 5.3 in this claim is verified with the standard Schauder technique and by means of ([10], Theorem 4.1) and bound Equation 5.2 in this art.

We focus on the proof of Equation 5.2 if DBC holds, the case of 3BC is analyzed by collecting the similar arguments with techniques leading to ([15], Lemma 5.3). We preliminary observe that verification of Equation 5.2 in the case of the absence of ${\text{D}}_t^{{\mu }_j}({\gamma }_{j}u)$, j = 1, 2, ..., N, (i.e., N = 0) is simpler and recasts the main steps (with minor changes) in arguments related with N ≥ 1. Thus, here, we assume the presence of at least one fractional derivative ${\text{D}}_t^{{\mu }_j}({\gamma }_{j}u)$ in the operator D_tu. Then, we will exploit the following strategy. Keeping in mind assumption h3, the homogeneous initial condition and relation Equation 4.1, we rewrite D_tu in the more suitable form:

$\begin{array}{ll}\begin{array}{c}{D}_{t}u{=}_1{D}_{t}u{+}_2{D}_{t}u,{}_1{D}_{t}u={\partial }_t^{\nu }({\rho }_{0}u)+{\displaystyle {\sum }_{i=1}^M{\partial }_t^{{\nu }_i}}({\rho }_{i}u),\\ {}_2{D}_{t}u=\sum _{j=1}^N\frac{\partial }{\partial t}({\mathcal{N}}_j*({\gamma }_{j}u)),\end{array}& (5.4)\end{array}$

where

$\begin{array}{l}{\mathcal{N}}_j={\mathcal{N}}_j(t;\nu ,{\mu }_j)={\omega }_{1-\nu }(t)-{\omega }_{1-{\mu }_j}(t).\end{array}$

Appealing to (ii) in Proposition 4.2, we introduce

$\begin{array}{l}{T}_j^{*}=T^{*}({\mu }_j)>0,j=1,2,...,N,\end{array}$

such that the function ${\mathcal{N}}_j$ is strictly positive for all $t\in \left[0,T_j^{ }\right]$.

After that, for each fixed T₀:

$\begin{array}{l}0<T_0<\min \{T,T_1^{*},\mathrm{...},T_N^{*},(\nu {\mu }_1^{-1}\Gamma (1+\nu -{\mu }_1){)}^{\frac{1}{\nu -{\mu }_1}},\mathrm{...},\\ (\nu {\mu }_N^{-1}\Gamma (1+\nu -{\mu }_N){)}^{\frac{1}{\nu -{\mu }_N}}\},& (5.5)\end{array}$

we obtain the estimates

$\begin{array}{ll}\Vert u{\Vert }_{\mathcal{C}({\overline{\Omega }}_{T_0})}\le C[1+\Vert g{\Vert }_{\mathcal{C}({\overline{\Omega }}_{T_0})}]\le C[1+\Vert g{\Vert }_{\mathcal{C}({\overline{\Omega }}_T)}]& (5.6)\end{array}$

with the positive constant being independent of λ and T₀.

Then, we discuss the extension of these bounds to the interval (T₀, T] and reach the estimate Equation 5.2. It is worth noting that this step is absent in the case of N = 0, due to the proof of Equation 5.2 and consequently Equation 5.6 are carried out immediately on the entire time interval [0, T].

Step 1: Verification of Equation 5.6. Here, we focus on the obtaining of Equation 5.6 if h3(i) holds, the case of h3(ii) is analyzed with the similar arguments and is left to the interested readers.

Let $\bar{\mathcal{K}}$ be the conjugate kernel to $\mathcal{K}$, its properties are described in ([15], Proposition 4.4), in particular,

$\begin{array}{ll}||\bar{\mathcal{K}}|{|}_{{\mathcal{C}}^1(\left[0,T\right])}\le C||\mathcal{K}|{|}_{{\mathcal{C}}^1(\left[0,T\right])}(1+e^{T||\mathcal{K}|{|}_{{\mathcal{C}}^1(\left[0,T\right])}}).& (5.7)\end{array}$

Setting

$\begin{array}{l}{ℒ}_0={\displaystyle \sum _{ij=1}^n\frac{\partial }{\partial x_i}}(a_{ij}(x,t)\frac{\partial }{\partial x_j}),\overline{w}=-{ℒ}_{0}u,\\ w=-\lambda f(u)+g-D_{t}u+({ℒ}_1-{ℒ}_0)u+\mathcal{K}*({ℒ}_2-{ℒ}_0)u,\end{array}$

and exploiting [15, Proposition 4.4] and Proposition 4.1, we rewrite Equation 5.1 in more suitable form

$\begin{array}{ll}{}_1{D}_{t}u{+}_2{D}_{t}u-{ℒ}_{0}u={\displaystyle \sum _{l=1}^7{\mathcal{𝔉}}_l},& (5.8)\end{array}$

where

$\begin{array}{ll}{𝔉}_1& =-\lambda f(u)+\bar{\mathcal{K}}*f(u)+g-\bar{\mathcal{K}}*g,{𝔉}_2=({\mathcal{L}}_1-{\mathcal{L}}_0)u,\\ {𝔉}_3& =-\bar{\mathcal{K}}*({\mathcal{L}}_1-{\mathcal{L}}_0)u,{𝔉}_4=\bar{\mathcal{K}}*({\mathcal{L}}_2-{\mathcal{L}}_0)u,\\ {𝔉}_5& =\bar{\mathcal{K}}(0)({\omega }_{1-\nu }*(\rho u))+{\bar{\mathcal{K}}}^{\prime }*{\omega }_{1-\nu }*(\rho u),\\ {𝔉}_6& =-{\displaystyle \sum _{j=1}^N}\left[\bar{\mathcal{K}}(0)({\omega }_{1-{\mu }_j}*({\gamma }_{j}u))+{\bar{\mathcal{K}}}^{\prime }*{\omega }_{1-{\mu }_j}*({\gamma }_{j}u)\right],\\ {𝔉}_7& ={\displaystyle \sum _{j=1}^M}\left[\bar{\mathcal{K}}(0)({\omega }_{1-{\nu }_i}*({\rho }_{i}u))+{\bar{\mathcal{K}}}^{\prime }*{\omega }_{1-{\nu }_i}*({\rho }_{i}u)\right].\end{array}$

After that, multiplying equality (Equation 5.8) by pu^p−1 with p = 2^m, m≥1, and then integrating over Ω, we end up with the inequality (after standard technical calculations with exploiting h2)

$\begin{array}{l}{\displaystyle {\int }_{\Omega }p}{u}^{p-1}(x,\tau {)}_1{D}_{\tau }udx+{\displaystyle {\int }_{\Omega }p}{u}^{p-1}(x,\tau {)}_2{D}_{\tau }udx\\ +p(p-1){\delta }_2{\displaystyle {\int }_{\Omega }{u}^{p-2}}(x,\tau )|\nabla u(x,\tau ){|}^{2}dx\\ \le {\displaystyle \sum _{l=1}^7{\displaystyle {\int }_{\Omega }p}}{u}^{p-1}(x,\tau ){𝔉}_{l}dx.\end{array}$

It is worth noting that in the case of h3(ii), one should multiply Equation 5.8 by $p{({\rho }_{0}u)}^{p-1}$.

Computing the fractional integral $I_t^{\nu }$ of both sides in this inequality, we arrive at the bound

$\begin{array}{ll}{ℛ}_{0,1}(t)+{ℛ}_{0,2}(t)+p(p-1){\delta }_2{I}_t^{\nu }\left({\displaystyle {\int }_{\Omega }{u}^{p-2}}\right|\nabla u{|}^{2}dx)(t)\le {\displaystyle \sum _{l=1}^7{ℛ}_l}(t),& (5.9)\end{array}$

where we put

$\begin{array}{l}{ℛ}_{0,1}(t)=I_t^{\nu }\left({\displaystyle {\int }_{\Omega }p}{u^{p-1}}_1{D}_{\tau }udx\right)(t),\\ {ℛ}_{0,2}(t)=I_t^{\nu }\left({\displaystyle {\int }_{\Omega }p}{u^{p-1}}_2{D}_{\tau }udx\right)(t),\\ {ℛ}_l(t)=I_t^{\nu }\left({\displaystyle {\int }_{\Omega }{F}_l}p{u}^{p-1}dx\right)(t),l=1,\mathrm{...},7.\end{array}$

At this point, we evaluate each term ${\mathcal{R}}_l,$ ${\mathcal{R}}_{0,1}$, and ${\mathcal{R}}_{0,2}$.

• First, we notice that the terms ${\mathcal{R}}_l,$ l = 1, 2, 3, 4, are evaluated with the arguments providing the estimates of ${\mathcal{D}}_l,$ l = 1, 2, 3, 4, in ([11], Section 7.1). Thus, we immediately have

${\displaystyle \sum _{l=1}^4|}{ℛ}_l(t)|\le Cp[1+\Vert g{\Vert }_{\mathcal{C}([0,T_0])}^p]+\frac{p(p-1){\delta }_2}{2}{I}_t^{\nu }{\displaystyle {\int }_{\Omega }{u}^{p-2}}|\nabla u{|}^{2}dx,$

where the positive C is independent of λ, p, and T₀, and depends only on the structural parameters of the model.

• Coming to ${\mathcal{R}}_l$, l = 5, 6, 7, we pre-observe that ${\mathcal{R}}_6$ and ${\mathcal{R}}_7$ are evaluated with the same arguments which provide the bound of ${\mathcal{R}}_5$. Hence, here, we tackle only ${\mathcal{R}}_5$. Applying the Young inequality to the function u(x, s)u^p−1(x, τ) and then employing Proposition 4.1, estimate Equation 5.7, and assumptions h3 and h5, we get the inequality

${\displaystyle \sum _{l=5}^7|}{ℛ}_l(t)|\le Cp{I}_t^{\nu }({\displaystyle {\int }_{\Omega }|}u{|}^{p}dx)(t)$

with the positive constant C depending only on T, and the norms of γ_j, ρ_i, ρ, $\mathcal{K}$, and being independent of p, T₀, and λ.

• Now, we are left to evaluate ${\mathcal{R}}_{0,1}$ and ${\mathcal{R}}_{0,2}$. First, denoting

${\rho }_{\theta }=\{\begin{array}{l}{\rho }_0,\text{ if }\theta =0,\hfill \\ {\rho }_i,\text{ if }{\theta }_i={\nu }_i,i=1,2,\mathrm{...},M,\hfill \end{array}$

and performing technical calculations and using Propositions 4.2, 4.3, the homogeneous initial condition to u and assumption h3, we end up with the inequalities

${\displaystyle {\int }_{\Omega }p}{u}^{p-1}{\partial }_t^{\theta }({\rho }_{\theta }u)dx\ge {\displaystyle {\int }_{\Omega }{\rho }_{\theta }^{1-p}}{\partial }_t^{\theta }{({\rho }_{\theta }u)}^{p}dx,$

$\begin{array}{l}{I}_t^{\nu }({\displaystyle {\int }_{\Omega }{\rho }_{\theta }^{1-p}}{\partial }_t^{\theta }({\rho }_{\theta }^p{u}^p)dx)(t)\\ \ge \{\begin{array}{l}{\displaystyle {\int }_{\Omega }{\rho }_{\theta }}{u}^{p}dx-\nu {\displaystyle {\int }_{\Omega }{I}_t^1}(W({\rho }_{\theta }^{1-p}){\rho }_{\theta }^p{u}^p)(t)dx,\hfill \\ \text{if }\theta =\nu ,\hfill \\ I_t^{\nu -\theta }({\displaystyle {\int }_{\Omega }{\rho }_{\theta }}{u}^{p}dx)(t)-\theta {\displaystyle {\int }_{\Omega }{I}_t^{1+\nu -\theta }}(W({\rho }_{\theta }^{1-p}){\rho }_{\theta }^p{u}^p)(t)dx,\hfill \\ \text{if }\theta ={\nu }_i,i=1,\mathrm{...},M\hfill \end{array}\\ \ge \{\begin{array}{l}{\displaystyle {\int }_{\Omega }{\rho }_0}{u}^{p}dx,\text{ if }\theta =\nu ,\\ I_t^{\nu -\theta }({\displaystyle {\int }_{\Omega }{\rho }_{\theta }}{u}^{p}dx)(t),\text{ if }\theta ={\nu }_i,i=1,\mathrm{...},M.\end{array}\end{array}$

Here, to reach the last inequalities, we appeal to the definition of $\mathcal{W}$ and to assumption h3(i) (meaning the non-negativity of $\frac{\partial {\rho }_{\theta }}{\partial t}$) and taking into account the non-negativity of ${({\rho }_{\theta }u)}^p$ (since p = 2^m).

Bearing in mind these inequalities and the non-negativity of the term $I_t^{\nu -\theta }({\displaystyle {\int }_{\Omega }{\rho }_{\theta }{u}^{p}dx)(t)}$, we arrive at the desired bound

${ℛ}_{0,1}(t)\ge {\displaystyle {\int }_{\Omega }{\rho }_0}(x,t)u^p(x,t)dx.$

Concerning the term ${\mathcal{R}}_{0,2}(t)$, we will use the analogous arguments. Namely, Proposition 4.2 provides the estimate

${\displaystyle {\int }_{\Omega }p}{u}^{p-1}\frac{\partial }{\partial t}({\mathcal{N}}_j*{\gamma }_{j}u)dx\ge {\displaystyle {\int }_{\Omega }{\gamma }_j^{1-p}}\frac{\partial }{\partial t}({\mathcal{N}}_j*{({\gamma }_{j}u)}^p)dx.$

Then, collecting this bound with Proposition 4.3 arrives at inequalities:

$\begin{array}{l}{I}_t^{\nu }({\displaystyle {\int }_{\Omega }p}{u}^{p-1}\frac{\partial }{\partial t}({\mathcal{N}}_j*{\gamma }_{j}u)dx)(t)\\ \ge I_t^{\nu }({\displaystyle {\int }_{\Omega }{\gamma }_j^{1-p}}{\partial }_t^{\nu }({\gamma }_{j}u{)}^{p}dx)(t)-I_t^{\nu }({\displaystyle {\int }_{\Omega }{\gamma }_j^{1-p}}{\partial }_t^{{\mu }_j}({\gamma }_{j}u{)}^{p}dx)(t)\\ ={\displaystyle {\int }_{\Omega }{\gamma }_j}(x,t)u^p(x,t)dx-I_t^{\nu -{\mu }_j}({\displaystyle {\int }_{\Omega }{\gamma }_j}{u}^{p}dx)(t)\\ +(\nu I_t^1-{\mu }_j{I}_t^{1+\nu -{\mu }_j})({\displaystyle {\int }_{\Omega }\mathcal{W}}(-{\gamma }_j^{1-p}){\gamma }_j^p{u}^{p}dx)(t).& (5.10)\end{array}$

First, we notice that h3(i) provides the non-negativity of $\mathcal{W}(-{\gamma }_j^{1-p})$. Hence, ([12], Corollary 5.4) (where we put $w=\mathcal{W}(-{\gamma }_j^{1-p}){\gamma }_j^p{u}^p$) tells us that

$(\nu I_t^1-{\mu }_j{I}_t^{1+\nu -{\mu }_j})({\displaystyle {\int }_{\omega }\mathcal{W}}(-{\gamma }_j^{1-p}){\gamma }_j^p{u}^{p}dx)(t)\ge 0.$

After that, this bound and Equation 5.10 lead to the inequality

$\begin{array}{l}{I}_t^{\nu }({\displaystyle {\int }_{\Omega }p}{u}^{p-1}\frac{\partial }{\partial t}({\mathcal{N}}_j*{\gamma }_{j}u)dx)(t)\ge {\displaystyle {\int }_{\Omega }{\gamma }_j}(x,t)u^p(x,t)dx\\ -I_t^{\nu -{\mu }_j}\left({\displaystyle {\int }_{\Omega }{\gamma }_j}{u}^{p}dx\right)(t),\end{array}$

which in turn leads to the inequality

${ℛ}_{0,2}(t)\ge {\displaystyle {\int }_{\Omega }{\displaystyle \sum _{j=1}^N{\gamma }_j}}(x,t)u^p(x,t)dx-{\displaystyle \sum _{j=1}^N{I}_t^{\nu -{\mu }_j}}\left({\displaystyle {\int }_{\Omega }{\gamma }_j}{u}^{p}dx\right)(t).$

At last, collecting all estimates of ${\mathcal{R}}_l,$ ${\mathcal{R}}_{0,1},$ ${\mathcal{R}}_{0,2}$ with Equation 5.9, and taking into account the representation of ρ(x, t) in the case of N≥1, we arrive at the bound

$\begin{array}{l}{\displaystyle {\int }_{\Omega }\rho }(x,t)u^p(x,t)dx+\frac{p(p-1){\delta }_2}{2}{I}_t^{\nu }\left({\displaystyle {\int }_{\Omega }{u}^{p-2}}\right|\nabla u{|}^{2}dx)(t)\\ \le {\displaystyle \sum _{j=1}^N{I}_t^{\nu -{\mu }_j}}\left({\displaystyle {\int }_{\Omega }{\gamma }_j}{u}^{p}dx\right)(t)\\ +Cp(1+\Vert g{\Vert }_{\mathcal{C}({\overline{\Omega }}_{T_0})}^p)+Cp{I}_t^{\nu }\left({\displaystyle {\int }_{\Omega }{u}^p}dx\right)(t)\end{array}$

with C being independent of p, T₀, and λ.

Then, keeping in mind the restriction on ρ (see h3) to handle the first term in the left-hand side, and exploiting the easily verified relation

$|\nabla u^{p/2}{|}^2\le p(p-1)u^{p-2}|\nabla u{|}^2$

to manage the second term there, we have

$\begin{array}{ll}\begin{array}{l}{\displaystyle {\int }_{\Omega }{u}^p(x,t)dx+I_t^{\nu }}\left({\displaystyle {\int }_{\Omega }|\nabla u^{p/2}{|}^{2}dx}\right)(t)\hfill \\ \le C\underset{j}{\max }\Vert {\gamma }_j{\Vert }_{\mathcal{C}({\overline{\Omega }}_T)}{\displaystyle \sum _{j=1}^N{I}_t^{\nu -{\mu }_j}}\left({\displaystyle {\int }_{\Omega }{u}^{p}dx}\right)(t)\hfill \\ +Cp[1+\Vert g{\Vert }_{\mathcal{C}({\overline{\Omega }}_{T_0})}^p]+Cp{I}_t^{\nu }\left({\displaystyle {\int }_{\Omega }{u}^{p}dx}\right)(t).\hfill \end{array}& (5.11)\end{array}$

To handle the last term in the right-hand side, we employ the first interpolation inequality in ([15], Proposition 4.6) with $\epsilon =\frac{1}{2Cp(p-1)}$. Thus, we get

$\begin{array}{l}{\displaystyle {\int }_{\Omega }{u}^p}(x,t)dx+I_t^{\nu }({\displaystyle {\int }_{\Omega }|}\nabla u^{p/2}{|}^{2}dx)(t)\le Cp[1+\Vert g{\Vert }_{\mathcal{C}({\overline{\Omega }}_{T_0})}^p]\\ +[Cp(p-1){]}^{\frac{n+2}{2}}\Vert {\displaystyle {\int }_{\Omega }{u}^{p/2}}dx{\Vert }_{\mathcal{C}([0,T_0])}^2+C{\displaystyle \sum _{j=1}^N{I}_t^{\nu -{\mu }_j}}({\displaystyle {\int }_{\Omega }{u}^p}dx)(t).\end{array}$

Finally, taking advantage of the easily verified estimate

$\begin{array}{l}{\omega }_{\nu -{\mu }_j}(t)\le \frac{\Gamma (\nu -{\mu }_N)}{\Gamma (\nu -{\mu }_j)}{T}^{{\mu }_N-{\mu }_j}{\omega }_{\nu -{\mu }_N}(t),\\ j=1,2,...,N-1,t\in \left[0,T\right],\end{array}$

we deduce

$\begin{array}{l}{\displaystyle {\int }_{\Omega }{u}^p}(x,t)dx\text{ }\le Cp[1+\Vert g{\Vert }_{\mathcal{C}({\overline{\Omega }}_{T_0})}^p]\\ +[Cp(p-1){]}^{\frac{n+2}{2}}{\Vert {\displaystyle {\int }_{\Omega }{u}^{p/2}}dx\Vert }_{C([0,T_0])}^2\\ +C^{*}{I}_t^{\nu -{\mu }_N}\left({\displaystyle {\int }_{\Omega }{u}^p}dx\right)(t),& (5.12)\end{array}$

where

$C^{*}=C\left[1+{\displaystyle \sum _{j=1}^{N-1}\frac{\Gamma (\nu -{\mu }_N)}{\Gamma (\nu -{\mu }_j)}{T}^{{\mu }_N-{\mu }_{j]}}}\right]$

being independent of T₀, p, and λ.

To control the last term in the right-hand side, we apply the Gronwall-type inequality [15, Proposition 4.3] and then use formula (3.7.43) in [19]. Thus, we have

$\begin{array}{l}{C}^{*}{I}_t^{\nu -{\mu }_N}\left({\displaystyle {\int }_{\Omega }{u}^p}dx\right)(t)\le C^{*}A{I}_t^{\nu -{\mu }_N}(E_{\nu -{\mu }_N}(C^{*}{t}^{\nu -{\mu }_N}))(t)\\ =A[E_{\nu -{\mu }_N}(C^{*}{t}^{\nu -{\mu }_N})-1]\\ \le A[E_{\nu -{\mu }_N}(C^{*}{T}^{\nu -{\mu }_N})-1]\\ \text{ for all }t\in [0,T],\end{array}$

where we put

$A=Cp[1+\Vert g{\Vert }_{\mathcal{C}({\overline{\Omega }}_{T_0})}^p]+{[Cp(p-1)]}^{\frac{n+2}{2}}{\Vert {\displaystyle {\int }_{\Omega }{u}^{p/2}}dx\Vert }_{\mathcal{C}([0,T_0])}^2,$

and $E_{\theta }(t)=\sum _{k=0}^{+\infty }\frac{z^k}{\Gamma (k\theta +1)}$ is the classical Mittag-Leffler function of the order θ.

Taking into account this estimate to evaluate the last term in the right-hand side of Equation 5.12, we achieve

${\displaystyle {\int }_{\Omega }{u}^p(x,t)dx\le A{E}_{\nu -{\mu }_N}(C^{*}{T}^{\nu -{\mu }_N}).}$

In fine, denoting

$\begin{array}{l}ℬ=4C{E}_{\nu -{\mu }_N}(C^{*}{T}^{\nu -{\mu }_N}),A_m=\underset{t\in [0,T_0]}{\sup }({\displaystyle {\int }_{\Omega }{u}^p}dx{)}^{1/p}\\ \text{ with }p=2^m,\end{array}$

we derive the bound

$\begin{array}{ll}{\mathcal{A}}_m\le {\mathcal{B}}^{m{2}^{-m}}[1+||g|{|}_{\mathcal{C}({\bar{\Omega }}_{T_0})}]+{\mathcal{B}}^{mn{2}^{-m}}{\mathcal{A}}_{m-1}.& (5.13)\end{array}$

Then, two possibilities occur:

(i) either $max\left\{{\mathcal{A}}_{m-1},1+\left|\right|g|{|}_{\mathcal{C}({\bar{\Omega }}_{T_0})}\right\}=1+\left|\right|g|{|}_{\mathcal{C}({\bar{\Omega }}_{T_0})}$,

(ii) or $max\left\{{\mathcal{A}}_{m-1},1+\left|\right|g|{|}_{\mathcal{C}({\bar{\Omega }}_{T_0})}\right\}={\mathcal{A}}_{m-1}$.

Clearly, in the case of (i), passing to the limit as m → +∞ in Equation 5.13, we end up with the desired estimate for t∈[0, T₀].

If (ii) holds, then the standard technical calculations arrive at the inequality

$\begin{array}{l}{\mathcal{A}}_m\le [B^{m{2}^{-m}}+B^{nm{2}^{-m}}]{\mathcal{A}}_{m-1}<C{\displaystyle \prod _{k=1}^m{[B+B^n]}^{k{2}^{-k}}{\mathcal{A}}_1}\\ <C\exp \left\{|\ln [B+B^n]|{\displaystyle \sum _{k=1}^{+\infty }\frac{kn}{2^k}}\right\}{\mathcal{A}}_1.\end{array}$

Letting m → +∞ in this estimates and bearing in mind the convergence of the series, we have

$||u|{|}_{\mathcal{C}({\bar{\Omega }}_{T_0})}\le C{\mathcal{A}}_1,$

where the positive constant C is independent of T₀ and λ.

Finally, to manage the term ${\mathcal{A}}_1$, we first put p = 2 in Equation 5.11 and then apply Gronwall inequality [15, Proposition 4.3], where we set $k={\omega }_{\nu }(t)+C\underset{j}{\max }\left|\right|{\gamma }_j\left|\right|{\mathcal{C}}_{({\overline{\Omega }}_T}{\displaystyle {\sum }_{j=1}^N{\omega }_{\nu -{\mu }_j}(t)}$. Thus, we end up with bound Equation 5.6 and as a consequence with Equation 5.3 where T = T₀.

Step 2: Extension of Equation 5.6 to the whole time interval. Actually, we only need in the technique which allows us to extend Equation 5.6 to the interval [T₀, 3T₀/2]. Then, repeating this procedure a finite number of times, we exhaust the entire [T₀, T] and hence complete the proof of Equation 5.2.

Denoting

$\Phi (x,t)=\{\begin{array}{l}-\lambda f(u)+g(x,t),\text{ if }(x,t)\in {\overline{\Omega }}_{T_0/2},\hfill \\ [-\lambda f(u)+g(x,t)]{|}_{t=T_0/2},\text{ if }x\in \overline{\Omega },t>T_0/2,\hfill \end{array}$

we designate 𝔘(x, t) as a solution to the linear problem

$\begin{array}{ll}\{\begin{array}{l}{D}_t\mathcal{U}-{ℒ}_1\mathcal{U}-\mathcal{K}*{ℒ}_2\mathcal{U}=\Phi (x,t)\text{ in }{\Omega }_{3T_0/2},\hfill \\ \mathcal{U}(x,0)=0\text{ in }\overline{\Omega },\hfill \\ \mathcal{U}(x,t)=0\text{ on }\partial {\Omega }_{3T_0/2}.\hfill \end{array}& (5.14)\end{array}$

Thanks to Equations 5.3, 5.6 (with T = T₀) and assumptions h6, h7.I, we get

$\begin{array}{l}\Vert \Phi {\Vert }_{{\mathcal{C}}^{\alpha ,\frac{\alpha {\nu }_0}{2}}({\overline{\Omega }}_{3T_0/2})}\\ \le C[\Vert u{\Vert }_{\mathcal{C}({\overline{\Omega }}_{T_0})}+1+\Vert g{\Vert }_{{\mathcal{C}}^{\alpha ,\frac{\alpha \nu }{2}}({\overline{\Omega }}_{T_0})}]\le C[1+\Vert g{\Vert }_{{\mathcal{C}}^{\alpha ,\frac{\alpha \nu }{2}}({\overline{\Omega }}_{T_0})}],\\ \Vert \Phi {\Vert }_{\mathcal{C}({\overline{\Omega }}_{3T_0/2})}\le C[1+\Vert g{\Vert }_{\mathcal{C}({\overline{\Omega }}_{T_0})}],\\ \Phi (x,0)=0\text{ if }x\in \partial \Omega ,& (5.15)\end{array}$

where the positive value C is independent of T₀, λ and the right-hand side of Equation 5.14.

Keeping in mind these properties of Φ, we can apply [10, Theorem 4.1] to Equation 5.14 and obtain the unique classical solution 𝔘 satisfying the following relations:

$\begin{array}{l}\Vert \mathcal{𝔘}{\Vert }_{{\mathcal{C}}^{2+\alpha ,\frac{2+\alpha }{2}\nu }({\overline{\Omega }}_{3T_0/2})}+{\displaystyle \sum _{i=1}^N\Vert D_t^{{\nu }_i}\mathcal{𝔘}{\Vert }_{{\mathcal{C}}^{\alpha ,\frac{\alpha \nu }{2}}({\overline{\Omega }}_{3T_0/2})}}\\ +{\displaystyle \sum _{j=1}^M\Vert D_t^{{\mu }_j}\mathcal{𝔘}{\Vert }_{{\mathcal{C}}^{\alpha ,\frac{\alpha \nu }{2}}({\overline{\Omega }}_{3T_0/2})}}\\ \le C[1+\Vert g{\Vert }_{{\mathcal{C}}^{\alpha ,\frac{\alpha \nu }{2}}({\overline{\Omega }}_{T_0})}],\\ \Vert \mathcal{𝔘}{\Vert }_{\mathcal{C}({\overline{\Omega }}_{3T_0/2})}\le C[1+\Vert g{\Vert }_{\mathcal{C}({\overline{\Omega }}_{T_0})}],\\ \mathcal{𝔘}(x,t)=u(x,t)\text{ if }(x,t)\in {\overline{\Omega }}_{T_0/2}.\end{array}$

In fine, we introduce new unknown function

$𝔳(x,t)=u(x,t)-𝔘(x,t)$

solving the problem

$\begin{array}{ll}\{\begin{array}{l}{D}_{t}v-{ℒ}_{1}v-\mathcal{K}*{ℒ}_{2}v=-\lambda f^{\star }(v)+g^{\star }(x,t)\text{ in }{\Omega }_{3T_0/2},\hfill \\ v(x,0)=0\text{ in }\overline{\Omega },\hfill \\ v(x,t)=0\text{ on }\partial {\Omega }_{3T_0/2}.\hfill \end{array}& (5.16)\end{array}$

Here, we set

$f^{\star }(𝔳)=f(𝔳+𝔘),g^{\star }(x,t)=g(x,t)-\Phi (x,t).$

By virtue of Equation 5.15 and representation of the right-hand sides in Equation 5.16, we deduce that f^*(𝔳) has all properties of f(u), and

$\begin{array}{l}{g}^{*}-\lambda f^{*}=0\text{ if }x\in \overline{\Omega },t\in [0,T_0/2],\\ \Vert g^{*}{\Vert }_{{\mathcal{C}}^{\alpha ,\frac{\alpha \nu }{2}}({\overline{\Omega }}_{3T_0/2})}\le C[1+\Vert g{\Vert }_{{\mathcal{C}}^{\alpha ,\frac{\alpha \nu }{2}}({\overline{\Omega }}_{T_0})}],\\ \Vert g^{*}{\Vert }_{\mathcal{C}({\overline{\Omega }}_{3T_0/2})}\le C[1+\Vert g{\Vert }_{\mathcal{C}({\overline{\Omega }}_{T_0})}],\end{array}$

where the constant C is independent of λ and T₀.

Finally, introducing the new time-variable

$\sigma =t-\frac{T_0}{2},\sigma \in \left[-\frac{T_0}{2},T_0\right],$

and repeating arguments of the end of Section 6.3 in [10], we arrive at the problem

$\begin{array}{ll}\{\begin{array}{l}{\overline{D}}_{\sigma }\overline{v}-{\overline{ℒ}}_1\overline{v}-\mathcal{K}*{\overline{ℒ}}_2\overline{v}=-\lambda {\overline{f}}^{\star }(\overline{v})+{\overline{g}}^{\star }\text{ in }{\Omega }_{T_0},\hfill \\ \overline{v}(x,0)=0\text{ in }\overline{\Omega },\hfill \\ \overline{v}(x,\sigma )=0\text{ on }\partial {\Omega }_{T_0},\hfill \end{array}& (5.17)\end{array}$

besides,

$\overline{v}(x,\sigma )=0\text{ if }\sigma \in [-\frac{T_0}{2},0],x\in \overline{\Omega }.$

Here, we put

$\begin{array}{l}\overline{v}(x,\sigma )=v(x,\sigma +T_0/2),{\overline{g}}^{*}(x,\sigma )=g^{*}(x,\sigma +T_0/2),\\ {\overline{f}}^{*}(\overline{v})=f^{*}(v){|}_{t=\sigma +T_0/2},\end{array}$

and we call ${\overline{\mathcal{L}}}_k,$ ${\overline{\text{D}}}_{\sigma }$ the operators ${\mathcal{L}}_k$ and D_σ, respectively, with the bar coefficients. It is easy to check that the coefficients of these operators and the functions ḡ^* and ${\bar{f}}^{*}$ meet the requirements of Lemma 5.2. Then, arguing as Step 1, we end up with estimates Equations 5.2, 5.3, 5.6 to the function 𝔳. Collecting the obtained results with the properties of the function 𝔘, we extend the desired estimates to the whole segment [0, 3T₀/2]. This completes the proof of Lemma 5.2

Remark 5.1. Collecting estimate Equation 5.2 with Lemma 5.1 provides the following a priori estimate to solution of Equation 5.1 satisfying homogeneous boundary and initial conditions:

$\begin{array}{l}\Vert u{\Vert }_{ℌ{}_p^{\nu ,2}({\Omega }_T)}+\Vert u{\Vert }_{{\mathcal{C}}^{\alpha ,\frac{\alpha \nu }{2}}({\overline{\Omega }}_T)}+{\displaystyle \sum _{i=1}^M\Vert D_t^{{\nu }_i}u{\Vert }_{L_p({\Omega }_T)}}\\ +{\displaystyle \sum _{j=1}^N\Vert D_t^{{\mu }_j}u{\Vert }_{L_p({\Omega }_T)}}\\ \le C[1+\Vert g{\Vert }_{L_p({\Omega }_T)}+\Vert g{\Vert }_{\mathcal{C}({\overline{\Omega }}_T)}]\end{array}$

with C being independent of λ.

6 Proof of Theorems 3.1, 3.2

Here, we will exploit the continuation approach based on the a priori estimates in the fractional Hölder spaces. It is worth noting that this technique has been utilized in [11] to prove the well-posedness of Equations 1.1, 1.2 with two-term fractional derivatives in the operator Equation 1.3 in the DBC case. Hence, in our arguments, we focus on only main difficulties connected with multi-term fractional derivatives in Equation 1.3.

Concerning the proof of Theorem 3.2, we will exploit the technique leading to Theorem 4.4. in [12]. This approach includes a priori estimates of Equations 1.1, 1.2 in the fractional Sobolev spaces and the construction of the corresponding solutions via consideration of approximated problems.

6.1 Conclusion of the proof of Theorem 3.1

First, we prove Theorem 3.1 in the case of homogeneous boundary and initial conditions and then we remove this restriction.

To this end, we rely on the so-called continuation arguments. For λ∈[0, 1], we consider the family of problem

$\begin{array}{ll}\{\begin{array}{l}{D}_{t}u-{ℒ}_{1}u-\mathcal{K}*{ℒ}_{2}u+\lambda f(u)=g(x,t)\text{ in }{\Omega }_T,\hfill \\ u(x,0)=0\text{ in }\overline{\Omega },\hfill \\ u(x,t)=0\text{ or }ℳu+{\mathcal{K}}_1*ℳu-c_{0}u=0\text{ on }\partial {\Omega }_T.\hfill \end{array}& (6.1)\end{array}$

Denoting Λ as the set of those λ for which Equation 6.1 is solvable on [0, T]. Obviously, if λ = 0, then Equation 6.1 transforms to the linear problem analyzed in [10]. Hence, assumptions h1–h6 allow us to apply Theorem 4.1 and Remark 4.4 from [10] and obtain the global classical solvability. Thus, 0∈Λ. Then, we are left to examine if the set Λ is open and closed at the same time. To this end, exploiting Lemmas 5.1, 5.2 (in particular, the estimate of $\left|\right|u|{|}_{{\mathcal{C}}^{\alpha ,\alpha \nu /2}({\bar{\Omega }}_T)}$ via $\left|\right|g|{|}_{\mathcal{C}({\bar{\Omega }}_T)}$) and recasting step-by-step the arguments of ([15], Section 5.2), we complete the proof of Theorem 3.1 in the case of homogeneous initial and boundary conditions.

To remove this restriction, we consider the following linear problem with the unknown function w = w(x, t)

$\{\begin{array}{l}{D}_{t}w-{ℒ}_{1}w-\mathcal{K}*{ℒ}_{2}w=g(x,t)-f(u_0)\text{ in }{\Omega }_T,\hfill \\ w(x,0)=u_0(x)\text{ in }\overline{\Omega },\hfill \\ w(x,t)={\phi }_1(x,t)\text{ or }ℳw+{\mathcal{K}}_1*ℳw-c_{0}w={\phi }_2(x,t)\hfill \\ \text{on }\partial {\Omega }_T.\hfill \end{array}$

Applying ([15], Remark 3.1) and ([10], Remark 4.4) arrives at the one-valued classical solvability of this linear problem and, besides, at the bound

$\begin{array}{l}||w|{|}_{{\mathcal{C}}^{2+\alpha ,\frac{2+\alpha }{2}\nu }({\bar{\Omega }}_T)}+{\displaystyle \sum _{i=1}^M}||{\text{D}}_t^{{\nu }_i}w|{|}_{{\mathcal{C}}^{\alpha ,\frac{\nu \alpha }{2}}({\bar{\Omega }}_T)}+{\displaystyle \sum _{j=1}^N}||{\text{D}}_t^{{\mu }_j}w|{|}_{{\mathcal{C}}^{\alpha ,\frac{\nu \alpha }{2}}({\bar{\Omega }}_T)}\\ \le CG(u_0,g,\phi ),\end{array}$

where

$\begin{array}{l}G(u_0,g,\phi )=1+\Vert g{\Vert }_{{\mathcal{C}}^{\alpha ,\frac{\nu \alpha }{2}}({\overline{\Omega }}_T)}+\Vert u_0{\Vert }_{{\mathcal{C}}^{2+\alpha }(\overline{\Omega })}+|\phi {|}_{\mathcal{C}},\\ |\phi {|}_{\mathcal{C}}=\{\begin{array}{l}\Vert {\phi }_1{\Vert }_{{\mathcal{C}}^{2+\alpha ,\frac{2+\alpha }{2}\nu }(\partial {\Omega }_T)}\text{ in }DBC\text{ case},\hfill \\ \Vert {\phi }_2{\Vert }_{{\mathcal{C}}^{1+\alpha ,\frac{1+\alpha }{2}\nu }(\partial {\Omega }_T)}\text{ in }3BC\text{ case}.\hfill \end{array}\end{array}$

Here, we exploited assumption h7 and ([15], Remark 3.1) to handle the term $||f|{|}_{{\mathcal{C}}^{\alpha ,\frac{\nu \alpha }{2}}({\bar{\Omega }}_T)}$.

After that, we look for a solution to the original problem Equations 1.1, 1.2 in the form

$u(x,t)=w(x,t)+W(x,t),$

where the new unknown function W solves the problem Equation 6.1 with λ = 1 and the new right-hand sides:

$\bar{f}(W)=f(W+w)-f(w),ḡ=f(u_0)-f(w).$

Remark 6.1. Assumption h4 and the estimate of w provide the inequality

$\Vert \overline{g}{\Vert }_{{\mathcal{C}}^{\alpha ,\frac{\nu \alpha }{2}}({\overline{\Omega }}_T)}\le CG(u_0,g,\phi ).$

In addition, the function $\bar{f}(W)$ satisfies assumption h7 with the constant depending only on L or L_i and $\mathcal{G}(u_0,g,\phi )$. Moreover, the straightforward calculations and the definition of w arrive at the relations

$\overline{g}(x,0)=0\text{ for each }x\in \overline{\Omega },\overline{f}(0)=0\text{ for each }(x,t)\in {\overline{\Omega }}_T.$

The last equalities in Remark 6.1 tell us that the consistency conditions in the non-linear problem for the function W are satisfied. In summary, we reduce problem Equation 1.1, 1.2 to Equation 6.1 with the right-hand sides satisfying the assumptions of Theorem 3.1. Hence, this completes the proof of this theorem in the general case.

6.2 Proof of Theorem 3.2

Actually, the verification of Theorem 3.2 is a simple consequence of Theorem 3.1 and a priori estimates obtained in Section 5 and repeats the arguments leading to ([12], Theorem 4.4). Indeed, thanks to Theorem 3.1 in the case of homogeneous initial and boundary conditions in Equation 1.2, we construct an approximate solution u_n. Then, exploiting uniform estimates in Lemma 5.1 and Remark 5.1 and passing to the limit via standard arguments, we obtain a strong solution to Equations 1.1, 1.2 satisfying the regularity stated in Theorem 3.2. Finally, to reach the uniqueness of this solution, we assume the existence of two solutions u₁ and u₂ satisfying Equations 1.1, 1.2 with the same right-hand sides. Clearly, the difference ū = u₁−u₂ solves the problem Equation 6.1 with λ = 1, g = 0 and f(ū) = f(u₁)−f(u₂), where

$|f(\overline{u})|\le C|u_1-u_2|,C=\{\begin{array}{l}L,\text{ ifh7}.\text{Iholds},\hfill \\ |f^{\prime }(\xi )|,\text{ if h7.II holds},\hfill \end{array}$

where ξ is a middle point lying between u₁ and u₂.

Finally, recasting the arguments leading to the estimate Equation 5.2, we obtain the equality

$\overline{u}=0,(x,t)\in {\overline{\Omega }}_T,$

which finishes the verification of Theorem 3.2.

7 Conclusion

In this study, we propose a technique to study the well-posedness (for each fixed T) of initial-boundary value problems to semilinear multi-term time-fractional diffusion equations with memory. The particular case of the problems analyzed models the oxygen transport through capillaries [6]. The introduction of fractional calculus in the model of the evolution of the oxygen density is well-presented with some interesting details. Our approach is particularly efficient when the multi-term derivatives can be represented in the form $\frac{\partial }{\partial t}(\mathcal{N} \rho u)$ with a some non-positive kernel $\mathcal{N}$ and given coefficient ρ = ρ(x, t).

Our analytical technique and ideas can be incorporated to study the corresponding inverse problems concerning the reconstruction of unknown parameters (e.g., the time lag in concentration of oxygen along capillaries; the order of oxygen subdiffusion; and so on). Moreover, our investigation can be employed to analyze the corresponding initial-boundary value problems to fully non-linear equations containing a term $\frac{\partial }{\partial t}(\mathcal{N} f(u))$ and to the equations with degenerate coefficients in the fractional operator. These issues will be addressed with a possible further research.

Data availability statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

Author contributions

NV: Conceptualization, Investigation, Methodology, Supervision, Writing – original draft, Writing – review & editing.

Funding

The author(s) declare that no financial support was received for the research, authorship, and/or publication of this article.

Conflict of interest

The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

1. Heymans N, Bauwens JC. Fractal rheological models and fractional differential equations for viscoelastic behavior. Rheol Acta. (1994) 83:210–9. doi: 10.1007/BF00437306

Crossref Full Text | Google Scholar

2. Matsuzaki T, Nakagawa M. A chaos neuron model with fractional differential equation. J Phys Soc Japan. (2003) 72:2678–84. doi: 10.1143/JPSJ.72.2678

Crossref Full Text | Google Scholar

3. Picozzi S, West BJ. Fractional Langevan model of memory in financial markets. Phys Rev E. (2002) 66:46–118. doi: 10.1103/PhysRevE.66.046118

PubMed Abstract | Crossref Full Text | Google Scholar

4. Fritz M, Kuttler M, Rajendran MZ, Wohlmuth B, Scarabosio L. On a subdiffusive tumor growth model with fractional time derivative. IMA J Appl Math. (2021) 86:688–729. doi: 10.1093/imamat/hxab009

Crossref Full Text | Google Scholar

5. Morales-Delgado VF, Gómez-Aguilar JF, Saad KM, Khan MA, Agarwal P. Analytic solution for oxygen diffusion from capillary to tissues involving external force effects: a fractional calculus approach. Physica A. (2019) 523:48–65. doi: 10.1016/j.physa.2019.02.018

Crossref Full Text | Google Scholar

6. Srivastava V, Rai KN. A multi-term fractional diffusion equation for oxygen delivery through a capillary to tissues. Math Comput Model. (2010) 51:616–24. doi: 10.1016/j.mcm.2009.11.002

Crossref Full Text | Google Scholar

7. Daftardar-Gejji V, Bhalekar S. Boundary value problems for multi-term fractional differential equations. J Math Anal Appl. (2008) 345:754–65. doi: 10.1016/j.jmaa.2008.04.065

Crossref Full Text | Google Scholar

8. Kassymov A, Tokmagambetov N, Torebek BT. Multi-term time-fractional diffusion equation and system: mild solutions and critical exponents. Publ Math Debrecen. (2022) 100:295–321. doi: 10.5486/PMD.2022.8922

Crossref Full Text | Google Scholar

9. Janno J, Kinash N. Reconstruction of an order of derivative and a source term in a fractional diffusion equation from final measurements. Inverse Probl. (2018) 34:02507. doi: 10.1088/1361-6420/aaa0f0

Crossref Full Text | Google Scholar

10. Pata V, Siryk S, Vasylyeva N. Multi-term fractional linear equations modeling oxygen subdiffusion through capillaries. arXiv preprint arXiv:2210.05009 (2022).

Google Scholar

11. Vasylyeva N. Cauchy-Dirichlet problem to semilinear multi-term fractional differential equations. Fractal Fract. (2023) 7:249. doi: 10.3390/fractalfract7030249

Crossref Full Text | Google Scholar

12. Siryk SV, Vasylyeva N. Initial-boundary value problems to semilinear multi-term fractional differential equations. Comm Pure Appl Anal. (2023) 22:2321–64. doi: 10.3934/cpaa.2023068

Crossref Full Text | Google Scholar

13. Giorgi C, Pata V, Marzocchi A. Asymptotic behavior of a semilinear problem in heat conduction with memory. NoDEA. (1998) 5:333–54. doi: 10.1007/s000300050049

Crossref Full Text | Google Scholar

14. Adams RA, Fournier JJF. Sobolev Spaces. 2nd edition. New York: Academic Press (2003).

Google Scholar

15. Krasnoschok M, Pata V, Vasylyeva N. Semilinear subdiffusion with memory in multidimensional domains. Mathematische Nachrichten. (2019) 292:1490–513. doi: 10.1002/mana.201700405

Crossref Full Text | Google Scholar

16. Krasnoschok M, Pata V, Vasylyeva N. Semilinear subdiffusion with memory in the one-dimensional case. Nonlinear Anal. (2017) 165:1–17. doi: 10.1016/j.na.2017.09.004

Crossref Full Text | Google Scholar

17. Zacher R. Quasilinear parabolic integro-differential equations with nonlinear boundary conditions. Differ Integral Equa. (2006) 19:1129–56. doi: 10.57262/die/1356050312

Crossref Full Text | Google Scholar

18. Zacher R. Global strong solvability of a quasilinear subdiffusion problem. J Evol Equ. (2012) 12:813–31. doi: 10.1007/s00028-012-0156-0

Crossref Full Text | Google Scholar

19. Gorenflo R, Kilbas A, Mainardi F, Rogosin S. Mittag-Leffler Functions: Related Topics ad Applications. Berlin-Heidelberg: Springer. (2014). doi: 10.1007/978-3-662-43930-2

Crossref Full Text | Google Scholar