# Multi-term fractional linear equations
modeling oxygen subdiffusion through capillaries

Vittorino Pata, Sergii V. Siryk and Nataliya Vasylyeva  Dipartimento di
Matematica, Politecnico di Milano
Via Bonardi 9, 20133 Milano, Italy<EMAIL_ADDRESS>CONCEPT Lab,
Istituto Italiano di Tecnologia
Via Morego 30, 16163, Genova, Italy<EMAIL_ADDRESS>Institute of Applied
Mathematics and Mechanics of NASU
G.Batyuka st. 19, 84100 Sloviansk, Ukraine;
and Dipartimento di Matematica, Politecnico di Milano
Via Bonardi 9, 20133 Milano, Italy  nataliy<EMAIL_ADDRESS>
###### Abstract.

For $0<\nu_{2}<\nu_{1}\leq 1$, we analyze a linear integro-differential
equation on the space-time cylinder $\Omega\times(0,T)$ in the unknown
$u=u(x,t)$

$\mathbf{D}_{t}^{\nu_{1}}(\varrho_{1}u)-\mathbf{D}_{t}^{\nu_{2}}(\varrho_{2}u)-\mathcal{L}_{1}u-\mathcal{K}*\mathcal{L}_{2}u=f$

where $\mathbf{D}_{t}^{\nu_{i}}$ are the Caputo fractional derivatives,
$\varrho_{i}=\varrho_{i}(x,t)$ with $\varrho_{1}\geq\mu_{0}>0$,
$\mathcal{L}_{i}$ are uniform elliptic operators with time-dependent smooth
coefficients, $\mathcal{K}$ is a summable convolution kernel, and $f$ is an
external force. Particular cases of this equation are the recently proposed
advanced models of oxygen transport through capillaries. Under suitable
conditions on the given data, the global classical solvability of the
associated initial-boundary value problems is addressed. To this end, a
special technique is needed, adapting the concept of a regularizer from the
theory of parabolic equations. This allows us to remove the usual assumption
about the nonnegativity of the kernel representing fractional derivatives. The
problem is also investigated from the numerical point of view.

###### Key words and phrases:

Oxygen subdiffusion, Caputo derivative, a priori estimates, regularizer,
classical solvability, finite-difference schemes

###### 2000 Mathematics Subject Classification:

Primary 35R11, 35C15, 65M06; Secondary 45N05, 26A33, 35B30

##  1\. Introduction

Oxygen transport is a complex phenomenon including chemical reactions with
hemoglobin, convective transport in red blood cells, diffusion and metabolic
consumption [38, 39]. Convective oxygen in blood depends on active energy
consuming processes generating flow in the circulation. Diffusion transport
refers to the passive movement of oxygen down its concentration gradient
across tissue barriers, including the alveolar-capillary membrane, and across
the extracellular matrix between the tissue capillaries and individual cells
to mitochondria. The amount of diffusive oxygen movement depends on the oxygen
tension gradient and the diffusion distance, which is related to the tissue
capillary density. The greater is the difference between capillary and
cellular oxygen concentration and the shorter is the distance, the faster is
the rate of diffusion [29, 50]. In abnormal body circulation, cells closer to
the capillary at the venous end begin to suffer from hypoxia when perfusion
levels drop to critically low values [6, 50]. The mechanisms controlling
oxygen distribution involving a series of convective and diffusive processes
are not yet completely understood [39].

There are actually some methods to measure the oxygen level, such as two-
photon phosphorescence lifetime microscopy, that can be applied in vivo [30].
However, the existing techniques can hardly offer a complete spatial-temporal
picture of the oxygen field on microscopic scales. Thus, analytic
studying/theoretical modeling [6, 27, 38, 39] and numerical simulation [8, 33,
42] are widely utilized in evaluating oxygen level and angiogenesis research.
The classical Krogh cylinder model roughly describes the oxygen transport from
blood vessels to tissues [27]. In particular, Krogh proposed that oxygen is
transported in the tissue by passive diffusion driven by gradients of oxygen
tension, and gave a simple geometrical model of the elementary tissue unit
supplied by a single capillary. Coupled models for oxygen delivery, even in
presence of a relatively complex vessel network structure, and a detailed
description of the blood flow in the vessel network were also proposed in [7,
39]. Go [6] used a mathematical model for oxygen delivery through capillaries
where the longitudinal diffusion of solute is neglected, and the diffusion and
the consumption rate of oxygen are assumed to be the same everywhere, which is
not the case in real situations [29]. The further paper [46] introduces a new
advanced mathematical model for oxygen delivery through a capillary to tissue
in both (transverse and longitudinal) directions. In this work, conveying
oxygen from the capillary to the surrounding tissue is described by means of a
subdiffusion equation containing two fractional derivatives in time, that is,

$\mathbf{D}_{t}^{\nu_{1}}\mathfrak{C}-\tau\mathbf{D}_{t}^{\nu_{2}}\mathfrak{C}={\rm
div}(\varrho\nabla\mathfrak{C})-k,$

with $0<\nu_{2}<\nu_{1}\leq 1$. The equation can also exhibit extra terms
accounting for the presence of external forces, even in convolution form [37].
Here, $\mathfrak{C}$ is a function of space and time, representing the
concentration of oxygen, $\tau$ is the time lag in concentration of oxygen
along the capillary, $k$ is the rate of consumption per volume of tissue, and
$\varrho$ is the diffusion coefficient of oxygen, which possibly dependent on
$\mathfrak{C}$. In particular, the term
$\mathbf{D}_{t}^{\nu_{1}}\mathfrak{C}-\tau\mathbf{D}_{t}^{\nu_{2}}\mathfrak{C}$
details the net diffusion of oxygen to all tissues. In the equation, the
symbol $\mathbf{D}_{t}^{\theta}$ stands for the usual Caputo fractional
derivative of order $\theta\in(0,1)$ with respect to time, defined as

$\mathbf{D}_{t}^{\theta}\mathfrak{h}(t)=\frac{1}{\Gamma(1-\theta)}\frac{d}{dt}\int\limits_{0}^{t}\frac{\mathfrak{h}(s)-\mathfrak{h}(0)}{(t-s)^{\theta}}ds=\frac{1}{\Gamma(1-\theta)}\int\limits_{0}^{t}(t-s)^{-\theta}\frac{d\mathfrak{h}}{ds}(s)ds,$

where $\Gamma$ is the Euler Gamma function, and the latter equality holds if
$\mathfrak{h}$ is an absolutely continuous function. In the limit cases
$\theta=0$ and $\theta=1$, the Caputo fractional derivatives of
$\mathfrak{h}(t)$ boil down to $[\mathfrak{h}(t)-\mathfrak{h}(0)]$ and
$\frac{d\mathfrak{h}}{dt}(t)$, respectively.

In this paper, motivated by the discussion above, we focus on the analytical
and the numerical study of initial-boundary value problems for evolution
equations with multi-term fractional derivatives. Let
$\Omega\subset\mathbb{R}^{n},$ with $n\geq 1$, be a bounded domain with smooth
boundary $\partial\Omega$, and let $T>0$ be an arbitrarily fixed final time.
We denote

$\Omega_{T}=\Omega\times(0,T)\qquad\text{and}\qquad\partial\Omega_{T}=\partial\Omega\times[0,T].$

For $0<\nu_{2}<\nu_{1}\leq 1$, we consider the following non-autonomous multi-
term subdiffusion equation with memory terms in the unknown function
$u=u(x,t):\Omega_{T}\to\mathbb{R},$

$\mathbf{D}_{t}^{\nu_{1}}(\varrho_{1}u)-\mathbf{D}_{t}^{\nu_{2}}(\varrho_{2}u)-\mathcal{L}_{1}u-\mathcal{K}*\mathcal{L}_{2}u=f,$
(1.1)

where the symbol $*$ stands for the usual time-convolution product

$(\mathfrak{h}_{1}*\mathfrak{h}_{2})(t)=\int\limits_{0}^{t}\mathfrak{h}_{1}(t-s)\mathfrak{h}_{2}(s)ds.$

Here, $\varrho_{i}=\varrho_{i}(x,t)$ and $f=f(x,t)$ are given functions,
$\mathcal{K}$ is a summable convolution kernel, and $\mathcal{L}_{i}$ are
linear elliptic operators of the second order with time-dependent
coefficients, whose precise form will be given in Section 3, where we detail
the general assumptions of our problem. The equation is supplemented with the
initial condition

$u(x,0)=u_{0}(x)\quad\text{in}\quad\bar{\Omega},$ (1.2)

and subject to the one of the following boundary conditions on
$\partial\Omega_{T}$:

  * (i)

Dirichlet boundary condition (DBC)

$u(x,t)=\psi_{1},$ (1.3)

  * (ii)

Boundary condition of the third kind (3BC)

$\mathcal{M}_{1}u+\mathcal{K}_{0}*\mathcal{M}_{2}u=\psi_{2},$ (1.4)

  * (iii)

Fractional dynamic boundary condition (FDBC)

$\mathbf{D}_{t}^{\nu_{1}}(\varrho_{1}u)-\mathbf{D}_{t}^{\nu_{2}}(\varrho_{2}u)-\mathcal{M}_{1}u+\mathcal{K}_{0}*\mathcal{M}_{2}u=\psi_{3}.$
(1.5)

The functions $\psi_{i}=\psi_{i}(x,t)$ are prescribed, as well as the summable
memory kernel $\mathcal{K}_{0}$, while $\mathcal{M}_{i}$ are first-order
differential operators, whose precise form, again, will be given in Section 3.
It is then apparent that the aforementioned advanced models of oxygen
transport through capillaries are just particular cases of our problem.

For last few decades, initial and initial-boundary value problems governed by
subdiffusion with and without memory terms (i.e., (1.1) with $\varrho_{1}=1$
and $\varrho_{2}=0$) have been extensively studied via various approaches of
contemporary analysis, such as the qualitative theory of differential
equations and numerical calculus. With no claim of completeness, we recall a
number of published results. Existence, uniqueness, regularity, longtime
behavior of mild, weak and strong solutions to linear and nonlinear initial-
boundary value problems subject to Dirichlet or Neumann boundary conditions
for evolution equations with single-term fractional derivatives in time were
discussed in [1, 10, 15, 34, 49] and references therein. The $L_{p}$-theory
for linear and semilinear subdiffusion equations was analyzed in [16, 51, 52],
whereas for the solvability of the corresponding problems in smooth classes,
we refer to [20, 25, 26, 21, 22, 14, 32]. Concerning the mathematical
treatment of fractional dynamic boundary conditions (with $\mathcal{K}_{0}=0$,
$\varrho_{1}=1$, $\varrho_{2}=0$), global and local solvability, regularity of
solutions to linear and semilinear elliptic and parabolic operators were
discussed in [9, 5, 19, 17, 25, 47]. The physical interpretation of boundary
conditions of this kind can be found in [9, 47].

Evolution equations containing the general integro-differential operator

$\frac{\partial}{\partial t}(\mathcal{N}*u)(\cdot,t),$ (1.6)

where $\mathcal{N}(t)$ is a nonnegative kernel, are studied in the papers [40,
41]. A particular case of this operator is the multi-term fractional
derivative in time

$\sum_{i=1}^{M}q_{i}\mathbf{D}_{t}^{\nu_{i}}u,$

with $0<\nu_{M}<\ldots<\nu_{1}<1$ and $q_{i}\geq 0$. The Cauchy problem for a
general diffusion equation on unbounded domains was discussed in detail in
[18]. Existence and uniqueness along with a maximal principle for initial-
boundary value problems were studied in [32, 21, 35, 36, 31]. Optimal decay
estimates for equations on bounded domains and subject to the homogenous
Dirichlet boundary condition were examined in [48], which shows in particular
that the decay pattern (e.g., exponential, algebraic or logarithmic) depends
on the (positive) kernel $\mathcal{N}$. An initial value problem for a
semilinear differential equation with a fractional operator of the form (1.6)
was examined in [43], where local/global existence and uniqueness of solutions
were established by exploiting the Schauder fixed point theorem. Finally, we
quote [13, 37, 46, 53], where certain explicit and numerical solutions were
constructed to the corresponding initial-boundary value problems to evolution
equations with multi-term fractional derivatives with $q_{i}>0$.

Coming to equation (1.1) and related problems, we point out two main
differences with respect to the previous literature. The first is related to
the presence of Caputo fractional derivatives of the product of two functions,
that is, $\varrho_{1}u$ and $\varrho_{2}u$. Incidentally, we recall that the
well-known Leibniz rule does not work in the case of fractional Caputo
derivatives. The second difference is that the fractional derivatives
appearing in (1.1), under certain assumptions on $\varrho_{1}$ and
$\varrho_{2}$, can be represented in the form (1.6), but with a _negative_
kernel. Indeed, [11, Lemma 4] tells us that, if $0<\nu_{2}<\nu_{1}<1$,

$\frac{t^{-\nu_{1}}}{\Gamma(1-\nu_{1})}-\frac{t^{-\nu_{2}}}{\Gamma(1-\nu_{2})}<0\qquad\text{whenever}\qquad
t>e^{-\gamma},$

$\gamma$ being the Euler-Mascheroni constant, which in turn provides the
relation

$\mathbf{D}_{t}^{\nu_{1}}(\varrho_{1}u)-\mathbf{D}_{t}^{\nu_{2}}(\varrho_{2}u)=\frac{\partial}{\partial
t}(\mathcal{N}*u)$

for the kernel

$\mathcal{N}(t)=\varrho_{1}\frac{t^{-\nu_{1}}}{\Gamma(1-\nu_{1})}-\varrho_{2}\frac{t^{-\nu_{2}}}{\Gamma(1-\nu_{2})},$

which is negative for $t>e^{-\gamma}$ and
$0<\varrho_{1}(x)\leq\varrho_{2}(x)$, with $\varrho_{i}$ time-independent. In
fact, the nonnegativity of the kernel $\mathcal{N}$ is a key assumption in the
previous works which is removed in our investigation.

The main goal of the present paper is the proof of the well-posedness and the
regularity of a global classical solution to problems (1.1)-(1.5) in smooth
classes for any fixed time $T$, without the assumption on the sign of the
function $\varrho_{2}=\varrho_{2}(x,t)$. This will be obtained by adapting the
technique of a regularizer for parabolic equations [28] to the subdiffusion
equation, so to establish the one-valued global classical solvability of
(1.1)-(1.5). Our analysis is complemented by numerical simulations. It is also
worth observing that, once the linear case is fully understood, it is then
possible to tackle the global classical solvability of boundary-value problems
for nonlinear extensions of (1.1). This will be possibly the subject of future
investigations.

### Outline of the paper

In the next Section 2 we introduce the functional spaces and notations. The
general assumptions are presented in Section 3. The main Theorem 4.1 is stated
in Section 4. Section 5 is devoted to some auxiliary results concerning the
properties of solutions to subdiffusion equations, which will play a key role
in the investigation. In Section 6 we provide the proof of Theorem 4.1,
combining some ideas from [28] with a priory estimates of the solutions. In
the final Section 7 we study the equation from the numerical side.

##  2\. Functional Spaces and Notation

Throughout this work, the symbol $C$ will denote a generic positive constant,
depending only on the structural quantities of the problem. We will carry out
our analysis in the framework of the fractional Hölder spaces. To this end, in
what follows we take two arbitrary (but fixed) parameters

$\alpha\in(0,1)\qquad\text{and}\qquad\nu\in(0,1].$

For any non-negative integer $l$, any Banach space
$(\mathbf{X},\|\cdot\|_{\mathbf{X}}),$ and any $p\geq 1$ and $s\geq 0$, we
consider the usual spaces

$\mathcal{C}^{s}([0,T],\mathbf{X}),\qquad\mathcal{C}^{l+\alpha}(\bar{\Omega}),\qquad
W^{l,p}(\Omega),\qquad L_{p}(\Omega).$

Denoting for $\beta\in(0,1)$

$\displaystyle\langle v\rangle^{(\beta)}_{x,\Omega_{T}}$
$\displaystyle=\sup\Big{\\{}\frac{|v(x_{1},t)-v(x_{2},t)|}{|x_{1}-x_{2}|^{\beta}}:\quad
x_{2}\neq x_{1},\quad x_{1},x_{2}\in\bar{\Omega},\quad t\in[0,T]\Big{\\}},$
$\displaystyle\langle v\rangle^{(\beta)}_{t,\Omega_{T}}$
$\displaystyle=\sup\Big{\\{}\frac{|v(x,t_{1})-v(x,t_{2})|}{|t_{1}-t_{2}|^{\beta}}:\quad
t_{2}\neq t_{1},\quad x\in\bar{\Omega},\quad t_{1},t_{2}\in[0,T]\Big{\\}},$

we have the following definitions.

######  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$ together with its corresponding derivatives are
continuous and the norms here below are finite:

$\displaystyle\|v\|_{\mathcal{C}^{l+\alpha,\frac{l+\alpha}{2}\nu}(\bar{\Omega}_{T})}$
$\displaystyle=\|v\|_{\mathcal{C}([0,T],\mathcal{C}^{l+\alpha}(\bar{\Omega}))}+\sum_{|j|=0}^{l}\langle
D_{x}^{j}v\rangle^{(\frac{l+\alpha-|j|}{2}\nu)}_{t,{\Omega}_{T}},\quad l=0,1,$
$\displaystyle\|v\|_{\mathcal{C}^{2+\alpha,\frac{2+\alpha}{2}\nu}(\bar{\Omega}_{T})}$
$\displaystyle=\|v\|_{\mathcal{C}([0,T],\mathcal{C}^{2+\alpha}(\bar{\Omega}))}+\|\mathbf{D}_{t}^{\nu}v\|_{\mathcal{C}^{\alpha,\frac{\alpha}{2}\nu}(\bar{\Omega}_{T})}+\sum_{|j|=1}^{2}\langle
D_{x}^{j}v\rangle^{(\frac{2+\alpha-|j|}{2}\nu)}_{t,{\Omega}_{T}}.$

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 [22, Section 2]. It is worth
noting that, in the limiting case $\nu=1$, the class
$\mathcal{C}^{l+\alpha,\frac{l+\alpha}{2}\nu}$ coincides with the usual
parabolic Hölder space $H^{l+\alpha,\frac{l+\alpha}{2}}$ (see e.g.[28,
(1.10)-(1.12)]).

######  Definition 2.2.

For $l=0,1,2,$ we define
$\mathcal{C}^{l+\alpha,\frac{l+\alpha}{2}\nu}_{0}(\bar{\Omega}_{T})$ to be the
space consisting of those functions
$v\in\mathcal{C}^{l+\alpha,\frac{l+\alpha}{2}\nu}(\bar{\Omega}_{T})$
satisfying the zero initial conditions:

$v|_{t=0}=0\qquad\text{and}\qquad\underbrace{\mathbf{D}_{t}^{\nu}...\mathbf{D}_{t}^{\nu}}_{m-times}v|_{t=0}=0,\quad
m=0,\ldots,\Big{\lfloor}\frac{l}{2}\Big{\rfloor},$

where $\lfloor\cdot\rfloor$ denotes the floor function.

In a similar manner we introduce the space
$\mathcal{C}^{l+\alpha,\frac{l+\alpha}{2}\nu}_{0}(\partial\Omega_{T})$.

##  3\. General Assumptions

We begin to state our general hypothesis on the structural terms appearing in
the equation and in the boundary conditions.

H1. Conditions on the fractional order of the derivatives:

We assume that

$\nu_{1}\in(0,1]\qquad\text{and}\qquad\nu_{2}\in\begin{cases}\Big{(}0,\frac{\nu_{1}(2-\alpha)}{2}\Big{)},\quad\text{if
either {DBC} or {3BC} hold},\\\ \\\
\Big{(}0,\frac{\nu_{1}(1-\alpha)}{2}\Big{)},\quad\text{if {FDBC}
holds}.\end{cases}$

H2. Conditions on the operators:

The operators appearing in (1.1), (1.4) and (1.5) read

$\mathcal{L}_{1}=\sum_{ij=1}^{n}a_{ij}(x,t)\frac{\partial^{2}}{\partial
x_{i}\partial x_{j}}+\sum_{i=1}^{n}a_{i}(x,t)\frac{\partial}{\partial
x_{i}}+a_{0}(x,t),$
$\mathcal{L}_{2}=\sum_{ij=1}^{n}b_{ij}(x,t)\frac{\partial^{2}}{\partial
x_{i}\partial x_{j}}+\sum_{i=1}^{n}b_{i}(x,t)\frac{\partial u}{\partial
x_{i}}+b_{0}(x,t),$

and

$\mathcal{M}_{1}=\sum_{i=1}^{n}c_{i}(x,t)\frac{\partial}{\partial
x_{i}}+c_{0}(x,t),$
$\mathcal{M}_{2}=\sum_{i=1}^{n}d_{i}(x,t)\frac{\partial}{\partial
x_{i}}+d_{0}(x,t).$

There exist constants $\mu_{2}>\mu_{1}>0,$ $\mu_{3}>0$ and $\mu_{0}>0$ such
that

$\mu_{1}|\xi|^{2}\leq\sum_{ij=1}^{n}a_{ij}(x,t)\xi_{i}\xi_{j}\leq\mu_{2}|\xi|^{2},$

for any $(x,t,\xi)\in\bar{\Omega}_{T}\times\mathbb{R}^{n}$;

$\varrho_{1}(x,t)\geq\mu_{0}>0,$

for any $(x,t)\in\bar{\Omega}_{T}$;

$\sum_{i=1}^{n}c_{i}(x,t)N_{i}(x)\geq\mu_{3}>0,$

for any $(x,t)\in\partial\Omega_{T}$, where $N=\\{N_{1}(x),...,N_{n}(x)\\}$ is
the unit outward normal vector to $\Omega$.

H3. Conditions on the coefficients:

For $i,j=1,\ldots,n$,

$a_{ij}(x,t),a_{i}(x,t),a_{0}(x,t),b_{ij}(x,t),b_{i}(x,t),b_{0}(x,t)\in\mathcal{C}^{\alpha,\frac{\alpha\nu_{1}}{2}}(\bar{\Omega}_{T}),$

and

$c_{i}(x,t),c_{0}(x,t),d_{i}(x,t),d_{0}(x,t)\in\mathcal{C}^{1+\alpha,\frac{1+\alpha}{2}\nu_{1}}(\partial\Omega_{T}).$

We assume that

$\varrho_{1}\in\begin{cases}\mathcal{C}^{\gamma_{0}}([0,T],\mathcal{C}^{1}(\bar{\Omega}))\qquad\qquad\qquad\qquad\qquad\text{in
the {DBC} and {3BC} cases,}\\\
\mathcal{C}^{\gamma_{0}}([0,T],\mathcal{C}^{1}(\bar{\Omega}))\cap\mathcal{C}^{\gamma_{3}}([0,T],\mathcal{C}^{2}(\partial\Omega))\quad\text{in
the {FDBC} case,}\end{cases}$
$\varrho_{2}\in\begin{cases}\mathcal{C}^{\gamma_{1}}([0,T],\mathcal{C}^{1}(\bar{\Omega}))\qquad\qquad\qquad\qquad\qquad\text{in
the {DBC} and {3BC} cases,}\\\
\mathcal{C}^{\gamma_{1}}([0,T],\mathcal{C}^{1}(\bar{\Omega}))\cap\mathcal{C}^{\gamma_{4}}([0,T],\mathcal{C}^{2}(\partial\Omega))\quad\text{in
the {FDBC} case,}\end{cases}$

where

$\gamma_{0}>\frac{\nu_{1}(2+\alpha)}{2},\qquad\gamma_{1}>\frac{\nu_{2}(2+\alpha)}{2},\qquad\gamma_{3}>\frac{\nu_{1}(3+\alpha)}{2},\qquad\gamma_{4}>\frac{\nu_{2}(3+\alpha)}{2}.$

Besides, if $\gamma_{0}$ or/and $\gamma_{1}<1$; and $\gamma_{3}$ or/and
$\gamma_{4}<1$, then we additionally require that
$\mathbf{D}_{t}^{\nu_{1}}\varrho_{1}$ or/and
$\mathbf{D}_{t}^{\nu_{2}}\varrho_{2}\in\mathcal{C}^{\alpha,\frac{\alpha}{2}\nu_{1}}(\bar{\Omega}_{T})$;
and $\mathbf{D}_{t}^{\nu_{1}}\varrho_{1}$ or/and
$\mathbf{D}_{t}^{\nu_{2}}\varrho_{2}\in\mathcal{C}^{1+\alpha,\frac{1+\alpha}{2}\nu_{1}}(\partial\Omega_{T})\cap\mathcal{C}^{\alpha,\frac{\alpha}{2}\nu_{1}}(\bar{\Omega}_{T})$,
respectively.

H4. Conditions on the given functions:

$\mathcal{K}(t),\mathcal{K}_{0}(t)\in L_{1}(0,T),$ $u_{0}(x)\in
C^{2+\alpha}(\bar{\Omega}),\quad
f(x,t)\in\mathcal{C}^{\alpha,\frac{\alpha\nu_{1}}{2}}(\bar{\Omega}_{T}),$
$\psi_{1}(x,t)\in\mathcal{C}^{2+\alpha,\frac{2+\alpha}{2}\nu_{1}}(\partial\Omega_{T}),\quad\psi_{2}(x,t),\psi_{3}(x,t)\in\mathcal{C}^{1+\alpha,\frac{1+\alpha}{2}\nu_{1}}(\partial\Omega_{T}).$

H5. Compatibility conditions:

The following compatibility conditions hold for every $x\in\partial\Omega$ at
the initial time $t=0$:

$\psi_{1}(x,0)=u_{0}(x),\quad\mathbf{D}_{t}^{\nu_{1}}(\varrho_{1}\psi_{1})|_{t=0}-\mathbf{D}_{t}^{\nu_{2}}(\varrho_{2}\psi_{1})|_{t=0}=\mathcal{L}_{1}u_{0}(x)|_{t=0}+f(x,0),$

if the DBC holds; and

$\mathcal{M}_{1}u_{0}(x)|_{t=0}=\psi_{2}(x,0),$

if the 3BC takes place; and in the case of FDBC

$\mathcal{L}_{1}u_{0}(x)|_{t=0}+f(x,0)=\psi_{2}(x,0)+\mathcal{M}_{1}u_{0}(x)|_{t=0}.$

Assumption H2 on the coefficients $c_{i}$ means that the vector
$c=\\{c_{1}(x,t),...,$ $c_{n}(x,t)\\}$ does not lie in the tangent plane to
$\partial\Omega$ at any point.

######  Remark 3.1.

Thanks to Lemma 4.1 in [21], the equalities

$(\mathcal{K}*\mathcal{L}_{2}u)(x,0)=0\qquad\text{and}\qquad(\mathcal{K}_{0}*\mathcal{M}_{2}u)(x,0)=0$

hold for any
$u\in\mathcal{C}^{2+\alpha,\frac{2+\alpha}{2}\nu_{1}}(\partial\Omega_{T})$ and
any $x\in\partial\Omega$. This explains the absence of the memory terms
$(\mathcal{K}*\mathcal{L}_{2}u)$ and $(\mathcal{K}_{0}*\mathcal{M}_{2}u)$ in
the compatibility condition H5.

######  Remark 3.2.

In the case of FDBC, the compatibility H5 and assumptions H3 and H4 provide
the regularity

$(\mathcal{L}_{1}u_{0}|_{t=0}+f(x,0))\in\mathcal{C}^{1+\alpha}(\partial\Omega).$

##  4\. Main Results

We are now ready to state our main result related to the global classical
solvability of (1.1)-(1.5).

######  Theorem 4.1.

Let $T>0$ be fixed, $\partial\Omega\in\mathcal{C}^{2+\alpha}$, and let
assumptions H1-H5 hold. Then equation (1.1) with the initial condition (1.2),
subject to either boundary condition DBC, 3BC or FDBC admits a unique
classical solution $u=u(x,t)$ on $\bar{\Omega}_{T}$, satisfying the regularity
$u\in\mathcal{C}^{2+\alpha,\frac{2+\alpha}{2}\nu_{1}}(\bar{\Omega}_{T})$.
Besides, this solution fulfills the estimate

$\displaystyle\|u\|_{\mathcal{C}^{2+\alpha,\frac{2+\alpha}{2}\nu_{1}}(\bar{\Omega}_{T})}+\|\mathbf{D}_{t}^{\nu_{2}}u\|_{\mathcal{C}^{\alpha,\frac{\alpha}{2}\nu_{1}}(\bar{\Omega}_{T})}$
$\displaystyle\leq
C\begin{cases}\|f\|_{\mathcal{C}^{\alpha,\frac{\alpha}{2}\nu_{1}}(\bar{\Omega}_{T})}+\|u_{0}\|_{\mathcal{C}^{2+\alpha}(\bar{\Omega})}+\|\psi_{1}\|_{\mathcal{C}^{2+\alpha,\frac{2+\alpha}{2}\nu_{1}}(\partial\Omega_{T})}\quad\text{in
the {\rm{DBC}} case},\\\
\|f\|_{\mathcal{C}^{\alpha,\frac{\alpha}{2}\nu_{1}}(\bar{\Omega}_{T})}+\|u_{0}\|_{\mathcal{C}^{2+\alpha}(\bar{\Omega})}+\|\psi_{2}\|_{\mathcal{C}^{1+\alpha,\frac{1+\alpha}{2}\nu_{1}}(\partial\Omega_{T})}\quad\text{in
the {\rm{3BC}} case},\end{cases}$

while if the FDBC case holds then

$\displaystyle\|u\|_{\mathcal{C}^{2+\alpha,\frac{2+\alpha}{2}\nu_{1}}(\bar{\Omega}_{T})}+\|\mathbf{D}_{t}^{\nu_{2}}u\|_{\mathcal{C}^{\alpha,\frac{\alpha}{2}\nu_{1}}(\bar{\Omega}_{T})\cap\mathcal{C}^{1+\alpha,\frac{1+\alpha}{2}\nu_{1}}(\partial\Omega_{T})}+\|\mathbf{D}_{t}^{\nu_{1}}u\|_{\mathcal{C}^{1+\alpha,\frac{1+\alpha}{2}\nu_{1}}(\partial\Omega_{T})}$
$\displaystyle\leq
C[\|f\|_{\mathcal{C}^{\alpha,\frac{\alpha}{2}\nu_{1}}(\bar{\Omega}_{T})}+\|u_{0}\|_{\mathcal{C}^{2+\alpha}(\bar{\Omega})}+\|\psi_{3}\|_{\mathcal{C}^{1+\alpha,\frac{1+\alpha}{2}\nu_{1}}(\partial\Omega_{T})}].$

The generic constant $C$ is independent of the right-hand sides of
(1.1)-(1.5).

Indeed the positive constant $C$ depends only on the Lebesgue measures of
$\Omega$ and its boundary $\partial\Omega$, on the norm
$\|\mathcal{K}\|_{L_{1}(0,T)}$, and on the norms of the coefficients of the
operators $\mathcal{L}_{i}$ (as well as $\mathcal{M}_{1}$, $\mathcal{M}_{2}$
and $\|\mathcal{K}_{0}\|_{L_{1}(0,T)}$ in the case of 3BC and FDBC), and the
corresponding norms of $\varrho_{1}$ and $\varrho_{2}$.

######  Remark 4.2.

It is worth noting that our assumptions on the kernels $\mathcal{K}$ and
$\mathcal{K}_{0}$ include the case $\mathcal{K}=\mathcal{K}_{0}\equiv 0$,
meaning that the multi-term subdiffusion equation:

$\mathbf{D}_{t}^{\nu_{1}}(\varrho_{1}u)-\mathbf{D}_{t}^{\nu_{2}}(\varrho_{2}u)-\mathcal{L}_{1}u=f$

fits in our analysis and is described by the theorem above.

######  Remark 4.3.

Actually, assumptions H2, H3 on the coefficients $c_{i}$, $c_{0}$ and
condition H4 on the right-hand side $\psi_{2}$ tell us that initial-value
problem (1.1)-(1.2) subject to the Neumann boundary condition (NBC), that is,

$\frac{\partial u}{\partial N}+\mathcal{K}_{0}*\frac{\partial u}{\partial
N}=\psi_{2}\quad\text{on}\quad\partial\Omega_{T},$

is just a particular case of problem (1.1), (1.2), (1.4) with
$c_{0}=d_{0}\equiv 0,$ $c_{i}=d_{i}=N_{i}(x),$ $i=1,2,...,n$. Thus, results of
Theorem 4.1 extend to the NBC.

######  Remark 4.4.

With inessential modification in the proofs, the very same results hold for
the $M$-term fractional equations:

$\mathbf{D}_{t}^{\nu_{1}}(\varrho_{1}u)-\sum_{i=2}^{M}\mathbf{D}_{t}^{\nu_{i}}(\varrho_{i}u)-\mathcal{L}_{1}u-\int_{0}^{t}\mathcal{K}(t-s)\mathcal{L}_{2}u(\cdot,s)ds=f(x,t),$
$\varrho_{1}\mathbf{D}_{t}^{\nu_{1}}u-\sum_{i=2}^{M}\varrho_{i}\mathbf{D}_{t}^{\nu_{i}}u-\mathcal{L}_{1}u-\int_{0}^{t}\mathcal{K}(t-s)\mathcal{L}_{2}u(\cdot,s)ds=f(x,t).$

In these cases, we additionally assume that all $\nu_{i}$ and $\varrho_{i}$
$i=3,...M,$ have the properties of $\nu_{2}$ and $\varrho_{2}$ (see
assumptions H1,H3), besides the second equality in compatibility conditions in
the DBC case takes the form

$\displaystyle\mathbf{D}_{t}^{\nu_{1}}(\varrho_{1}\psi_{1})|_{t=0}$
$\displaystyle-\sum_{i=2}^{M}\mathbf{D}_{t}^{\nu_{i}}(\varrho_{i}\psi_{1})|_{t=0}=\mathcal{L}_{1}u_{0}(x)|_{t=0}+f(x,0)\quad\text{or}$
$\displaystyle\varrho_{1}(x,0)\mathbf{D}_{t}^{\nu_{1}}(\psi_{1})|_{t=0}$
$\displaystyle-\sum_{i=2}^{M}\varrho_{i}(x,0)\mathbf{D}_{t}^{\nu_{i}}(\psi_{1})|_{t=0}=\mathcal{L}_{1}u_{0}(x)|_{t=0}+f(x,0),$

respectively.

Finally, we remark, that in the case of the second equation here, the
regularity of the functions $\varrho_{i},$ $i=1,...,M$, can be relaxed.
Namely, we need in
$\varrho_{i}\in\mathcal{C}^{\alpha,\alpha\nu_{1}/2}(\bar{\Omega}_{T})$ in the
case of DBC or 3BC cases, while
$\varrho_{i}\in\mathcal{C}^{\alpha,\alpha\nu_{1}/2}(\bar{\Omega}_{T})\cap\mathcal{C}^{1+\alpha,(1+\alpha)\nu_{1}/2}(\partial\Omega_{T})$
in FDBC case.

##  5\. Technical Results

We recall some properties of fractional derivatives and integrals, along with
several technical results that will be used in this article. In what follows,
for any $\theta>0$ we denote

$\omega_{\theta}(t)=\frac{t^{\theta-1}}{\Gamma(\theta)}.$

We define the fractional Riemann-Liouville integral and the derivative of
order $\theta$ of a function $v=v(t)$ (possibly also depending on other
variables) as

$I_{t}^{\theta}v(t)=(\omega_{\theta}*v)(t),\qquad\partial_{t}^{\theta}v(t)=\frac{\partial^{\lceil\theta\rceil}}{\partial
t^{\lceil\theta\rceil}}(\omega_{\lceil\theta\rceil-\theta}*v)(t),$

respectively, where $\lceil\theta\rceil$ is the ceiling function of $\theta$
(i.e. the smallest integer greater than or equal to $\theta$). In particular,
for $\theta\in(0,1)$

$\partial_{t}^{\theta}v(t)=\frac{\partial}{\partial
t}(\omega_{1-\theta}*v)(t).$

Accordingly, the Caputo fractional derivative of the order $\theta\in(0,1)$
reads

$\mathbf{D}_{t}^{\theta}v(t)=\frac{\partial}{\partial
t}(\omega_{1-\theta}*v)(t)-\omega_{1-\theta}(t)v(0)=\partial_{t}^{\theta}v(t)-\omega_{1-\theta}(t)v(0),$

provided that both derivatives exist.

Our first assertion playing a key point in the proof of Theorem 4.1 describes
the regularity of lower fractional derivatives in time
$\mathbf{D}_{t}^{\beta}w$, with $0<\beta<\nu\leq 1,$ in the case when
$w\in\mathcal{C}_{0}^{2+\alpha,\frac{2+\alpha}{2}\nu}(\bar{\Omega}_{T})$. To
this end, we define

$\Omega^{r}=\Omega\cap
B_{r}(x^{0}),\quad\Omega_{T}^{r}=\Omega^{r}\times(0,T),\quad\partial\Omega^{\star}=\partial\Omega\cap\partial\Omega^{r},$

where $B_{r}(x^{0})$ is the ball centered at a point $x^{0}\in\bar{\Omega}$ of
radius $r$. Then we introduce the functions

$\xi=\xi(x)\in\mathcal{C}_{0}^{\infty}(\mathbb{R}^{n}),\qquad\xi\in[0,1],\qquad\xi=\begin{cases}1,\quad
x\in\Omega^{r},\\\ 0,\quad
x\in\mathbb{R}^{n}\backslash\bar{\Omega}^{2r},\end{cases}$

and

$\mathfrak{J}_{\theta}(t)=\mathfrak{J}_{\theta}(t;w_{1},w_{2})=\int_{0}^{t}\frac{[w_{2}(\cdot,t)-w_{2}(\cdot,s)][w_{1}(\cdot,s)-w_{1}(\cdot,0)]}{(t-s)^{\theta+1}}ds$

where $\theta\in(0,1)$, $w_{1}$ and $w_{2}$ are some given smooth functions.

######  Lemma 5.1.

Let $x^{0}\in\bar{\Omega}$ be arbitrarily fixed, let $0<\beta<\nu\leq 1.$ We
assume that

$w_{1}\in\mathcal{C}_{0}^{2+\alpha,\frac{2+\alpha}{2}\nu}(\bar{\Omega}_{T})\qquad\text{and}\qquad
w_{2}\in\mathcal{C}^{\gamma}([0,T],\mathcal{C}^{1}(\bar{\Omega})),$

with $\gamma\geq\frac{2+\alpha}{2}\nu$, and we set

$\delta=\min\\{1,\gamma\\}.$

If $\gamma<1$ we additionally require
$\mathbf{D}_{t}^{\nu}w_{2}\in\mathcal{C}^{\alpha,\alpha\nu/2}(\bar{\Omega}_{T})$.
Then, for any $\beta\in(0,\frac{2-\alpha}{2}\nu)$ and any $\tau\in(0,T]$, the
following estimates hold:

i:

$\displaystyle\|w_{2}\xi\mathbf{D}_{t}^{\beta}w_{1}\|_{\mathcal{C}^{\alpha,\alpha\nu/2}(\bar{\Omega}^{r}_{\tau})}$
$\displaystyle\leq
C[\tau^{\nu-\beta+\alpha\nu/2}+\tau^{\nu-\beta+\nu\alpha/2}r^{-\alpha}+\tau^{\nu-\beta}+\tau^{\nu-\beta}r^{-\alpha}]\|\mathbf{D}_{t}^{\nu}w_{1}\|_{\mathcal{C}^{\alpha,\alpha\nu/2}(\bar{\Omega}^{r}_{\tau})}.$

ii:

If $w_{2}(x,0)=0$ then

$\|w_{2}\xi\mathbf{D}_{t}^{\nu}w_{1}\|_{\mathcal{C}^{\alpha,\alpha\nu/2}(\bar{\Omega}^{r}_{\tau})}\leq
C[\tau^{\delta-\alpha\nu/2}+\tau^{\delta}r^{-\alpha}+\tau^{\delta}]\|\mathbf{D}_{t}^{\nu}w_{1}\|_{\mathcal{C}^{\alpha,\alpha\nu/2}(\bar{\Omega}^{r}_{\tau})}.$

iii:

$\|\xi\mathfrak{J}_{\beta}(t)\|_{\mathcal{C}^{\alpha,\alpha\nu/2}(\bar{\Omega}^{r}_{\tau})}\leq
C\tau^{\delta-\beta+\nu(1-\alpha)/2}[\tau^{\nu/2}+\tau^{\nu\alpha}+\tau^{\nu\frac{(1+\alpha)}{2}}(r^{-\alpha}+1)]\|w_{1}\|_{\mathcal{C}^{2+\alpha,\frac{2+\alpha}{2}\nu}(\bar{\Omega}^{r}_{\tau})},$

and

$\|\xi\mathfrak{J}_{\nu}(t)\|_{\mathcal{C}^{\alpha,\alpha\nu/2}(\bar{\Omega}^{r}_{\tau})}\leq
C\tau^{\delta-\nu(1+\alpha)/2}[\tau^{\nu/2}+\tau^{\nu\alpha}+\tau^{\nu\frac{(1+\alpha)}{2}}(r^{-\alpha}+1)]\|w_{1}\|_{\mathcal{C}^{2+\alpha,\frac{2+\alpha}{2}\nu}(\bar{\Omega}^{r}_{\tau})}.$

The positive quantity $C$ depends only on $\nu,\beta,T$, the Lebesgue measure
of $\Omega$ and the norm of $w_{2}$.

###### Proof.

We start with the evaluation of the term $\mathbf{D}_{t}^{\beta}w_{1}$. To
this end, appealing to representation (10.34) in [24], we deduce that

$w_{2}\xi\mathbf{D}_{t}^{\beta}w_{1}=i_{1}+i_{2},$

where we put

$\displaystyle i_{1}$
$\displaystyle=i_{1}(x,t)=\xi(x)w_{2}(x,t)\frac{{t}^{\nu-\beta}}{\Gamma(1+\nu-\beta)}\mathbf{D}_{t}^{\nu}w_{1}(x,t),$
$\displaystyle i_{2}$
$\displaystyle=i_{2}(x,t)=\xi(x)w_{2}(x,t)\int_{0}^{t}\frac{(t-s)^{\nu-\beta-1}}{\Gamma(\nu-\beta)}[\mathbf{D}_{t}^{\nu}w_{1}(x,t)-\mathbf{D}_{s}^{\nu}w_{1}(x,s)]ds.$

We estimate the norms of $i_{1}$ and $i_{2}$ separately.

As for
$\|i_{1}\|_{\mathcal{C}^{\alpha,\alpha\nu/2}(\bar{\Omega}^{r}_{\tau})},$ the
desired bound is a simple consequence of the following easily verified
relations:

$\displaystyle\underset{\bar{\Omega}^{r}_{\tau}}{\sup}\,|i_{1}|$
$\displaystyle\leq
C\tau^{\nu-\beta}\underset{\bar{\Omega}^{r}_{\tau}}{\sup}\,|w_{2}|\,\underset{\bar{\Omega}^{r}_{\tau}}{\sup}\,|\mathbf{D}_{t}^{\nu}w_{1}|,$
$\displaystyle\langle i_{1}\rangle_{x,\Omega^{r}_{\tau}}^{(\alpha)}$
$\displaystyle\leq\tau^{\nu-\beta}\Big{[}\langle
w_{2}\rangle_{x,\bar{\Omega}^{r}_{\tau}}^{(\alpha)}\underset{\bar{\Omega}^{r}_{\tau}}{\sup}\,|\mathbf{D}_{t}^{\nu}w_{1}|+\underset{\bar{\Omega}^{r}_{\tau}}{\sup}\,|w_{2}|\langle\xi\mathbf{D}_{t}^{\nu}w_{1}\rangle_{x,\bar{\Omega}^{r}_{\tau}}^{(\alpha)}\Big{]}$
$\displaystyle\leq
C[r^{1-\alpha}\tau^{\nu-\beta}+\tau^{\nu-\beta}r^{-\alpha}+\tau^{\nu-\beta}]\|\mathbf{D}_{t}^{\nu}w_{1}\|_{\mathcal{C}([0,\tau],\mathcal{C}^{\alpha}(\bar{\Omega}))},$

and

$\displaystyle\langle i_{1}\rangle_{t,\Omega^{r}_{\tau}}^{(\alpha\nu/2)}$
$\displaystyle\leq C\Big{[}\tau^{\nu-\beta+\delta-\alpha\nu/2}\langle
w_{2}\rangle_{t,\Omega^{r}_{\tau}}^{(\delta)}\|\mathbf{D}_{t}^{\nu}w_{1}\|_{\mathcal{C}(\bar{\Omega}_{\tau}^{r})}+\langle
t^{\nu-\beta}\mathbf{D}_{t}^{\nu}w_{1}\rangle_{t,\Omega^{r}_{\tau}}^{(\alpha\nu/2)}\|w_{2}\|_{\mathcal{C}(\bar{\Omega}_{T})}\Big{]}$
$\displaystyle\leq
C\tau^{\nu-\beta-\alpha\nu/2}[1+\tau^{\alpha\nu/2}+\tau^{\delta}]\|w_{2}\|_{\mathcal{C}^{\gamma}([0,T],\mathcal{C}^{1}(\bar{\Omega}))}\|\mathbf{D}_{t}^{\nu}w_{1}\|_{\mathcal{C}^{\alpha\nu/2}([0,\tau],\mathcal{C}(\bar{\Omega}^{r}))}.$

Coming to
$\|i_{2}\|_{\mathcal{C}^{\alpha,\alpha\nu/2}(\bar{\Omega}_{\tau}^{r})}$, we
easily find

$\|i_{2}\|_{\mathcal{C}(\bar{\Omega}_{\tau}^{r})}+\langle
i_{2}\rangle_{x,\Omega_{\tau}^{r}}^{(\alpha)}\leq
C\tau^{\nu-\beta}\Big{[}\tau^{\alpha\nu/2}(1+r^{-\alpha})\langle\mathbf{D}_{t}^{\nu}w_{1}\rangle_{t,\Omega_{\tau}^{r}}^{(\alpha\nu/2)}+\langle\mathbf{D}_{t}^{\nu}w_{1}\rangle_{x,\Omega_{\tau}^{r}}^{(\alpha)}\Big{]}\|w_{2}\|_{\mathcal{C}([0,T],\mathcal{C}^{\alpha}(\bar{\Omega}))}.$

In order to complete the estimate of $i_{2}$, hence establishing point (i), we
are left to examine the difference $|i_{2}(x,t_{2})-i_{2}(x,t_{1})|$. To this
end, assuming $t_{2}>t_{1}$ and setting $\Delta t=t_{2}-t_{1}$, we have

$|i_{2}(x,t_{2})-i_{2}(x,t_{1})|\leq\sum_{j=1}^{4}i_{2,j},$

where

$\displaystyle i_{2,1}$
$\displaystyle=\xi|w_{2}(x,t_{2})-w_{2}(x,t_{1})|\Big{|}\int_{0}^{t_{2}}\frac{s^{\nu-\beta-1}}{\Gamma(\nu-\beta)}[\mathbf{D}_{t}^{\nu}w_{1}(x,t_{2}-s)-\mathbf{D}_{t}^{\nu}w_{1}(x,t_{2})]ds\Big{|},$
$\displaystyle i_{2,2}$
$\displaystyle=\xi|w_{2}(x,t_{1})|\Big{|}\int_{0}^{t_{1}}\frac{s^{\nu-\beta-1}}{\Gamma(\nu-\beta)}[\mathbf{D}_{t}^{\nu}w_{1}(x,t_{2}-s)-\mathbf{D}_{t}^{\nu}w_{1}(x,t_{1}-s)]ds\Big{|},$
$\displaystyle i_{2,3}$
$\displaystyle=\xi|w_{2}(x,t_{1})||\mathbf{D}_{t}^{\nu}w_{1}(x,t_{2})-\mathbf{D}_{t}^{\nu}w_{1}(x,t_{1})|\bigg{|}\int_{0}^{t_{1}}\frac{s^{\nu-\beta-1}}{\Gamma(\nu-\beta)}ds\bigg{|},$
$\displaystyle i_{2,4}$
$\displaystyle=\xi|w_{2}(x,t_{1})|\bigg{|}\int_{t_{1}}^{t_{2}}\frac{s^{\nu-\beta-1}}{\Gamma(\nu-\beta)}[\mathbf{D}_{t}^{\nu}w_{1}(x,t_{2}-s)-\mathbf{D}_{t}^{\nu}w_{1}(x,t_{2})]ds\bigg{|}.$

Exploiting the smoothness of the functions $w_{2}$ and
$\mathbf{D}_{t}^{\nu}w_{1}$, and taking into account of the relation between
$\nu$ and $\beta$, we arrive at the sought estimate for $\langle
i_{2}\rangle_{t,\Omega_{\tau}^{r}}^{(\alpha\nu/2)}$.

To verify statement (ii), it is worth noting that the assumptions on $w_{2}$
(i.e. $w_{2}(x,0)=0$) provide the estimate

$\|w_{2}\|_{\mathcal{C}^{\alpha,\alpha\nu/2}(\bar{\Omega}_{\tau}^{r})}\leq
C\tau^{\delta-\alpha\nu/2}[\tau^{\alpha\nu/2}+1+r^{1-\alpha}\tau^{\alpha\nu/2}]\|w_{2}\|_{\mathcal{C}^{\delta}([0,T],\mathcal{C}^{1}(\bar{\Omega}))}.$

Collecting this bound with the regularity of $\mathbf{D}_{t}^{\nu}w_{1}$, the
desired claim follows.

Coming to (iii), we restrict ourselves to the verification of the first
inequality, for the second one is deduced in a similar manner. Straightforward
calculations lead to the relations

$\underset{\bar{\Omega}_{\tau}^{r}}{\sup}\,|\xi\mathfrak{J}_{\beta}(t)|\leq
C\tau^{\delta-\beta+\nu}\langle
w_{2}\rangle_{t,\Omega_{\tau}^{r}}^{(\delta)}\|\mathbf{D}_{t}^{\nu}w_{1}\|_{\mathcal{C}(\bar{\Omega}_{\tau}^{r})},$
$\displaystyle\Big{|}\xi(x_{2})\mathfrak{J}_{\beta}(t)|_{x=x_{1}}-\xi(x_{1})\mathfrak{J}_{\beta}(t)|_{x=x_{2}}\Big{|}$
$\displaystyle\leq
Cr^{-\alpha}\tau^{\delta-\beta+\nu}\|\mathbf{D}_{t}^{\nu}w_{1}\|_{\mathcal{C}(\bar{\Omega}_{\tau}^{r})}\langle
w_{2}\rangle_{t,\Omega_{\tau}^{r}}^{(\delta)}|x_{1}-x_{2}|^{\alpha}+C\Big{|}\mathfrak{J}_{\beta}(t)|_{x=x_{1}}-\mathfrak{J}_{\beta}(t)|_{x=x_{2}}\Big{|},$

and

$\displaystyle\Big{|}\mathfrak{J}_{\beta}(t)|_{x=x_{1}}-\mathfrak{J}_{\beta}(t)|_{x=x_{2}}\Big{|}$
$\displaystyle\leq C\Big{\\{}|x_{2}-x_{1}|\langle
D_{x}w_{1}\rangle_{t,\Omega_{\tau}^{r}}^{(\frac{1+\alpha}{2}\nu)}\int\limits_{0}^{t}\frac{|w_{2}(x_{2},s)-w_{2}(x_{2},t)|s^{\frac{1+\alpha}{2}\nu}}{(t-s)^{1+\beta}}ds$
$\displaystyle\quad+\int\limits_{0}^{t}\frac{|w_{2}(x_{1},t)-w_{2}(x_{2},t)+w_{2}(x_{2},s)-w_{2}(x_{1},s)|s^{\nu}\|\mathbf{D}_{t}^{\nu}w_{1}\|_{\mathcal{C}(\bar{\Omega}_{\tau}^{r})}}{(t-s)^{1+\beta}}ds\Big{\\}}$
$\displaystyle\leq
C|x_{1}-x_{2}|^{\alpha}\tau^{\delta-\beta+\frac{\nu(1+\alpha)}{2}}(1+r^{1-\alpha})\|w_{2}\|_{\mathcal{C}^{\delta}([0,T],\mathcal{C}^{1}(\bar{\Omega}))}\|w_{1}\|_{\mathcal{C}^{2+\alpha,\frac{2+\alpha}{2}\nu}(\bar{\Omega}_{\tau}^{r})},$

which in turn entail

$\displaystyle\|\xi\mathfrak{J}_{\beta}(t)\|_{\mathcal{C}([0,\tau],\mathcal{C}^{\alpha}(\bar{\Omega}^{r}))}$
$\displaystyle\leq
C\tau^{\delta-\beta+\frac{\nu(1+\alpha)}{2}}(r^{1-\alpha}+\tau^{(1-\alpha)\nu/2}r^{-\alpha}+1)\|w_{1}\|_{\mathcal{C}^{2+\alpha,\frac{2+\alpha}{2}\nu}(\bar{\Omega}^{r}_{\tau})}\|w_{2}\|_{\mathcal{C}^{\gamma}([0,T],\mathcal{C}^{1}(\bar{\Omega}))}.$
(5.1)

Thus, taking into account Definition 2.1, we will complete the proof of the
first bound in (iii) if we obtain the corresponding estimate of the seminorm
$\langle\xi\mathfrak{J}_{\beta}(t)\rangle_{t,\Omega_{\tau}^{r}}^{(\alpha\nu/2)}$.
To this end, assuming $T\geq t_{2}>t_{1}\geq 0$, let us define

$\Delta t=t_{2}-t_{1}.$

If $\Delta t\geq t_{1}/2$, it is apparent that

$|\mathfrak{J}_{\beta}(t_{2})-\mathfrak{J}_{\beta}(t_{1})|\leq
C\sum_{j=1}^{3}b_{j},$

where

$\displaystyle b_{1}$
$\displaystyle=\int_{0}^{t_{1}}\frac{|w_{2}(x,t_{1})-w_{2}(x,s)||w_{1}(x,s)-w_{1}(x,0)|}{(t_{1}-s)^{1+\beta}}ds,$
$\displaystyle b_{2}$
$\displaystyle=\int_{0}^{t_{1}}\frac{|w_{2}(x,t_{2})-w_{2}(x,s)||w_{1}(x,s)-w_{1}(x,0)|}{(t_{2}-s)^{1+\beta}}ds,$
$\displaystyle b_{3}$
$\displaystyle=\int^{t_{2}}_{t_{1}}\frac{|w_{2}(x,t_{2})-w_{2}(x,s)||w_{1}(x,s)-w_{1}(x,0)|}{(t_{2}-s)^{1+\beta}}ds.$

Then, appealing to the smoothness of the functions $w_{1}$ and $w_{2}$, we
conclude that

$\displaystyle b_{1}$ $\displaystyle\leq
C\int_{0}^{t_{1}}(t_{1}-s)^{\delta-\beta-1}s^{\nu}ds\langle
w_{2}\rangle_{t,\Omega_{T}}^{(\delta)}\|\mathbf{D}^{\nu}_{t}w_{1}\|_{\mathcal{C}(\bar{\Omega}_{\tau}^{r})}$
$\displaystyle\leq C(\Delta
t)^{\alpha\nu/2}\tau^{\delta+\nu-\beta-\alpha\nu/2}\|\mathbf{D}^{\nu}_{t}w_{1}\|_{\mathcal{C}(\bar{\Omega}_{\tau}^{r})}\langle
w_{2}\rangle_{t,\Omega_{T}}^{(\delta)},$

and

$\displaystyle b_{3}$
$\displaystyle\leq\|\mathbf{D}^{\nu}_{t}w_{1}\|_{\mathcal{C}(\bar{\Omega}_{\tau}^{r})}\langle
w_{2}\rangle_{t,\Omega_{T}}^{(\delta)}t_{2}^{\nu}\int_{t_{1}}^{t_{2}}(t_{2}-s)^{\delta-1-\beta}ds$
$\displaystyle\leq C\tau^{\delta+\nu-\beta-\alpha\nu/2}(\Delta
t)^{\alpha\nu/2}\|\mathbf{D}^{\nu}_{t}w_{1}\|_{\mathcal{C}(\bar{\Omega}_{\tau}^{r})}\langle
w_{2}\rangle_{t,\Omega_{T}}^{(\gamma^{\star})}.$

The estimate for $b_{2}$ is analogous to the one of $b_{1}$. Taking into
account (5), this yields the desired bound in (iii) when $\Delta t\geq
t_{1}/2$. If instead $\Delta t<t_{1}/2$, we rewrite the difference as

$|\mathfrak{J}_{\beta}(t_{2})-\mathfrak{J}_{\beta}(t_{1})|\leq
C\sum_{j=1}^{4}a_{j},$ (5.2)

where

$\displaystyle a_{1}$ $\displaystyle=\int_{t_{1}-2\Delta
t}^{t_{1}}\frac{|w_{2}(x,t_{2})-w_{2}(x,s)||w_{1}(x,s)-w_{1}(x,0)|}{(t_{2}-s)^{1+\beta}}ds,$
$\displaystyle a_{2}$ $\displaystyle=\int_{t_{1}-2\Delta
t}^{t_{1}}\frac{|w_{2}(x,t_{1})-w_{2}(x,s)||w_{1}(x,s)-w_{1}(x,0)|}{(t_{1}-s)^{1+\beta}}ds,$
$\displaystyle a_{3}$ $\displaystyle=\int_{0}^{t_{1}-2\Delta
t}|w_{2}(x,t_{2})-w_{2}(x,s)||w_{1}(x,s)-w_{1}(x,0)|\Big{|}(t_{2}-s)^{-1-\beta}-(t_{1}-s)^{-1-\beta}\Big{|}ds,$
$\displaystyle a_{4}$ $\displaystyle=\int^{t_{1}-2\Delta
t}_{0}\frac{|w_{2}(x,t_{2})-w_{2}(x,t_{1})||w_{1}(x,s)-w_{1}(x,0)|}{(t_{1}-s)^{1+\beta}}ds.$

Due to the properties of the functions $w_{2}$ and $w_{1}$, and exploiting the
mean-value theorem in the evaluation of the term $a_{3}$, we end up with

$|\mathfrak{J}_{\beta}(t_{2})-\mathfrak{J}_{\beta}(t_{1})|\leq C\langle
w_{2}\rangle_{t,\Omega_{T}}^{(\delta)}\|\mathbf{D}_{t}^{\nu}w_{1}\|_{\mathcal{C}(\bar{\Omega}_{\tau}^{r})}(\Delta
t)^{\alpha\nu/2}\tau^{\delta+\nu-\beta-\alpha\nu/2},$

which completes the argument. ∎

Recasting the same proof, we immediately obtain

######  Lemma 5.2.

Let the assumptions of Lemma 5.1 hold. Besides, let
$\gamma>\frac{(3+\alpha)\nu}{2}$ and

$\mathbf{D}_{t}^{\nu}w_{1}\in\mathcal{C}_{0}^{1+\alpha,\frac{1+\alpha}{2}\nu}(\partial\Omega_{T})\qquad\text{and}\qquad
w_{2}\in\mathcal{C}^{\gamma}([0,T],\mathcal{C}^{2}(\partial\Omega)).$

If $\gamma<1$ we also require
$\mathbf{D}_{t}^{\nu}w_{2}\in\mathcal{C}_{0}^{1+\alpha,\frac{1+\alpha}{2}\nu}(\partial\Omega_{T})$.
Then, for $\beta\in(0,\frac{\nu(1-\alpha)}{2})$, with $\delta$ as above, the
following estimates hold:

i:

$\displaystyle\|w_{2}\xi\mathbf{D}_{t}^{\beta}w_{1}\|_{\mathcal{C}^{1+\alpha,\frac{1+\alpha}{2}\nu}(\partial\Omega^{\star}_{\tau})}$
$\displaystyle\leq
C[\tau^{\nu-\beta}+\tau^{\nu\frac{(1-\alpha)}{2}-\beta}+\tau^{\nu-\beta}r^{-1-\alpha}+\tau^{\nu-\beta}r^{-\alpha}]\|\mathbf{D}_{t}^{\nu}w_{1}\|_{\mathcal{C}^{1+\alpha,\frac{1+\alpha}{2}\nu}(\partial\Omega^{\star}_{\tau})}.$

ii:

If $w_{2}(x,0)=0$ then

$\displaystyle\|w_{2}\xi\mathbf{D}_{t}^{\nu}w_{1}\|_{\mathcal{C}^{1+\alpha,\frac{1+\alpha}{2}\nu}(\partial\Omega^{\star}_{\tau})}$
$\displaystyle\leq
C[\tau^{\delta-\frac{(\alpha+1)}{2}\nu}+\tau^{\delta}r^{-\alpha-1}+\tau^{\delta-\frac{\alpha\nu}{2}}r^{-1}][\|\mathbf{D}_{t}^{\nu}w_{1}\|_{\mathcal{C}^{1+\alpha,\frac{1+\alpha}{2}\nu}(\partial\Omega^{\star}_{\tau})}+\|w_{1}\|_{\mathcal{C}^{2+\alpha,\frac{2+\alpha}{2}\nu}(\bar{\Omega}^{r}_{\tau})}].$

iii:

$\displaystyle\|\xi\mathfrak{J}_{\beta}(t)\|_{\mathcal{C}^{1+\alpha,\frac{1+\alpha}{2}\nu}(\partial\Omega^{\star}_{\tau})}$
$\displaystyle\leq
C\frac{\tau^{\delta-\beta}}{r^{1+\alpha}}[\tau^{\nu}+\tau^{\beta}+r^{\alpha}(\tau^{\beta-\nu\alpha/2}+\tau^{\nu\frac{(1+\alpha)}{2}}+\tau^{\nu\frac{2-\alpha}{2}}+r\tau^{\nu\frac{1-\alpha}{2}}]$
$\displaystyle\quad\cdot[\|\mathbf{D}_{t}^{\nu}w_{1}\|_{\mathcal{C}^{1+\alpha,\frac{1+\alpha}{2}\nu}(\partial\Omega^{\star}_{\tau})}+\|w_{1}\|_{\mathcal{C}^{2+\alpha,\frac{2+\alpha}{2}\nu}(\bar{\Omega}^{r}_{\tau})}],$
$\displaystyle\|\xi\mathfrak{J}_{\nu}(t)\|_{\mathcal{C}^{1+\alpha,\frac{1+\alpha}{2}\nu}(\partial\Omega^{\star}_{\tau})}$
$\displaystyle\leq
C\frac{\tau^{\delta-\nu}}{r^{1+\alpha}}[\tau^{\nu}+r^{\alpha}(\tau^{\nu-\nu\alpha/2}+\tau^{\nu\frac{(1+\alpha)}{2}}+\tau^{\nu\frac{2-\alpha}{2}}+r\tau^{\nu\frac{1-\alpha}{2}}]$
$\displaystyle\quad\cdot[\|\mathbf{D}_{t}^{\nu}w_{1}\|_{\mathcal{C}^{1+\alpha,\frac{1+\alpha}{2}\nu}(\partial\Omega^{\star}_{\tau})}+\|w_{1}\|_{\mathcal{C}^{2+\alpha,\frac{2+\alpha}{2}\nu}(\bar{\Omega}^{r}_{\tau})}].$

Again, $C$ depends only on $\nu,\beta,T$, the Lebesgue measure of $\Omega$ and
the norm of $w_{2}$.

######  Remark 5.3.

It is apparent that the estimates of the terms
$\|\xi\mathfrak{J}_{\beta}\|_{\mathcal{C}^{\alpha,\alpha\nu/2}(\bar{\Omega}_{\tau}^{r})}$
and
$\|\xi\mathfrak{J}_{\beta}\|_{\mathcal{C}^{1+\alpha,(1+\alpha)\nu/2}(\partial\Omega^{\star}_{\tau})}$
in points (iii) of both lemmas above hold within weaker assumptions on
$\gamma$, namely, $\gamma>\frac{2+\alpha}{2}\beta$ and
$\gamma>\frac{3+\alpha}{2}\beta$, respectively.

######  Remark 5.4.

The following estimates are simple consequences of Lemma 5.1:

$\displaystyle\|\mathbf{D}_{t}^{\beta}(w_{2}w_{1})\|_{\mathcal{C}^{\alpha,\frac{\alpha\nu}{2}}(\bar{\Omega}_{T})}+\|\mathbf{D}_{t}^{\nu}(w_{2}w_{1})\|_{\mathcal{C}^{\alpha,\frac{\alpha\nu}{2}}(\bar{\Omega}_{T})}+\|\mathbf{D}_{t}^{\beta}w_{1}\|_{\mathcal{C}^{\alpha,\frac{\alpha\nu}{2}}(\bar{\Omega}_{T})}$
$\displaystyle\leq
C\|\mathbf{D}_{t}^{\nu}w_{1}\|_{\mathcal{C}^{\alpha,\frac{\alpha\nu}{2}}(\bar{\Omega}_{T})},$
$\displaystyle\|\mathbf{D}_{t}^{\beta}(w_{2}w_{1})\|_{\mathcal{C}^{1+\alpha,\frac{(1+\alpha)\nu}{2}}(\partial\Omega_{T})}+\|\mathbf{D}_{t}^{\nu}(w_{2}w_{1})\|_{\mathcal{C}^{1+\alpha,\frac{(1+\alpha)\nu}{2}}(\partial\Omega_{T})}+\|\mathbf{D}_{t}^{\beta}w_{1}\|_{\mathcal{C}^{1+\alpha,\frac{(1+\alpha)\nu}{2}}(\partial\Omega_{T})}$
$\displaystyle\leq
C\big{[}\|\mathbf{D}_{t}^{\nu}w_{1}\|_{\mathcal{C}^{1+\alpha,\frac{(1+\alpha)\nu}{2}}(\partial\Omega_{T})}+\|w_{1}\|_{\mathcal{C}^{2+\alpha,\frac{(2+\alpha)\nu}{2}}(\bar{\Omega}_{T})}\big{]},$

and

$\mathfrak{J}_{\beta}(0;w_{1},w_{2})=\mathfrak{J}_{\nu}(0;w_{1},w_{2})=0.$

Here the positive constant $C$ depends only on $T$, the Lebesgue measure of
$\Omega$, and the norm of $w_{2}$.

We complete this section by discussing the properties of the solution to
initial and initial-boundary value problems for a certain subdiffusion
equation, which will be the key point in the construction of a regularizer to
the linear problems (1.1)-(1.5). To this end, we denote

$\mathbb{R}^{n}_{+}=\\{x:(x_{1},...,x_{n-1})\in\mathbb{R}^{n-1},\,x_{n}>0\\}\qquad\text{and}\qquad\mathbb{R}_{+,T}^{n}=\mathbb{R}^{n}_{+}\times(0,T).$

Let the function $v_{i}=v_{i}(x,t)$ solve the problems

$\begin{cases}\mathbf{D}_{t}^{\nu}v_{1}-\Delta
v_{1}=F_{0}(x,t)\quad\text{in}\quad\mathbb{R}^{n}_{T},\\\
v_{1}(x,0)=v_{1,0}(x)\quad\text{in}\quad\mathbb{R}^{n},\end{cases}$ (5.3)

where $F_{0}$ and $v_{1,0}$ are some given functions; and for $i=2,3,4,$

$\begin{cases}\mathbf{D}_{t}^{\nu}v_{i}-\Delta
v_{i}=0\quad\text{in}\quad\mathbb{R}_{+,T}^{n},\\\
v_{i}(x,0)=0\quad\text{in}\quad\mathbb{R}_{+}^{n},\\\ v_{i}(x,t)\to
0\quad\text{if}\quad|x|\to+\infty,\end{cases}$ (5.4)

with one of the following boundary conditions:

$\displaystyle v_{2}(x,t)$
$\displaystyle=F_{1}(x,t)\quad\text{on}\quad\partial\mathbb{R}_{+,T}^{n},$
(5.5) $\displaystyle\sum_{i=1}^{n}c_{i}\frac{\partial v_{3}}{\partial x_{i}}$
$\displaystyle=F_{2}(x,t)\quad\text{on}\quad\partial\mathbb{R}_{+,T}^{n},$
(5.6)
$\displaystyle\mathbf{D}_{t}^{\nu}v_{4}-\sum_{i=1}^{n}c_{i}\frac{\partial
v_{4}}{\partial x_{i}}$
$\displaystyle=F_{3}(x,t)\quad\text{on}\quad\partial\mathbb{R}_{+,T}^{n},$
(5.7)

where $F_{i}$ are given functions, and $c_{1},\ldots,c_{n}$ are constants with
$c_{n}\neq 0$. The classical solvability of problems (5.3)-(5.7) with
$\nu\in(0,1)$ has been discussed in the one-dimensional case in [25, 26], and
in the multi-dimensional case in [19, 20]. As for $\nu=1$, these problems are
analyzed in [28, Section 4]. We subsume these results in a lemma.

######  Lemma 5.5.

Let $c_{n}\neq 0$, or $c_{n}>0$ in the case of the fractional dynamic boundary
condition (5.7), let $v_{1,0}\in\mathcal{C}^{2+\alpha}(\mathbb{R}^{n})$, and
let

$F_{0}\in\mathcal{C}^{\alpha,\alpha\nu/2}(\bar{\mathbb{R}}_{T}^{n}),\quad
F_{1}\in\mathcal{C}_{0}^{2+\alpha,\frac{2+\alpha}{2}\nu}(\partial\mathbb{R}_{+,T}^{n}),\quad
F_{2},F_{3}\in\mathcal{C}_{0}^{1+\alpha,\frac{1+\alpha}{2}\nu}(\partial\mathbb{R}_{+,T}^{n}).$

Assume also that there exists a positive number $r_{0}$ such that

$v_{1,0}(x),F_{0}(x,t),F_{i}(x,t)\equiv 0,\quad\text{if}\quad|x|>r_{0},\quad
t\in[0,T].$

Then, there are unique classical solutions $v_{i}(x,t)$ to problems
(5.3)-(5.7). In addition, the following estimates hold:

$\displaystyle\|v_{1}\|_{\mathcal{C}^{2+\alpha,\frac{2+\alpha}{2}\nu}(\bar{\mathbb{R}}_{T}^{n})}$
$\displaystyle\leq
C[\|v_{1,0}\|_{\mathcal{C}^{2+\alpha}(\bar{\mathbb{R}}^{n})}+\|F_{0}\|_{\mathcal{C}^{\alpha,\frac{\alpha}{2}\nu}(\bar{\mathbb{R}}_{T}^{n})}],$
$\displaystyle\|v_{2}\|_{\mathcal{C}^{2+\alpha,\frac{2+\alpha}{2}\nu}(\bar{\mathbb{R}}_{T}^{n})}$
$\displaystyle\leq
C\|F_{1}\|_{\mathcal{C}^{2+\alpha,\frac{2+\alpha}{2}\nu}(\partial\mathbb{R}_{+,T}^{n})},$
$\displaystyle\|v_{3}\|_{\mathcal{C}^{2+\alpha,\frac{2+\alpha}{2}\nu}(\bar{\mathbb{R}}_{T}^{n})}$
$\displaystyle\leq
C\|F_{2}\|_{\mathcal{C}^{1+\alpha,\frac{1+\alpha}{2}\nu}(\partial\mathbb{R}_{+,T}^{n})},$
$\displaystyle\|v_{4}\|_{\mathcal{C}^{2+\alpha,\frac{2+\alpha}{2}\nu}(\bar{\mathbb{R}}_{T}^{n})}$
$\displaystyle+\|\mathbf{D}_{t}^{\nu}v_{4}\|_{\mathcal{C}^{1+\alpha,\frac{1+\alpha}{2}\nu}(\partial\mathbb{R}_{+,T}^{n})}\leq
C\|F_{3}\|_{\mathcal{C}^{1+\alpha,\frac{1+\alpha}{2}\nu}(\partial\mathbb{R}_{+,T}^{n})}.$

Here the generic constant $C$ is independent of the right-hand sides in
(5.3)-(5.7).

##  6\. Proof of Theorem 4.1

The strategy of the proof is based on the construction of a regularizer (see
[28, Section 4]), and it consists of fourth main steps. In the first one, we
build a special covering of the domain $\Omega$. Next, assuming the additional
hypotheses on the right-hand sides

$u_{0}(x)=0,\,\,\,x\in\bar{\Omega},\qquad
f(x,0)=0,\,\,\,x\in\bar{\Omega},\qquad\psi_{i}(x,0)=0,\,\,\,x\in\partial\Omega,$
(6.1)

which, in particular, give

$u(x,0)=0,\quad x\in\bar{\Omega},$ (6.2)

we freeze the coefficients of the operators $\mathcal{L}_{1}$ and
$\mathcal{M}_{1}$, and, by exploiting the properties of the solutions of the
so-called model problems (5.3)-(5.7), we construct a regularizer, i.e., the
inverse operator of (1.1)-(1.5) in the case of a small time interval
$t\in[0,\tau]$. After that, we discuss how to extend the obtained results to
the whole time interval $[0,T]$. Finally, we show how to reduce (1.1)-(1.5) in
the general case to the special one related with assumption (6.1), in other
words, we discuss the reduction of problems (1.1)-(1.5) to the problems with
homogenous initial data (6.2). In our analysis, we focus on the case
$\nu_{1}\in(0,1)$, whereas the case $\nu_{1}=1$ is examined either with
similar or simpler arguments, due to the equivalent definitions of Caputo
fractional derivatives.

###  6.1. Step I: Covering of the domain $\Omega$.

For an arbitrarily fixed $\lambda>0$, it is always possible to find a finite
collection of points $x^{m}\in\bar{\Omega}$ along with sets

$\omega^{m}=B_{\lambda/2}(x^{m})\cap\bar{\Omega}\qquad\text{and}\qquad\Omega^{m}=B_{\lambda}(x^{m})\cap\bar{\Omega},$

satisfying the following properties:

  * (i)

$\bigcup_{m}\omega^{m}=\bigcup_{m}\Omega^{m}=\bar{\Omega}$;

  * (ii)

there exists a number $\mathcal{N}_{0}$, independent of $\lambda$, such that
the intersection of any $\mathcal{N}_{0}+1$ distinct $\Omega^{m}$ (and
consequently any $\mathcal{N}_{0}+1$ distinct $\omega^{m}$) is empty.

Notice that, by construction,

$x^{m}\in\omega^{m}\subset\overline{\omega^{m}}\subset\Omega^{m}\subset\bar{\Omega}.$

Moreover, we partition the sets of indexes $m$ into the disjoint union
$\mathfrak{M}\cup\mathfrak{N}$, by setting

$m\in\mathfrak{M}\quad\text{if}\quad\overline{\Omega^{m}}\cap\partial\Omega=\emptyset\qquad\text{and}\qquad
m\in\mathfrak{N}\quad\text{if}\quad\overline{\omega^{m}}\cap\partial\Omega\neq\emptyset.$

In the sequel, let us denote $\partial\Omega^{m}=\partial\Omega\cap
B_{\lambda}(x^{m})$. Let $\xi^{m}=\xi^{m}(x):\Omega\to[0,1]$ be a smooth
function possessing the following properties: $\xi^{m}\in(0,1)$ if
$x\in\Omega^{m}\backslash\omega^{m}$ and

$\xi^{m}=\begin{cases}1,\quad\text{if}\quad x\in\overline{\omega^{m}},\\\
0,\quad\text{if}\quad
x\in\bar{\Omega}\backslash\overline{\Omega^{m}},\end{cases}\qquad|D_{x}^{j}\xi^{m}|\leq
C\lambda^{-|j|},\,|j|\geq 1,\qquad
1\leq\sum_{m}(\xi^{m})^{2}\leq\mathcal{N}_{0}.$

Then, we define the function

$\eta^{m}=\frac{\xi^{m}}{\sum_{j}(\xi^{j})^{2}}.$ (6.3)

Due to the properties of the function $\xi^{m}$, we see that $\eta^{m}$
vanishes for $x\in\bar{\Omega}\backslash\overline{\Omega^{m}}$, and
$|D_{x}^{j}\eta^{m}|\leq C\lambda^{-|j|}$. Thus, the product $\eta^{m}\xi^{m}$
defines a partition of unity via the formula

$\sum_{m}\eta^{m}\xi^{m}=1.$

At this point, we define the local coordinate systems connected with each
point $x^{m}$, $m\in\mathfrak{N}$. For each $m\in\mathfrak{N}$, the point
$x^{m}$ will be the origin of a local coordinate system. Let $\partial\Omega$
be described by $y_{n}=\mathfrak{F}^{m}(y_{1},...,y_{n-1})$ in a small
vicinity of each point $x^{m},$ $m\in\mathfrak{N}$, and

$y=\mathfrak{B}^{(m)}(x-x^{m}),\qquad\bigg{|}\frac{\partial\mathfrak{F}^{m}}{\partial
y_{i}}\bigg{|}\leq C\lambda,\,i=1,2,...,n-1,$

where $\mathfrak{B}^{(m)}=(\mathfrak{b}_{ij}^{m})_{i,j=1,..,n}$ is an
orthogonal matrix with elements $\mathfrak{b}_{ij}^{m}$, and
$(\mathfrak{b}_{ij}^{m})^{-1}$ is an element of the inverse matrix to
$\mathfrak{B}^{(m)}$. To obtain the local “flatness” of the boundary, we make
the change of variables

$z_{i}=y_{i},\quad z_{n}=y_{n}-\mathfrak{F}^{m}(y_{1},...,y_{n-1}),\quad
i=1,2,...,n-1,\quad m\in\mathfrak{N}.$

Hence, we have built the mapping $Z_{m}$ which connects the original variable
$x=(x_{1},...,x_{n})$ with the new variable $z=(z_{1},...,z_{n})$ in a
neighborhood of each point $x^{m},$ $m\in\mathfrak{N}$ via relations:

$x=Z_{m}(z)\qquad\text{and}\qquad z=Z^{-1}_{m}(x).$

Next, we introduce the following norms in the spaces
$\mathcal{C}_{0}^{l+\alpha,\frac{(l+\alpha)\nu_{1}}{2}}(\bar{\Omega}_{T}),$
$l=0,1,2,$ which are related with the covering $\\{\Omega^{m}\\}$:

$\\{v\\}_{\mathcal{C}^{l+\alpha,\frac{(l+\alpha)\nu_{1}}{2}}(\bar{\Omega}_{T})}=\sup_{m}\|v\|_{\mathcal{C}^{l+\alpha,\frac{(l+\alpha)\nu_{1}}{2}}(\bar{\Omega}^{m}_{T})}.$

We now state a lemma, which subsumes Propositions 4.5-4.7 in our previous work
[21], in order to describe the properties of these norms. To this end, for an
arbitrarily given $0<\kappa<1$, we define

$\tau=\lambda^{2/\nu_{1}}\kappa,$ (6.4)

such that $\tau\in(0,T]$. Then we consider (any) function $\Phi_{m}(x)$
(defined in $\Omega^{m}$) such that

$|D_{x}^{j}\Phi_{m}(x)|\leq C\lambda^{-|j|},\quad 0\leq|j|\leq 2,$

along with (any) function $\tilde{v}(x,t)$ of the form

$\tilde{v}(x,t)=\sum_{m\in\mathfrak{M}\,\bigcup\mathfrak{N}}v^{m}(x,t),$

for some
$v^{m}\in\mathcal{C}^{l+\alpha,\frac{l+\alpha}{2}\nu_{1}}(\bar{\Omega}_{\tau}^{m}),$
$l=0,1,2,$ with $v^{m}$ vanishing outside $\Omega^{m}$.

######  Lemma 6.1.

Let (6.4) hold. Then for any
$v\in\mathcal{C}_{0}^{l+\alpha,\frac{l+\alpha}{2}\nu_{1}}(\bar{\Omega}_{\tau})$,
$l=0,1,2,$ we have the following relations:

$\displaystyle\\{v\\}_{\mathcal{C}^{l+\alpha,\frac{l+\alpha}{2}\nu_{1}}(\bar{\Omega}_{\tau})}$
$\displaystyle\leq\|v\|_{\mathcal{C}^{l+\alpha,\frac{l+\alpha}{2}\nu_{1}}(\bar{\Omega}_{\tau})}\leq
C\\{v\\}_{\mathcal{C}^{l+\alpha,\frac{l+\alpha}{2}\nu_{1}}(\bar{\Omega}_{\tau})},$
$\displaystyle\|\Phi_{m}v\|_{\mathcal{C}^{l+\alpha,\frac{l+\alpha}{2}\nu_{1}}(\overline{\Omega^{m}}_{\tau})}$
$\displaystyle\leq
C\|v\|_{\mathcal{C}^{l+\alpha,\frac{l+\alpha}{2}\nu_{1}}(\overline{\Omega^{m}}_{\tau})},$
$\displaystyle\\{\tilde{v}\\}_{\mathcal{C}^{l+\alpha,\frac{l+\alpha}{2}\nu_{1}}(\bar{\Omega}_{\tau})}$
$\displaystyle\leq
C\underset{m\in\mathfrak{M}\,\bigcup\mathfrak{N}}{\sup}\|v^{m}\|_{\mathcal{C}^{l+\alpha,\frac{l+\alpha}{2}\nu_{1}}(\overline{\Omega^{m}}_{\tau})}.$

Here the positive constant $C$ is independent of $\lambda$ and $\tau$.

###  6.2. Step II: Construction of a regularizer for (6.1).

We aim to construct the inverse operator for problem (1.1), (6.2) and (1.5),
i.e., in FDBC case. The analysis of the remaining cases (1.3) and (1.4) are
performed in similar manner. First, we recall that assumption H4, H5 and (6.1)
imply

$f(x,t)\in\mathcal{C}_{0}^{\alpha,\frac{\alpha\nu_{1}}{2}}(\bar{\Omega}_{T}),\quad\psi_{3}\in\mathcal{C}_{0}^{1+\alpha,\frac{1+\alpha}{2}\nu_{1}}(\partial\Omega_{T}).$
(6.5)

For the sake of convenience, we rewrite problem (1.1), (6.2) and (1.5) in the
compact form

$\mathbb{L}u=\mathbf{F},\qquad\mathbf{F}=\\{f,\psi_{3}\\}.$ (6.6)

Here, $\mathbb{L}$ is the linear operator acting as

$\mathbb{L}u=\\{\mathcal{A}u,\mathcal{A}_{1}u|_{\partial\Omega_{\tau}}\\},$

where $\mathcal{A}$ is the left-hand sides of (1.1), while $\mathcal{A}_{1}$
is the left-hand side of (1.5). For $m\in\mathfrak{M}\cup\mathfrak{N}$, we set

$a_{ij}^{m}=a_{ij}(x^{m},0),\qquad\varrho_{1}^{m}=\varrho_{1}(x^{m},0),\qquad
c_{i}^{m}=c_{i}(x^{m},0),\qquad f^{m}(x,t)=\xi^{m}(x)f(x,t),$

and, for $m\in\mathfrak{N}$,

$\tilde{f}^{m}(z,t)=f^{m}(x,t)\big{|}_{x=Z_{m}(z)},\qquad\tilde{\psi}^{m}(z,t)=\xi^{m}(x)\psi_{3}(x,t)\big{|}_{x=Z_{m}(z)},$

with $\xi^{m},$ $Z_{m}(z)$, $\mathfrak{M}$, $\mathfrak{N}$ as in Subsection
6.1. For $m\in\mathfrak{M}\cup\mathfrak{N}$, $\tau\in(0,T]$ and $\lambda$ as
in (6.4), we define the functions $u^{m}(x,t)$ to be the solutions to the
following problems: if $m\in\mathfrak{M},$ then

$\begin{cases}\displaystyle\varrho_{1}^{m}\mathbf{D}_{t}^{\nu_{1}}u^{m}-\sum_{ij=1}^{n}a_{ij}^{m}\frac{\partial^{2}u^{m}}{\partial
x_{i}\partial x_{j}}=f^{m}(x,t)\quad\text{in}\ \mathbb{R}^{n}_{\tau},\\\
u^{m}(x,0)=0\quad\text{in }\mathbb{R}^{n},\end{cases}$ (6.7)

while, for $m\in\mathfrak{N}$,

$u^{m}(x,t)=\tilde{u}^{m}(z,t)\big{|}_{z=Z_{m}^{-1}(x)},$

where $\tilde{u}^{m}$ solves the initial-boundary value problem

$\begin{cases}\displaystyle\varrho_{1}^{m}\mathbf{D}_{t}^{\nu_{1}}\tilde{u}^{m}-\sum_{ij=1}^{n}a_{ij}^{m}\frac{\partial^{2}\tilde{u}^{m}}{\partial
z_{i}\partial z_{j}}=\tilde{f}^{m}(z,t)\quad\text{in
}\mathbb{R}^{n}_{+,\tau},\\\
\displaystyle\varrho_{1}^{m}\mathbf{D}_{t}^{\nu_{1}}\tilde{u}^{m}-\sum\limits_{i=1}^{n}c_{i}^{m}\frac{\partial\tilde{u}^{m}}{\partial
z_{i}}=\tilde{\psi}^{m}(z,t)\quad\text{on }\partial\mathbb{R}^{n-1}_{\tau},\\\
\vskip 8.53581pt\cr\tilde{u}^{m}(z,0)=0\quad\text{in
}\mathbb{R}^{n}_{+}.\end{cases}$ (6.8)

At this point, we define the space

$\mathcal{H}=\big{\\{}u:\,u\in\mathcal{C}_{0}^{2+\alpha,(2+\alpha)\nu_{1}/2}(\bar{\Omega}_{\tau}),\,\,\mathbf{D}^{\nu_{1}}_{t}u\in\mathcal{C}_{0}^{1+\alpha,(1+\alpha)\nu_{1}/2}(\partial\Omega_{\tau})\big{\\}},$

normed by

$\|u\|_{\mathcal{H}}=\|u\|_{\mathcal{C}_{0}^{2+\alpha,(2+\alpha)\nu_{1}/2}(\bar{\Omega}_{\tau})}+\|\mathbf{D}^{\nu_{1}}_{t}u\|_{\mathcal{C}_{0}^{1+\alpha,(1+\alpha)\nu_{1}/2}(\partial\Omega_{\tau})},$

together with the product space

$\mathcal{H}_{0}=\mathcal{C}_{0}^{\alpha,\alpha\nu_{1}/2}(\bar{\Omega}_{\tau})\times\mathcal{C}_{0}^{1+\alpha,(1+\alpha)\nu_{1}/2}(\partial\Omega_{\tau}),$

normed by

$\|(f,\psi_{3})\|_{\mathcal{H}_{0}}=\|f\|_{\mathcal{C}_{0}^{\alpha,\alpha\nu_{1}/2}(\bar{\Omega}_{\tau})}+\|\psi_{3}\|_{\mathcal{C}_{0}^{1+\alpha,(1+\alpha)\nu_{1}/2}(\partial\Omega_{\tau})}.$

We are now in the position to give the definition of a regularizer.

######  Definition 6.2.

Let $\tau\in(0,T]$. An operator $\mathfrak{R}:\mathcal{H}_{0}\to\mathcal{H}$
is called a regularizer on the time-interval $[0,\tau]$, if

$\mathfrak{R}(f,\psi_{3})=\sum_{m\in\mathfrak{M}\cup\mathfrak{N}}\eta^{m}(x)u^{m}(x,t),$

where the functions $\eta^{m}(x)$ and $u^{m}(x,t)$ are defined in (6.3) and
(6.7)-(6.8), respectively.

The following result details the main properties of $\mathfrak{R}$, allowing
us eventually to construct the inverse of $\mathbb{L}$.

######  Lemma 6.3.

Let $\tau\in(0,T]$ satisfy (6.4). We assume that the hypotheses of Theorem 4.1
and (6.1) hold. Then, for any $\mathbf{F}\in\mathcal{H}_{0}$ and
$u\in\mathcal{H}$ the following hold:

  * (i)

$\mathfrak{R}$ is a bounded operator:

$\|\mathfrak{R}\mathbf{F}\|_{\mathcal{H}}\leq
C\|\mathbf{F}\|_{\mathcal{H}_{0}},$ (6.9)

where the positive constant $C$ is independent of $\lambda$ and $\tau$.

  * (ii)

There exist operators $\mathfrak{T}_{1}:\mathcal{H}_{0}\to\mathcal{H}_{0}$ and
$\mathfrak{T}_{2}:\mathcal{H}\to\mathcal{H}$ such that the decompositions

$\mathbb{L}\mathfrak{R}\mathbf{F}=\mathbf{F}+\mathfrak{T}_{1}\mathbf{F}\qquad\text{and}\qquad\mathfrak{R}\mathbb{L}u=u+\mathfrak{T}_{2}u$

hold, and

$\|\mathfrak{T}_{1}\mathbf{F}\|_{\mathcal{H}_{0}}\leq\frac{1}{2}\|\mathbf{F}\|_{\mathcal{H}_{0}}\qquad\text{and}\qquad\|\mathfrak{T}_{2}u\|_{\mathcal{H}}\leq\frac{1}{2}\|u\|_{\mathcal{H}}.$

###### Proof.

It is worth noting that the results of Lemma 5.5 are valid in the case of
problems (6.7) and (6.8). Then, collecting [21, Proposition 4.4] with Lemmas
5.1-5.5, 6.1 and Remark 5.4, we end up with the estimates:

$\displaystyle\|\mathfrak{R}\mathbf{F}\|_{\mathcal{H}}$ $\displaystyle\leq
C\Big{[}\underset{m\in\mathfrak{N}\cup\mathfrak{M}}{\sup}\|u^{m}\|_{\mathcal{C}^{2+\alpha,\frac{2+\alpha}{2}\nu_{1}}(\bar{\Omega}_{\tau}^{m})}+\underset{m\in\mathfrak{N}}{\sup}\|\mathbf{D}_{t}^{\nu_{1}}u^{m}\|_{\mathcal{C}^{1+\alpha,\frac{1+\alpha}{2}\nu_{1}}(\partial\Omega_{\tau}^{m})}\Big{]}$
$\displaystyle\leq
C\Big{[}\underset{m\in\mathfrak{N}\cup\mathfrak{M}}{\sup}\|f\xi^{m}\|_{\mathcal{C}^{\alpha,\frac{\alpha}{2}\nu_{1}}(\bar{\Omega}_{\tau}^{m})}+\underset{m\in\mathfrak{N}}{\sup}\|\psi_{3}\xi^{m}\|_{\mathcal{C}^{1+\alpha,\frac{1+\alpha}{2}\nu_{1}}(\partial\Omega_{\tau}^{m})}\Big{]}$
$\displaystyle\leq C\|\mathbf{F}\|_{\mathcal{H}_{0}},$

where $C$ is independent of $\lambda$ and $\tau$. The last inequality is just
(i). Now we verify (ii). Here, we limit ourselves to deal with
$\mathfrak{T}_{1}$, being the other case completely analogous. The definition
of the operator $\mathbb{L}$ together with (6.3) allow us to conclude that

$\mathbb{L}\mathfrak{R}\mathbf{F}=\\{\mathcal{A}\mathfrak{R}\mathbf{F},\mathcal{A}_{1}\mathfrak{R}\mathbf{F}|_{\partial\Omega_{\tau}}\\},$

with

$\mathcal{A}\mathfrak{R}\mathbf{F}=\mathcal{A}^{0}\mathfrak{R}\mathbf{F}+\mathcal{A}^{1}\mathfrak{R}\mathbf{F}\qquad\text{and}\qquad\mathcal{A}_{1}\mathfrak{R}\mathbf{F}|_{\partial\Omega_{\tau}}=\mathcal{A}_{1}^{0}\mathfrak{R}\mathbf{F}+\mathcal{A}_{1}^{1}\mathfrak{R}\mathbf{F},$

where we set

$\mathcal{A}^{0}\mathfrak{R}\mathbf{F}=\begin{cases}\sum_{m}\varrho_{1}^{m}\mathbf{D}_{t}^{\nu_{1}}u^{m}\eta^{m}(x)-\mathcal{L}_{1}\mathfrak{R}\mathbf{F}-\mathcal{K}*\mathcal{L}_{2}\mathfrak{R}\mathbf{F},\quad
m\in\mathfrak{M},\\\
\sum_{m}\varrho_{1}^{m}\eta^{m}(x)\mathbf{D}_{t}^{\nu_{1}}\tilde{u}^{m}(z,t)|_{z=Z_{m}^{-1}(x)}-\mathcal{L}_{1}\mathfrak{R}\mathbf{F}-\mathcal{K}*\mathcal{L}_{2}\mathfrak{R}\mathbf{F},\quad
m\in\mathfrak{N},\end{cases}$
$\mathcal{A}^{1}\mathfrak{R}\mathbf{F}=\begin{cases}\mathbf{D}_{t}^{\nu_{1}}(\varrho_{1}\mathfrak{R}\mathbf{F})-\mathbf{D}_{t}^{\nu_{2}}(\varrho_{2}\mathfrak{R}\mathbf{F})-\sum_{m}\varrho_{1}^{m}\mathbf{D}_{t}^{\nu_{1}}u^{m}\eta^{m}(x),\quad
m\in\mathfrak{M},\\\ \\\
\mathbf{D}_{t}^{\nu_{1}}(\varrho_{1}\mathfrak{R}\mathbf{F})-\mathbf{D}_{t}^{\nu_{2}}(\varrho_{2}\mathfrak{R}\mathbf{F})-\sum_{m}\varrho_{1}^{m}\eta^{m}(x)\mathbf{D}_{t}^{\nu_{1}}\tilde{u}^{m}(z,t)|_{z=Z_{m}^{-1}(x)},\quad
m\in\mathfrak{N},\end{cases}$
$\displaystyle\mathcal{A}^{0}_{1}\mathfrak{R}\mathbf{F}$
$\displaystyle=\bigg{\\{}\sum_{m\in\mathfrak{N}}\varrho_{1}^{m}\mathbf{D}_{t}^{\nu_{1}}\tilde{u}^{m}(z,t)|_{z=Z_{m}^{-1}(x)}-\mathcal{M}_{1}\mathfrak{R}\mathbf{F}+\mathcal{K}_{0}*\mathcal{M}_{2}\mathfrak{R}\mathbf{F}\bigg{\\}}\bigg{|}_{\partial\Omega_{\tau}}$
$\displaystyle\mathcal{A}^{1}_{1}\mathfrak{R}\mathbf{F}$
$\displaystyle=\bigg{\\{}\mathbf{D}_{t}^{\nu_{1}}(\varrho_{1}\mathfrak{R}\mathbf{F})-\sum_{m\in\mathfrak{N}}\varrho_{1}^{m}\mathbf{D}_{t}^{\nu_{1}}\tilde{u}^{m}(z,t)|_{z=Z_{m}^{-1}(x)}-\mathbf{D}_{t}^{\nu_{2}}(\varrho_{2}\mathfrak{R}\mathbf{F})\bigg{\\}}\bigg{|}_{\partial\Omega_{\tau}},$

Then, Lemma 5.2 in [21] and Theorem 2 in [20] tell us that

$\mathcal{A}^{0}\mathfrak{R}\mathbf{F}=f+\mathfrak{T}_{1}^{1}\mathfrak{R}\mathbf{F}\qquad\text{and}\qquad\mathcal{A}_{1}^{0}\mathfrak{R}\mathbf{F}=\psi_{3}+\mathfrak{T}_{1}^{2}\mathfrak{R}\mathbf{F},$

where

$\|\mathfrak{T}_{1}^{1}\mathfrak{R}\mathbf{F}\|_{\mathcal{C}^{\alpha,\alpha\nu_{1}/2}(\bar{\Omega}_{\tau})}\leq\frac{1}{8}\|\mathbf{F}\|_{\mathcal{H}_{0}}\qquad\text{and}\qquad\|\mathfrak{T}_{1}^{2}\mathfrak{R}\mathbf{F}\|_{\mathcal{C}^{1+\alpha,\frac{1+\alpha}{2}\nu_{1}}(\partial\Omega_{\tau})}\leq\frac{1}{8}\|\mathbf{F}\|_{\mathcal{H}_{0}},$
(6.10)

provided that $\lambda$ and $\tau$ comply with (6.4). Hence, we are left to
prove the estimates

$\|\mathcal{A}^{1}\mathfrak{R}\mathbf{F}\|_{\mathcal{C}^{\alpha,\alpha\nu_{1}/2}(\bar{\Omega}_{\tau})}\leq\frac{1}{8}\|\mathbf{F}\|_{\mathcal{H}_{0}}\qquad\text{and}\qquad\|\mathcal{A}^{1}_{1}\mathfrak{R}\mathbf{F}\|_{\mathcal{C}^{1+\alpha,\frac{1+\alpha}{2}\nu_{1}}(\partial\Omega_{\tau})}\leq\frac{1}{8}\|\mathbf{F}\|_{\mathcal{H}_{0}}.$
(6.11)

Indeed, point (ii) for $\mathfrak{T}_{1}$ immediately follows from
representation of $\mathbb{L}\mathfrak{R}\mathbf{F}$ and estimates
(6.10)-(6.11), implying that

$\mathfrak{T}_{1}\mathfrak{R}\mathbf{F}=\\{\mathfrak{T}_{1}^{1}\mathfrak{R}\mathbf{F}+\mathcal{A}^{1}\mathfrak{R}\mathbf{F},\mathfrak{T}_{1}^{2}\mathfrak{R}\mathbf{F}+\mathcal{A}^{1}_{1}\mathfrak{R}\mathbf{F}\\}\qquad\text{and}\qquad\|\mathfrak{T}_{1}\mathfrak{R}\mathbf{F}\|_{\mathcal{H}_{0}}\leq\frac{1}{2}\|\mathbf{F}\|_{\mathcal{H}_{0}}.$

Concerning the first inequality in (6.11), we treat the case
$m\in\mathfrak{M}$ (the case $m\in\mathfrak{N}$ being similar). Appealing to
Corollary 3.1 in [23], and keeping in mind that we have null initial data, we
have

$\displaystyle\mathbf{D}_{t}^{\nu_{1}}(\varrho_{1}\mathfrak{R}\mathbf{F})-\mathbf{D}_{t}^{\nu_{2}}(\varrho_{2}\mathfrak{R}\mathbf{F})-\sum_{m\in\mathfrak{M}}\varrho_{1}^{m}\mathbf{D}_{t}^{\nu_{1}}u^{m}(x,t)\eta^{m}$
$\displaystyle=\sum_{m\in\mathfrak{M}}\\{[\varrho_{1}-\varrho_{1}^{m}]\eta^{m}\mathbf{D}_{t}^{\nu_{1}}u^{m}(x,t)+\frac{\nu}{\Gamma(1-\nu_{1})}\eta^{m}\mathfrak{J}_{\nu_{1}}(t;u^{m},\varrho_{1})$
$\displaystyle\quad-\frac{\nu_{2}}{\Gamma(1-\nu_{2})}\eta^{m}\mathfrak{J}_{\nu_{2}}(t;u^{m},\varrho_{2})-\varrho_{2}\eta^{m}\mathbf{D}_{t}^{\nu_{2}}u^{m}(x,t)\\}.$

On account of the properties of the functions $\varrho_{1}$ and $\varrho_{2}$
(see H3), and exploiting Lemmas 5.1 and 6.1 along with Remark 5.3 to evaluate
the terms in the right-hand sides of the equality above, we conclude that

$\displaystyle\|\mathbf{D}_{t}^{\nu_{1}}(\varrho_{1}\mathfrak{R}\mathbf{F})-\mathbf{D}_{t}^{\nu_{2}}(\varrho_{2}\mathfrak{R}\mathbf{F})-\sum_{m\in\mathfrak{M}}\varrho_{1}^{m}\mathbf{D}_{t}^{\nu_{1}}u^{m}(x,t)\eta^{m}\|_{\mathcal{C}^{\alpha,\alpha\nu_{1}/2}(\bar{\Omega}_{\tau})}$
$\displaystyle\leq
C[\tau^{\delta_{0}-\alpha\nu_{1}/2}+\kappa^{\alpha}+\kappa^{\alpha}\tau^{\nu_{1}/2-\nu_{2}}+\tau^{\delta_{1}-\nu_{2}+\nu_{1}(1-\alpha)/2}]\|\mathfrak{R}\mathbf{F}\|_{\mathcal{H}}$
$\displaystyle\leq
C[\tau^{\delta_{0}-\alpha\nu_{1}/2}+\kappa^{\alpha}+\kappa^{\alpha}\tau^{\nu_{1}/2-\nu_{2}}+\tau^{\delta_{1}-\nu_{2}+\nu_{1}(1-\alpha)/2}]\|\mathbf{F}\|_{\mathcal{H}_{0}},$

where $\delta_{0}=\min\\{1,\gamma_{0}\\}$ and
$\delta_{1}=\min\\{1,\gamma_{1}\\}$. The constant $C$ is independent of
$\lambda$ and $\tau$, and depends only on the norms of $\varrho_{1}$,
$\varrho_{2}$, the Lebesgue measure of $\Omega$ and $T$. Thanks to the
relation between $\nu_{1}$ and $\nu_{2}$ (see H1), and assumption H3 on
$\gamma_{0}$ and $\gamma_{1},$ the last two estimates provide the first
inequality in (6.11). The second one follows by recasting the arguments above,
but using Lemma 5.2 in place of Lemma 5.1. ∎

Coming to construction of the inverse of $\mathbb{L}$, we note that Lemma 6.3
ensures the existence of the bounded operators $(I+\mathfrak{T}_{1})^{-1}$ and
$(I+\mathfrak{T}_{2})^{-1}$ ($I$ is the identity in the respective spaces).
Therefore,

$\mathbb{L}\mathfrak{R}(I+\mathfrak{T}_{1})^{-1}\mathbf{F}=\mathbf{F}\qquad\text{and}\qquad(I+\mathfrak{T}_{2})^{-1}\mathfrak{R}\mathbb{L}u=u,$

namely, $\mathbb{L}$ has bounded right and left inverse operators, hence

$\mathfrak{R}(I+\mathfrak{T}_{1})^{-1}=(I+\mathfrak{T}_{2})^{-1}\mathfrak{R}=\mathbb{L}^{-1}:{\mathcal{H}}_{0}\to{\mathcal{H}}.$

Accordingly, the unique solution of (6.6) is given by

$u=\mathbb{L}^{-1}(f,\psi_{3})\quad\text{for}\quad t\in[0,\tau].$ (6.12)

The estimate of the norm $\mathbb{L}^{-1}$ follows from the estimates of Lemma
6.3. In summary, we have verified Theorem 4.1 (in the case of (6.1)) for a
small time interval $[0,\tau]$.

###  6.3. Step III: Extension of the solution to whole interval $[\tau,T]$.

The next goal is to extend the solution found in Step I to the intervals
$[\tau,2\tau],$ $[2\tau,3\tau]$ and so on, so to cover the whole $[\tau,T]$.
Again, we shall give the details only for the (most difficult) case FDBC.
First, we set

$\Phi(x,t)=\begin{cases}\mathbf{D}_{t}^{\nu_{1}}u(x,t)-\Delta
u(x,t),\qquad\quad t\in[0,\tau],\,x\in\bar{\Omega},\\\
[\mathbf{D}_{t}^{\nu_{1}}u(x,t)-\Delta u(x,t)]|_{t=\tau},\quad
t\in[\tau,2\tau],\,x\in\bar{\Omega},\end{cases}$

and

$\Psi(x,t)=\begin{cases}\mathbf{D}_{t}^{\nu_{2}}u(x,t)-\frac{\partial
u}{\partial N},\qquad\quad\,t\in[0,\tau],\,x\in\partial\Omega,\\\
[\mathbf{D}_{t}^{\nu_{2}}u(x,t)-\frac{\partial u}{\partial N}]|_{t=\tau},\quad
t\in[\tau,2\tau],\,x\in\partial\Omega.\end{cases}$

The results of Step II tell us that

$\|\Phi\|_{\mathcal{C}^{\alpha,\frac{\alpha\nu_{1}}{2}}(\bar{\Omega}_{2\tau})}\leq
C\|u\|_{\mathcal{C}^{2+\alpha,\frac{(2+\alpha)\nu_{1}}{2}}(\bar{\Omega}_{\tau})}\leq
C\big{[}\|f\|_{\mathcal{C}^{\alpha,\frac{\alpha\nu_{1}}{2}}(\bar{\Omega}_{T})}+\|\psi_{3}\|_{\mathcal{C}^{1+\alpha,\frac{(1+\alpha)\nu_{1}}{2}}(\partial\Omega_{T})}\big{]},$

and

$\displaystyle\|\Psi\|_{\mathcal{C}^{1+\alpha,\frac{(1+\alpha)\nu_{1}}{2}}(\partial\Omega_{2\tau})}$
$\displaystyle\leq
C[\|u\|_{\mathcal{C}^{2+\alpha,\frac{(2+\alpha)\nu_{1}}{2}}(\bar{\Omega}_{\tau})}+\|\mathbf{D}_{t}^{\nu_{1}}u\|_{\mathcal{C}^{1+\alpha,\frac{(1+\alpha)\nu_{1}}{2}}(\partial\Omega_{\tau})}]$
$\displaystyle\leq
C[\|f\|_{\mathcal{C}^{\alpha,\frac{\alpha\nu_{1}}{2}}(\bar{\Omega}_{T})}+\|\psi_{3}\|_{\mathcal{C}^{1+\alpha,\frac{(1+\alpha)\nu_{1}}{2}}(\partial\Omega_{T})}].$

After that, we define the function $v=v(x,t)$ to be the solution of the
initial-boundary value problem

$\begin{cases}\mathbf{D}_{t}^{\nu_{1}}v-\Delta
v=\Phi\quad\text{in}\quad\Omega_{2\tau},\\\
\mathbf{D}_{t}^{\nu_{1}}v-\frac{\partial v}{\partial
N}=\Psi\quad\text{on}\quad\partial\Omega_{2\tau},\\\
v(x,0)=0\quad\,\text{in}\quad\bar{\Omega}.\end{cases}$ (6.13)

In light of the regularity of the right-hand side in (6.13) and the
compatibility conditions H5, together with the requirement (6.1), we can apply
Theorem 2 in [19] to (6.13), so to get the existence of a unique classical
solution $v$ satisfying the properties:

$v\in\mathcal{C}^{2+\alpha,\frac{2+\alpha}{2}\nu_{1}}(\bar{\Omega}_{2\tau}),\qquad\mathbf{D}_{t}^{\nu_{1}}v\in\mathcal{C}^{1+\alpha,\frac{1+\alpha}{2}\nu_{1}}(\partial\Omega_{2\tau}),$

and

$v(x,t)=u(x,t)\quad\text{for}\quad(x,t)\in\bar{\Omega}_{\tau}.$

Now we are ready to look for the solution of (1.1), (6.2), (1.5) for
$t\in[0,2\tau]$ in the form

$u(x,t)=U(x,t)+v(x,t),$

where the unknown function $U(x,t)$ solves the problem

$\begin{cases}\mathbf{D}_{t}^{\nu_{1}}(\varrho_{1}U)-\mathbf{D}_{t}^{\nu_{2}}(\varrho_{2}U)-\mathcal{L}_{1}U-\mathcal{K}*\mathcal{L}_{2}U=f^{\star}\qquad\text{in}\quad\Omega_{2\tau},\\\
\vskip
5.69054pt\cr\mathbf{D}_{t}^{\nu_{1}}(\varrho_{1}U)-\mathbf{D}_{t}^{\nu_{2}}(\varrho_{2}U)-\mathcal{M}_{1}U+\mathcal{K}_{0}\star\mathcal{M}_{2}U=\psi^{\star}\quad\text{on}\quad\partial\Omega_{2\tau},\\\
\vskip 5.69054pt\cr U(x,0)=0\quad\text{in}\quad\bar{\Omega}.\end{cases}$
(6.14)

Here we set

$\displaystyle f^{\star}$
$\displaystyle=f-\mathbf{D}_{t}^{\nu_{1}}(\varrho_{1}v)+\mathbf{D}_{t}^{\nu_{2}}(\varrho_{2}v)+\mathcal{L}_{1}v+\mathcal{K}*\mathcal{L}_{2}v,$
$\displaystyle\psi^{\star}$
$\displaystyle=\psi_{3}-\mathbf{D}_{t}^{\nu_{1}}(\varrho_{1}v)+\mathbf{D}_{t}^{\nu_{2}}(\varrho_{2}v)+\mathcal{M}_{1}v-\mathcal{K}_{0}*\mathcal{M}_{2}v.$

Collecting properties of $v$ with assumptions H3, H4, (6.1), and exploiting
Remark 5.4 and [21, Lemma 4.1], we arrive at the relations:

$\psi^{\star}\in\mathcal{C}^{1+\alpha,\frac{1+\alpha}{2}\nu_{1}}(\partial\Omega_{2\tau}),\qquad
f^{\star}\in\mathcal{C}^{\alpha,\frac{\nu_{1}\alpha}{2}}(\bar{\Omega}_{2\tau}),$
(6.15)

and

$\psi^{\star}\equiv 0\quad x\in\partial\Omega,\,\,\qquad f^{\star}\equiv
0\quad x\in\bar{\Omega},\,\,\,t\in[0,\tau].$ (6.16)

In particular, appealing to the results of Step II, (6.16) tell us that

$U(x,t)=0\quad\text{for}\quad x\in\bar{\Omega}_{\tau}.$ (6.17)

Finally, let us introduce the new time-variable

$\sigma=t-\tau\in[-\tau,\tau]$

in problem (6.14), and for every function $\zeta$ appearing in the sequel we
denote

$\bar{\zeta}(x,\sigma)=\zeta(x,\sigma+\tau).$

and we call $\bar{\mathcal{L}}_{i}$ and $\bar{\mathcal{M}}_{i}$ the operators
${\mathcal{L}}_{i}$ and ${\mathcal{M}}_{i}$, respectively, with the _bar_
coefficients. It is easy to verify that the coefficients of
$\bar{\mathcal{L}}_{i},$ $\bar{\mathcal{M}}_{i},$ and the functions
$\bar{\psi},$ $\bar{f},$ $\bar{\varrho}_{1},$ $\bar{\varrho}_{2}$ meet the
requirements of Theorem 4.1. Besides, relations (6.15)-(6.17) provide

$\bar{U}=\bar{\psi}=\bar{f}=0\quad\text{if}\quad\sigma\in[-\tau,0],$

and

$\mathbf{D}_{\sigma}^{\nu_{1}}\bar{U}=\mathbf{D}_{t}^{\nu_{1}}U,\qquad\mathbf{D}_{\sigma}^{\nu_{2}}\bar{U}=\mathbf{D}_{t}^{\nu_{2}}U,\quad\text{if}\quad\sigma\in[-\tau,\tau],\,t\in[0,\tau].$

It is worth noting that the latter two equalities above are examined in [25,
(3.111)]. Moreover, recasting the arguments in [21, p.441], we conclude that

$(\mathcal{K}*\mathcal{L}_{2}U)(x,t)=(\mathcal{K}*\bar{\mathcal{L}_{2}}\bar{U})(x,\sigma)\qquad\text{and}\qquad(\mathcal{K}_{0}*\mathcal{M}_{2}U)(x,t)=(\mathcal{K}_{0}*\bar{\mathcal{M}_{2}}\bar{U})(x,\sigma).$

In order to rewrite problem (6.14) in the new variable, we are left to
recalculate the terms: $\mathbf{D}_{t}^{\nu_{1}}(\varrho_{1}U)$ and
$\mathbf{D}_{t}^{\nu_{2}}(\varrho_{2}U)$. Keeping in mind the homogenous
initial condition and equality (6.17), we deduce that

$\displaystyle\Gamma(1-\nu_{1})\mathbf{D}_{t}^{\nu_{1}}(\varrho_{1}U)(x,t)$
$\displaystyle=\frac{\partial}{\partial
t}\int_{0}^{t}\frac{\varrho_{1}(x,s)U(x,s)ds}{(t-s)^{\nu_{1}}}=\frac{\partial}{\partial\sigma}\int_{-\tau}^{\sigma}\frac{\varrho_{1}(x,z+\tau)U(x,z+\tau)dz}{(t-\tau-z)^{\nu_{1}}}$
$\displaystyle=\frac{\partial}{\partial\sigma}\int_{0}^{\sigma}\frac{\bar{\varrho}_{1}(x,z)\bar{U}(x,z)dz}{(\sigma-z)^{\nu_{1}}}=\Gamma(1-\nu_{1})\mathbf{D}_{\sigma}^{\nu_{1}}(\bar{\varrho}_{1}\bar{U})(x,\sigma).$

Similar calculations entail the equality

$\Gamma(1-\nu_{2})\mathbf{D}_{t}^{\nu_{2}}(\varrho_{2}U)(x,t)=\mathbf{D}_{\sigma}^{\nu_{2}}(\bar{\varrho_{2}}\bar{U})(x,\sigma).$

As a result, we can rewrite problem (6.14) in the variable $\sigma$ as

$\begin{cases}\mathbf{D}_{\sigma}^{\nu_{1}}(\bar{\varrho}_{1}\bar{U})-\mathbf{D}_{\sigma}^{\nu_{2}}(\bar{\varrho}_{2}\bar{U})-\bar{\mathcal{L}}_{1}\bar{U}-\mathcal{K}*\bar{\mathcal{L}}_{2}\bar{U}=\bar{f^{\star}}(x,\sigma)\quad\text{in}\quad\Omega_{\tau},\\\
\vskip
2.84526pt\cr\mathbf{D}_{\sigma}^{\nu_{1}}(\bar{\varrho}_{1}\bar{U})-\mathbf{D}_{\sigma}^{\nu_{2}}(\bar{\varrho}_{2}\bar{U})-\bar{\mathcal{M}}_{1}\bar{U}+\mathcal{K}_{0}*\bar{\mathcal{M}}_{2}\bar{U}=\bar{\psi^{\star}}\quad\text{on}\quad\partial\Omega_{\tau},\\\
\vskip 2.84526pt\cr\bar{U}(x,0)=0\quad\text{in}\quad\bar{\Omega}.\end{cases}$

Recasting the arguments of Step II for this problem, we immediately draw the
one-to-one classical solvability in
$\mathcal{C}^{2+\alpha,\frac{2+\alpha}{2}\nu_{1}}$ for $\sigma\in[0,\tau]$,
i.e., $t\in[0,2\tau]$. Other words, we have extended the solution $u(x,t)$
from $[0,\tau]$ to $[\tau,2\tau]$ if (6.4) holds. By the same token, we repeat
this procedure to continue the constructed solution on the intervals
$[i\tau,(i+1)\tau]$, $i=2,3,...,$ until the whole interval $[0,T]$ is
exhausted. This allows us to get the classical solution $u(x,t)$ on $[0,T]$,
satisfying the inequalities stated in Theorem 4.1. This completes the proof of
Theorem 4.1 under the additional assumption (6.1).

######  Remark 6.4.

In order to continue the solution $u$ from $[0,\tau]$ to $[\tau,T]$ in the DBC
or 3BC cases, the initial-boundary value problem (6.13) is replaced by the
initial-value problem:

$\begin{cases}\mathbf{D}_{t}^{\nu_{1}}v-\Delta
v=\tilde{\Phi}\quad\text{in}\quad\mathbb{R}_{2\tau}^{n},\\\
v(x,0)=0\quad\text{on}\quad\mathbb{R}^{n},\end{cases}$

and

$\tilde{\Phi}=\begin{cases}\mathbf{D}_{t}^{\nu_{1}}\tilde{u}-\Delta\tilde{u},\qquad\quad
x\in\mathbb{R}^{n},\quad t\in[0,\tau],\\\
[\mathbf{D}_{t}^{\nu_{1}}\tilde{u}-\Delta\tilde{u}]|_{t=\tau},\quad
x\in\mathbb{R}^{n},\quad t\in[\tau,2\tau].\end{cases}$

###  6.4. Step IV: Removing restriction (6.1)

To complete the proof of Theorem 4.1, we just need to remove the additional
assumption (6.1). Again, we shall only focus on the FDBC case. Define

$\mathcal{W}=\mathcal{W}(x)\quad\text{and}\quad\mathfrak{U}=\mathfrak{U}(x,t),$

and let $\mathfrak{U}=\mathfrak{U}(x,t)$ be the solution to the problem

$\begin{cases}\varrho_{1}\mathbf{D}_{t}^{\nu_{1}}\mathfrak{U}-\mathcal{L}_{1}\mathfrak{U}-\mathcal{W}=f(x,t)\qquad\quad(x,t)\in\Omega_{T},\\\
\varrho_{2}\mathbf{D}_{t}^{\nu_{2}}\mathfrak{U}-\mathcal{M}_{1}\mathfrak{U}-\mathcal{W}=\psi_{3}(x,t)\qquad(x,t)\in\partial\Omega_{T},\\\
\mathfrak{U}(x,0)=u_{0}(x)\qquad\qquad\qquad\quad\qquad\text{in}\quad\bar{\Omega}.\end{cases}$
(6.18)

Indeed, the assumptions of Theorem 4.1 (see H1-H5) allow us to exploit [19,
Theorem 2] and Remark 5.4, yielding the one-valued classical solvability of
(6.18), along with the inequality

$\displaystyle\|\mathfrak{U}\|_{\mathcal{C}^{2+\alpha,\frac{2+\alpha}{2}\nu_{1}}(\bar{\Omega}_{T})}+\|\mathbf{D}_{t}^{\nu_{1}}\mathfrak{U}\|_{\mathcal{C}^{1+\alpha,\frac{1+\alpha}{2}\nu_{1}}(\partial\Omega_{T})}+\|\mathbf{D}_{t}^{\nu_{2}}\mathfrak{U}\|_{\mathcal{C}^{1+\alpha,\frac{1+\alpha}{2}\nu_{1}}(\partial\Omega_{T})\cap\mathcal{C}^{\alpha,\frac{\nu_{1}\alpha}{2}}(\bar{\Omega}_{T})}$
$\displaystyle\quad+\|\mathbf{D}_{t}^{\nu_{1}}(\varrho_{1}\mathfrak{U})\|_{\mathcal{C}^{1+\alpha,\frac{1+\alpha}{2}\nu_{1}}(\partial\Omega_{T})\cap\mathcal{C}^{\alpha,\frac{\nu_{1}\alpha}{2}}(\bar{\Omega}_{T})}+\|\mathbf{D}_{t}^{\nu_{2}}(\varrho_{2}\mathfrak{U})\|_{\mathcal{C}^{1+\alpha,\frac{1+\alpha}{2}\nu_{1}}(\partial\Omega_{T})\cap\mathcal{C}^{\alpha,\frac{\nu_{1}\alpha}{2}}(\bar{\Omega}_{T})}$
$\displaystyle\leq
C\\{\|u_{0}\|_{\mathcal{C}^{2+\alpha}(\bar{\Omega})}+\|f\|_{\mathcal{C}^{\alpha,\frac{\nu_{1}\alpha}{2}}(\bar{\Omega}_{T})}+\|\psi_{3}\|_{\mathcal{C}^{1+\alpha,\frac{1+\alpha}{2}\nu_{1}}(\partial\Omega_{T})}\\}.$
(6.19)

Collecting this estimate with formula (10.34) in [24], and applying Corollary
3.1 in [23] and Remark 5.4, we obtain

$[\mathbf{D}_{t}^{\nu_{1}}(\varrho_{1}\mathfrak{U})-\mathbf{D}_{t}^{\nu_{2}}(\varrho_{2}\mathfrak{U})]|_{t=0}=\varrho_{1}(x,0)\mathbf{D}_{t}^{\nu_{1}}\mathfrak{U}|_{t=0}-\mathcal{W}(x).$
(6.20)

Then, coming to the original problem (1.1), (1.2), (1.5), we look for a
solution of the form

$u(x,t)=\mathcal{V}(x,t)+\mathfrak{U}(x,t),$

where the new unknown $\mathcal{V}=\mathcal{V}(x,t)$ solves the problem

$\begin{cases}\mathbf{D}_{t}^{\nu_{1}}(\varrho_{1}\mathcal{V})-\mathbf{D}_{t}^{\nu_{2}}(\varrho_{2}\mathcal{V})-\mathcal{L}_{1}\mathcal{V}-\mathcal{K}*\mathcal{L}_{2}\mathcal{V}=\mathfrak{F}\quad\text{in}\quad\Omega_{T},\\\
\vskip
2.84526pt\cr\mathbf{D}_{t}^{\nu_{1}}(\varrho_{1}\mathcal{V})-\mathbf{D}_{t}^{\nu_{2}}(\varrho_{2}\mathcal{V})-\mathcal{M}_{1}\mathcal{V}+\mathcal{K}_{0}*\mathcal{M}_{2}\mathcal{V}=\mathfrak{F}_{1}\quad\text{on}\quad\partial\Omega_{T},\\\
\vskip
2.84526pt\cr\mathcal{V}(x,0)=0\quad\text{in}\quad\bar{\Omega}.\end{cases}$
(6.21)

Here we set

$\displaystyle\mathfrak{F}$
$\displaystyle=f-\mathbf{D}_{t}^{\nu_{1}}(\varrho_{1}\mathfrak{U})+\mathbf{D}_{t}^{\nu_{2}}(\varrho_{2}\mathfrak{U})+\mathcal{L}_{1}\mathfrak{U}+\mathcal{K}*\mathcal{L}_{2}\mathfrak{U},$
$\displaystyle\mathfrak{F}_{1}$
$\displaystyle=\psi_{3}-\mathbf{D}_{t}^{\nu_{1}}(\varrho_{1}\mathfrak{U})+\mathbf{D}_{t}^{\nu_{2}}(\varrho_{2}\mathfrak{U})+\mathcal{M}_{1}\mathfrak{U}-\mathcal{K}_{0}*\mathcal{M}_{2}\mathfrak{U}.$

Relations (6.18)-(6.21) and Remarks 3.1 and 3.2 readily yield

$\|\mathfrak{F}\|_{\mathcal{C}^{\alpha,\alpha\nu_{1}/2}(\bar{\Omega}_{T})}+\|\mathfrak{F}_{1}\|_{\mathcal{C}^{1+\alpha,(1+\alpha)\nu_{1}/2}(\partial\Omega_{T})}\leq
C\big{[}\|u_{0}\|_{\mathcal{C}^{2+\alpha}(\bar{\Omega})}+\|f\|_{\mathcal{C}^{\alpha,\frac{\nu_{1}\alpha}{2}}(\bar{\Omega}_{T})}+\|\psi_{3}\|_{\mathcal{C}^{1+\alpha,\frac{1+\alpha}{2}\nu_{1}}(\partial\Omega_{T})}\big{]},$

and

$\mathfrak{F}(x,0)=0\quad
x\in\bar{\Omega}\qquad\text{and}\qquad\mathfrak{F}_{1}(x,0)=0\quad
x\in\partial\Omega,$

which tell us that the right-hand sides of (6.21) meet the additional
requirement (6.1). Thus, recasting the arguments of Steps I-III in the case of
problem (6.21), and taking into account the representation of $u(x,t)$ and
(6.4)-(6.20), we complete the proof of the theorem in the FDBC case, without
the restriction (6.1).

For the DBC or the 3BC cases, such a restriction is removed in a similar
manner, but replacing problem (6.18) by

$\begin{cases}\varrho_{1}\mathbf{D}_{t}^{\nu_{1}}\mathfrak{U}-\mathcal{L}_{1}\mathfrak{U}-\mathcal{W}=f(x,t)\quad(x,t)\in\Omega_{T},\\\
\mathfrak{U}(x,0)=u_{0}(x)\quad x\in\bar{\Omega},\end{cases}$

subject either to the Dirichlet or the Neumann boundary condition. The proof
of Theorem 4.1 is now finished. ∎

##  7\. Numerical Simulations

Once the well-posedness of the problem is established by our main Theorem 4.1
(see also Remark 4.4), one might like to find explicit solutions. To this aim,
we implement a numerical scheme, and we apply it to a number of cases. In the
forthcoming Examples 7.1-7.4 we consider a one dimensional domain $\Omega$,
whereas the last Example 7.5 is set in a two-dimensional domain. In
particular, we examine the case

$\mathbf{D}_{t}^{\nu_{1}}(\varrho_{1}u)-\mathbf{D}_{t}^{\nu_{2}}(\varrho_{2}u)=\frac{\partial}{\partial
t}(\mathcal{N}*u),$

where the kernel

$\mathcal{N}=\varrho_{1}(x)\omega_{1-\nu_{1}}(t)-\varrho_{2}(x)\omega_{1-\nu_{2}}(t)$

is possibly nonpositive, as in Example 7.2.

Coming to Examples 7.1-7.4, we focus on the initial-boundary value problem in
the one-dimensional domain $\Omega=(0,1)$

$\begin{cases}\varrho_{1}(x)\mathbf{D}_{t}^{\nu_{1}}u-\varrho_{2}(x,t)\mathbf{D}_{t}^{\nu_{2}}u-\mathfrak{a}(x,t)\frac{\partial^{2}u}{\partial
x^{2}}+\mathfrak{d}(x,t)\frac{\partial u}{\partial x}\
-(\mathcal{K}*b\frac{\partial^{2}u}{\partial x^{2}})=f(x,t)\quad\text{in
}\Omega_{T},\\\ \vskip 0.85358pt\cr u(x,0)=u_{0}(x),\qquad\qquad\qquad\quad
x\in[0,1],\\\ \vskip 5.69054pt\cr\mathfrak{c}_{1}\frac{\partial u}{\partial
x}(0,t)+\mathfrak{c}_{2}u(0,t)=\varphi_{1}(t),\quad t\in[0,T],\\\ \vskip
5.69054pt\cr\mathfrak{c}_{3}\frac{\partial u}{\partial
x}(1,t)+\mathfrak{c}_{4}u(1,t)=\varphi_{2}(t),\quad t\in[0,T].\end{cases}$
(7.1)

We introduce the space-time mesh with nodes

$x_{k}=kh,\quad\sigma_{j}=j\sigma,\quad k=0,1,\ldots,K,\quad
j=0,1,\ldots,J,\quad h=L/K,\quad\sigma=T/J.$

For these examples, we actually take $L=1$, $K=10^{3}$ and $J=10^{2}$.
Denoting the finite-difference approximation of the function $u$ at the point
$(x_{k},\sigma_{j})$ by $u^{j}_{k}$, and calling

$\displaystyle\mathfrak{a}^{j+1}_{k}=\mathfrak{a}(x_{k},\sigma_{j+1}),\qquad\mathfrak{d}^{j+1}_{k}=\mathfrak{d}(x_{k},\sigma_{j+1}),\qquad
b^{j}_{k}=b(x_{k},\sigma_{j}),$
$\displaystyle\mathcal{K}_{m,j}=\int_{\sigma_{m}}^{\sigma_{m+1}}{\mathcal{K}}(\sigma_{j+1}-s)ds,\qquad\rho_{m}=(-1)^{m}\binom{\nu_{1}}{m},\qquad\tilde{\rho}_{m}=(-1)^{m}\binom{\nu_{2}}{m},$
$\displaystyle\varrho_{1,k}=\varrho_{1}(x_{k}),\qquad\varrho^{j+1}_{2,k}=\varrho_{2}(x_{k},\sigma_{j+1}),$

we approximate the differential equation in (7.1) at each time level
$\sigma_{j+1}$, so to obtain the finite-difference scheme

$\displaystyle\varrho_{1,k}\sigma^{-\nu_{1}}\sum\limits_{m=0}^{j+1}(u^{j+1-m}_{k}-u_{0}(x_{k}))\rho_{m}-\varrho_{2,k}^{j+1}\sigma^{-\nu_{2}}\sum\limits_{m=0}^{j+1}(u^{j+1-m}_{k}-u_{0}(x_{k}))\tilde{\rho}_{m}$
$\displaystyle\quad-\frac{\mathfrak{a}^{j+1}_{k}}{h^{2}}(u^{j+1}_{k-1}-2u^{j+1}_{k}+u^{j+1}_{k+1})+\frac{\mathfrak{d}^{j+1}_{k}}{2h}(u^{j+1}_{k+1}-u^{j+1}_{k-1})$
$\displaystyle=\sum_{m=0}^{j}\left(b^{m}_{k}\frac{u^{m}_{k-1}-2u^{m}_{k}+u^{m}_{k+1}}{h^{2}}+b^{m+1}_{k}\frac{u^{m+1}_{k-1}-2u^{m+1}_{k}+u^{m+1}_{k+1}}{h^{2}}\right)\\!\frac{\mathcal{K}_{m,j}}{2}+f(x_{k},\sigma_{j+1}),$

for

$k=1,\ldots,K-1\qquad\text{and}\qquad j=0,1,\ldots,J-1.$

Here, the derivatives $u_{x}$ and $u_{xx}$ are approximated by the second-
order finite-difference formulas; the trapezoid-rule is employed to
approximate the integrals in the sum (see [24])

$\sum_{m=0}^{j}\int^{\sigma_{m+1}}_{\sigma_{m}}{\mathcal{K}}(\sigma_{j+1}-s)b(x,s)u_{xx}(x,s)ds;$

and the Grünwald-Letnikov formula [4] is applied to approximate the fractional
derivatives $\mathbf{D}_{t}^{\nu_{1}}u$ and $\mathbf{D}_{t}^{\nu_{2}}u$. It is
worth noting that an improvement in the accuracy of the approximation of the
fractional derivatives is achieved here by the Richardson extrapolation, see
[4]. Finally, two fictitious mesh points outside the spatial domain to
approximate the derivatives in the boundary conditions with the second order
of accuracy are exploited (see, e.g., [24]). Further improvement in the
accuracy of calculations may be reached by resorting to finite element methods
[12, 44, 45], albeit we do not have the possibility to pursue this direction
further here.

In all our examples, including in the 2-dimensional case treated later in
Example 7.5, we can exhibit the exact solution $u$, and the absolute error

$\gimel=\max|u-u_{\mathsf{N}}|$

between $u$ and the numerical solution $u_{\mathsf{N}}$, where the maximum is
taken over all the grid points in the space-time mesh, is listed in Tables
1-6.

######  Example 7.1.

Consider problem (7.1) with $T=0.1$ and

$\displaystyle\mathcal{K}(t)$
$\displaystyle=t^{-1/3},\qquad\mathfrak{a}(x,t)=\cos(\pi x/4)+t,$
$\displaystyle\mathfrak{d}(x,t)$ $\displaystyle=x+t,\qquad
b(x,t)=t^{1/3}+\sin(\pi x),$ $\displaystyle\varrho_{1}(x)$
$\displaystyle=1+x^{2},\qquad\mathfrak{c}_{1}=\mathfrak{c}_{3}=1,\qquad\mathfrak{c}_{2}=\mathfrak{c}_{4}=0,$
$\displaystyle\varphi_{1}(t)$
$\displaystyle=\varphi_{2}(t)=0,\qquad\,u_{0}(x)=\cos(\pi x),$ $\displaystyle
f(x,t)$ $\displaystyle=\pi^{2}\Big{(}\cos\frac{\pi
x}{4}+t+\frac{3t^{2/3}\sin(\pi x)}{2}+\frac{t\pi}{3\sin(\pi/3)}\Big{)}\cos(\pi
x)$ $\displaystyle\quad-(x+t)\pi\sin(\pi
x)-\frac{\varrho_{2}(x,t)t^{\nu_{1}-\nu_{2}}}{\Gamma(1+\nu_{1}-\nu_{2})}+1+x^{2}.$

As for the function $\varrho_{2}(x,t)$, we have two options:

  * (i)

$\varrho_{2}(x,t)=1+(t+1)(x+0.01)$ if $\nu_{2}=\nu_{1}/2$ with $\nu_{1}$
listed in Table 1,

  * (ii)

$\varrho_{2}(x,t)=(x-0.5)^{3}$ if $\nu_{2}=\nu_{1}/3$ with $\nu_{1}$ listed in
Table 2.

It is easy to verify that the function

$u(x,t)=\cos(\pi x)+\frac{t^{\nu_{1}}}{\Gamma(1+\nu_{1})}$

the solves initial-boundary value problem (7.1) with the parameters specified
above. The outcomes of this example (the absolute errors and the plot of
numerical and analytical solutions) are given in Figure 1, Tables 1 and 2.

Table 1. Values of $\gimel$ in Example 7.1; $\varrho_{2}(x,t)=1+(t+1)(x+0.01)$, $\nu_{2}=\nu_{1}/2$. $\nu_{1}$ | $\gimel$
---|---
$0.1$ | 1.6544e-02
$0.2$ | 4.2775e-03
$0.3$ | 2.1238e-03
$0.4$ | 1.0632e-03
$0.5$ | 5.1204e-04
$0.6$ | 2.3459e-04
$0.7$ | 9.9984e-05
$0.8$ | 3.8166e-05
$0.9$ | 2.1979e-05
Table 2. Values of $\gimel$ in Example 7.1; $\varrho_{2}(x,t)=(x-1/2)^{3}$, $\nu_{2}=\nu_{1}/3$. $\nu_{1}$ | $\gimel$
---|---
$0.1$ | 8.7910e-04
$0.2$ | 1.4009e-04
$0.3$ | 3.6891e-04
$0.4$ | 3.5521e-04
$0.5$ | 2.4190e-04
$0.6$ | 1.3600e-04
$0.7$ | 6.5636e-05
$0.8$ | 2.6783e-05
$0.9$ | 1.1683e-05
Figure 1. Exact and numerical solutions in Example 7.1 at $t=0.1$,
$\nu_{1}=0.1$, $\nu_{2}=0.05,$ $\varrho_{2}(x,t)=1+(t+1)(x+0.01)$.

######  Example 7.2.

In this test we examine (7.1) with $\varrho_{2}(x,t)=0.5$ and
$\varrho_{2}(x,t)=2.2$ for $T=0.1$ and $T=0.7$, the remaining parameters being
as in Example 7.1. The corresponding results are reported in Table 3. In
Figures 2 and 3 we plot the kernel $\mathcal{N}$ for the different choice of
parameters. Note that $\mathcal{N}$ changes its sign in the considered time
period.

(a) $T=0.1$ (b) $T=0.7$

Figure 2. Example 7.2: $\varrho_{2}(x,t)=0.5$, $\nu_{1}=0.90$, $\nu_{2}=0.45$.

(a) $T=0.1$ (b) $T=0.7$

Figure 3. Example 7.2: $\varrho_{2}(x,t)=2.2$, $\nu_{1}=0.80$, $\nu_{2}=0.40$. Table 3. Values of $\gimel$ in Example 7.2; $\nu_{2}=\nu_{1}/2$. $\nu_{1}$ | $\varrho_{2}=0.5$ | $\varrho_{2}=2.2$
---|---|---
| $\gimel$ | $\gimel$ | $\gimel$ | $\gimel$
$0.6$ | 1.7643e-04 | 8.6293e-04 | 3.5257e-04 | 4.0691e-02
$0.7$ | 8.0630e-05 | 4.9315e-04 | 1.3883e-04 | 1.1322e-02
$0.8$ | 3.1975e-05 | 2.5790e-04 | 5.7888e-05 | 4.3900e-03
$0.9$ | 1.4874e-05 | 2.1302e-04 | 3.8448e-05 | 2.0862e-03
| $T=0.1$ | $T=0.7$ | $T=0.1$ | $T=0.7$

######  Example 7.3.

Consider problem (7.1) with $T=1$ and

$\displaystyle\mathfrak{a}(x,t)$
$\displaystyle=1,\qquad\mathfrak{d}(x,t)=0,\qquad b(x,t)=1,$
$\displaystyle\varrho_{1}(x)$
$\displaystyle=1+x,\qquad\varrho_{2}(x,t)=t\sin(2\pi x),$
$\displaystyle\mathcal{K}(t)$
$\displaystyle=\frac{t^{-\nu_{1}}}{\Gamma(1-\nu_{1})},\qquad\mathfrak{c}_{1}=\mathfrak{c}_{3}=1,\qquad\mathfrak{c}_{2}=\mathfrak{c}_{4}=0,$
$\displaystyle u_{0}(x)$ $\displaystyle=\cos(\pi
x),\qquad\varphi_{1}(t)=\varphi_{2}(t)=0,$ $\displaystyle f(x,t)$
$\displaystyle=\cos(\pi
x)\bigg{[}(1+x)\Gamma(1+\nu_{1})+\pi^{2}(1+t^{\nu_{1}})+\pi^{2}t(1+\Gamma(1+\nu_{1}))+\frac{1+x+\pi^{2}}{\Gamma(2-\nu_{1})}t^{1-\nu_{1}}$
$\displaystyle\quad+\frac{\pi^{2}}{\Gamma(3-\nu_{1})}t^{2-\nu_{1}}-\left(\frac{t^{2-\nu_{2}}}{\Gamma(2-\nu_{2})}+\frac{\Gamma(1+\nu_{1})t^{1+\nu_{1}-\nu_{2}}}{\Gamma(1+\nu_{1}-\nu_{2})}\right)\sin(2\pi
x)\bigg{]}.$

Here, the analytic solution reads

$u(x,t)=[1+t+t^{\nu_{1}}]\cos(\pi x).$

The outcomes of this example are listed in Table 4.

Table 4. Values of $\gimel$ in Example 7.3; $\nu_{2}=\nu_{1}/2$. $\nu_{1}$ | $\gimel$
---|---
$0.15$ | 7.4473e-04
$0.25$ | 1.2041e-03
$0.35$ | 1.1158e-03
$0.45$ | 6.5545e-04
$0.55$ | 2.5780e-04
$0.65$ | 2.1305e-04
$0.75$ | 2.6327e-04
$0.85$ | 2.7676e-04
$0.95$ | 1.7288e-04

######  Example 7.4.

Consider problem (7.1) with $T=1$ and

$\displaystyle\mathfrak{a}(x,t)$
$\displaystyle=(x+1)(t+1),\qquad\mathfrak{d}(x,t)=x\sin t,\,b(x,t)=0,$
$\displaystyle\varrho_{2}(x,t)$ $\displaystyle=t\cos(2\pi
x),\qquad\varrho_{1}(x)=2+\sin(2\pi x),$ $\displaystyle\mathcal{K}(t)$
$\displaystyle=0,\qquad\mathfrak{c}_{1}=\mathfrak{c}_{3}=1,\qquad\mathfrak{c}_{2}=-2,\quad\mathfrak{c}_{4}=0,$
$\displaystyle u_{0}(x)$
$\displaystyle=2x-x^{2},\qquad\varphi_{1}(t)=2E_{\nu_{1}}(t^{\nu_{1}}),\qquad\varphi_{2}(t)=0,$
$\displaystyle f(x,t)$
$\displaystyle=E_{\nu_{1}}(t^{\nu_{1}})\bigl{[}(2x-x^{2})(2+\sin(2\pi
x))+2(x+1)(t+1)+x(2-2x)\sin t\bigr{]}$
$\displaystyle\quad-t^{1-\nu_{2}}\cos(2\pi
x)(2x-x^{2})\bigl{(}E_{\nu_{1},1-\nu_{2}}(t^{\nu_{1}})-1/\Gamma(1-\nu_{2})\bigr{)},$

whose exact solution is

$u(x,t)=[2x-x^{2}]E_{\nu_{1}}(t^{\nu_{1}}).$

The outcomes of this example are listed in Table 5.

Table 5. Values of $\gimel$ in Example 7.4; $\nu_{2}=\nu_{1}/2$. $\nu_{1}$ | $\gimel$
---|---
$0.1$ | 4.6741e-03
$0.2$ | 3.3408e-03
$0.3$ | 1.9065e-03
$0.4$ | 8.0956e-04
$0.5$ | 3.3009e-04
$0.6$ | 2.2661e-04
$0.7$ | 1.7038e-04
$0.8$ | 1.1417e-04
$0.9$ | 4.6430e-05
Figure 4. Exact and numerical solutions in Example 7.4 at $t=T$,
$\nu_{1}=0.9$, $\nu_{2}=0.45$.

Our last test is set in the two-dimensional domain
$\Omega=(0,L_{x})\times(0,L_{y})$. Let us briefly describe the finite-
difference scheme exploited in this case. We rewrite (1.1), (1.2), (1.4) in
the more suitable form

$\begin{cases}\varrho_{1}(x,y)\mathbf{D}_{t}^{\nu_{1}}u-\varrho_{2}(x,y,t)\mathbf{D}_{t}^{\nu_{2}}u-\mathfrak{a}^{1}(x,y,t)\frac{\partial^{2}u}{\partial
x^{2}}-\mathfrak{a}^{2}(x,y,t)\frac{\partial^{2}u}{\partial y^{2}}\\\ \vskip
2.84526pt\cr\quad+\mathfrak{d}^{1}(x,y,t)\frac{\partial u}{\partial
x}+\mathfrak{d}^{2}(x,y,t)\frac{\partial u}{\partial
y}-(\mathcal{K}*[b^{1}\frac{\partial^{2}u}{\partial
x^{2}}+b^{2}\frac{\partial^{2}u}{\partial y^{2}}])=f(x,y,t)\quad\text{in
}\Omega_{T},\\\ \vskip 5.69054pt\cr
u(x,y,0)=u_{0}(x,y),\quad(x,y)\in\bar{\Omega},\\\ \vskip 5.69054pt\cr
u(0,y,t)=u(L_{x},y,t)=0,\quad t\in[0,T],\,y\in[0,L_{y}],\\\ \vskip
5.69054pt\cr\frac{\partial u}{\partial y}(x,0,t)=\frac{\partial u}{\partial
y}(x,L_{y},t)=0,\quad t\in[0,T],\,x\in[0,L_{x}],\end{cases}$ (7.2)

and we introduce the space-time mesh with nodes

$\displaystyle x_{k}$ $\displaystyle=kh_{x},\quad
y_{l}=lh_{y},\quad\sigma_{j}=j\sigma,\quad k=0,1,\ldots,K_{x},\quad
l=0,1,\ldots,K_{y},\quad j=0,1,\ldots,J,$ $\displaystyle h_{x}$
$\displaystyle=L_{x}/K_{x},\quad h_{y}=L_{y}/K_{y},\quad\sigma=T/J.$

At each time level $\sigma_{j+1}$, we approximate the differential equation in
(7.2) via the finite-difference scheme

$\displaystyle\varrho_{1,k,l}\sigma^{-\nu_{1}}\sum\limits_{m=0}^{j+1}[u^{j+1-m}_{k,l}-u_{0}(x_{k},y_{l})]\rho_{m}-\varrho_{2,k,l}^{j+1}\sigma^{-\nu_{2}}\sum\limits_{m=0}^{j+1}[u^{j+1-m}_{k,l}-u_{0}(x_{k},y_{l})]\tilde{\rho}_{m}$
$\displaystyle\quad-\frac{\mathfrak{a}^{1}_{k,l,j+1}}{h_{x}^{2}}[u^{j+1}_{k-1,l}-2u^{j+1}_{k,l}+u^{j+1}_{k+1,l}]-\frac{\mathfrak{a}^{2}_{k,l,j+1}}{h_{y}^{2}}[u^{j+1}_{k,l-1}-2u^{j+1}_{k,l}+u^{j+1}_{k,l+1}]$
$\displaystyle\quad+\frac{\mathfrak{d}^{1}_{k,l,j+1}}{2h_{x}}(u^{j+1}_{k+1,l}-u^{j+1}_{k-1,l})+\frac{\mathfrak{d}^{2}_{k,l,j+1}}{2h_{y}}(u^{j+1}_{k,l+1}-u^{j+1}_{k,l-1})$
$\displaystyle=\sum_{m=0}^{j}\bigg{[}b^{1}_{k,l,m}\frac{u^{m}_{k-1,l}-2u^{m}_{k,l}+u^{m}_{k+1,l}}{h^{2}_{x}}+b^{1}_{k,l,m+1}\frac{u^{m+1}_{k-1,l}-2u^{m+1}_{k,l}+u^{m+1}_{k+1,l}}{h_{x}^{2}}$
$\displaystyle\quad+b^{2}_{k,l,m}\frac{u^{m}_{k,l-1}-2u^{m}_{k,l}+u^{m}_{k,l+1}}{h^{2}_{y}}+b^{2}_{k,l,m+1}\frac{u^{m+1}_{k,l-1}-2u^{m+1}_{k,l}+u^{m+1}_{k,l+1}}{h_{y}^{2}}\bigg{]}\,\frac{\mathcal{K}_{m,j}}{2}+f(x_{k},y_{l},\sigma_{j+1}),$

for

$k=1,\ldots,K_{x}-1,\qquad l=1,\ldots,K_{y}-1,\qquad j=0,1,\ldots,J-1.$

Here we called $u_{k,l}^{j}$ the finite-difference approximation of the
function $u$ at the point $(x_{k},y_{l},\sigma_{j})$, and

$\displaystyle\mathfrak{a}^{1}_{k,l,j+1}=\mathfrak{a}^{1}(x_{k},y_{l},\sigma_{j+1}),\qquad\mathfrak{a}^{2}_{k,l,j+1}=\mathfrak{a}^{2}(x_{k},y_{l},\sigma_{j+1}),\qquad\mathfrak{d}^{1}_{k,l,j+1}=\mathfrak{d}^{1}(x_{k},y_{l},\sigma_{j+1}),$
$\displaystyle\mathfrak{d}^{2}_{k,l,j+1}=\mathfrak{d}^{2}(x_{k},y_{l},\sigma_{j+1}),\qquad
b^{1}_{k,l,j}=b^{1}(x_{k},y_{l},\sigma_{j}),\qquad
b^{2}_{k,l,j}=b^{2}(x_{k},y_{l},\sigma_{j}),$
$\displaystyle\varrho_{1,k,l}=\varrho_{1}(x_{k},y_{l}),\qquad\varrho^{j+1}_{2,k,l}=\varrho_{2}(x_{k},y_{l},\sigma_{j+1}),$

while $\mathcal{K}_{m,j}$, $\rho_{m}$, $\tilde{\rho}_{m}$ are defined as in
the one-dimensional case.

######  Example 7.5.

We analyze (7.2) with $L_{x}=L_{y}=T=1$ and

$\displaystyle\mathfrak{a}^{1}$ $\displaystyle=\cos\frac{\pi
x}{4}\cos\frac{\pi y}{4}+t,\qquad\mathfrak{a}^{2}=2\cos\frac{\pi
x}{4}\cos\frac{\pi y}{4}+2t,$ $\displaystyle\mathfrak{b}^{1}$
$\displaystyle=x+y+1,\qquad\mathfrak{b}^{2}=3-x-y,\qquad\mathfrak{d}^{1}=x+y+t,\qquad\mathfrak{d}^{2}=x+y-t,$
$\displaystyle\mathcal{K}$
$\displaystyle=\frac{t^{-\nu_{1}}}{\Gamma(1-\nu_{1})},\qquad\varrho_{1}=1+x^{2}+y^{2},\qquad\varrho_{2}=1+(t+1)(x+y+0.01),\qquad
u_{0}=\sin(\pi x)\cos(\pi y),$ $\displaystyle f$
$\displaystyle=\bigg{\\{}(1+x^{2}+y^{2})\bigg{[}\Gamma(1+\nu_{1})+\frac{t^{1-\nu_{1}}}{\Gamma(2-\nu_{1})}\bigg{]}-[1+(t+1)(x+y+0.01)]\bigg{[}\frac{t^{1-\nu_{2}}}{\Gamma(2-\nu_{2})}+\frac{t^{\nu_{1}-\nu_{2}}\Gamma(1+\nu_{1})}{\Gamma(1+\nu_{1}-\nu_{2})}\bigg{]}$
$\displaystyle+3\pi^{2}(1+t+t^{\nu_{1}})\bigg{[}t+\cos\frac{\pi
x}{4}+\cos\frac{\pi
y}{4}\bigg{]}+4\pi^{2}\bigg{[}t\Gamma(1+\nu_{1})+\frac{t^{1-\nu_{1}}}{\Gamma(2-\nu_{1})}+\frac{t^{2-\nu_{1}}}{\Gamma(3-\nu_{1})}\bigg{]}\bigg{\\}}\sin(\pi
x)\cos(\pi y)$
$\displaystyle+\pi(1+t+t^{\nu_{1}})[(x+y)\cos\pi(x+y)+t\cos\pi(x-y)].$

The function

$u(x,y,t)=[1+t+t^{\nu_{1}}]\sin(\pi x)\cos(\pi y)$

solves the initial-boundary value problem (7.2) for this choice of parameters.
In our numerical calculations, we set $K_{x}=K_{y}=J=10^{2}$. Table 6 reports
the results for various values $\nu_{1}$, while Figure 5 plots the
corresponding numerical solution at $\nu_{1}=0.5$.

Table 6. Values of $\gimel$ in Example 7.5; $\nu_{2}=\nu_{1}/2$. $\nu_{1}$ | $\gimel$
---|---
$0.1$ | 6.4793e-04
$0.2$ | 7.8800e-04
$0.3$ | 5.6016e-04
$0.4$ | 3.4389e-04
$0.5$ | 3.1859e-04
$0.6$ | 3.2473e-04
$0.7$ | 3.3240e-04
$0.8$ | 3.0207e-04
$0.9$ | 1.8727e-04
Figure 5. Numerical solution in Example 7.5 at $t=T$, $\nu_{1}=0.5$,
$\nu_{2}=0.25$.

##  8\. Conclusion

In this paper, we propose an approach to study the well-posedness of initial-
boundary value problems subject to various type of boundary conditions for
multi-term fractional derivatives. Our method is particularly efficient when
the multi-term derivatives can be represented in the form
$\frac{\partial}{\partial t}(\mathcal{N}*u)$, for some nonpositive kernel
$\mathcal{N}$. We find sufficient conditions on the orders of the fractional
derivatives, providing the one-valued classical solvability in the smooth
classes. Our theoretical result are confirmed by the computational outcomes,
and the numerical examples witness the high accuracy and efficacy of the
proposed numerical schemes. A possible further development of this research
regards the inverse problem related with the identification of the parameters
in the model of oxygen subdiffusion through capillaries. Also, the complete
knowledge of the linear case is a starting point for the investigation of the
corresponding nonlinear equations, including equations with degenerate
coefficients.

## References

  * [1] P. Clément, G. Gripenberg, S.O. Londen, Schauder estimates for equations and continuous interpolation spaces, J. Differential Equa., 196 (2004) 418–447.
  * [2] V. Daftardar-Gejji, S. Bhalekar, Boundary value problems for multi-term fractional differential equations, J. Math. Anal. Appl., 345 (2008) 754–765.
  * [3] M. D’Ovidio, Fractional boundary value problems and elastic sticky Brownian motions, arXiv;2205.04162, 2022.
  * [4] K. Diethelm, N.J. Ford, A.D. Freed, Yu. Luchko, Algorithms for the fractional calculus: A selection of numerical methods, Comput. Methods Appl. Mech. Engrg., 194 (2005) 743–773; DOI: 10.1016/j.cma.2004.06.006.
  * [5] C.G. Gal, M. Warma, Elliptic and parabolic equations with fractional diffusion and dynamic boundary conditions, Evol. Equa., Control Theory, 5 (2016) 61–103.
  * [6] J.G. Go, Oxygen delivery through capillaries, Math. Biosci., 208 (2007) 166–176.
  * [7] D. Goldman, Theoretical models of microvascular oxygen transport to tissue, Microcirculation, 15 (2008) 795–811.
  * [8] D. Goldman, A.S. Popel, A computational study of the effect of vasomotion on oxygen transport from capillary networks, J. Theor. Biology, 209 (2001) 189–199.
  * [9] G.R. Goldstein, Derivation and physical interpretation of general boundary conditions, Adv. Differenti. Equa., 11 (2006) 457–480.
  * [10] J. Janno, Determination of the order of fractional derivatives and a kernel in an inverse problem for a generalized time fractional diffusion equation, Electron. J. Differntial Equa., 2016 (2016) 1–28.
  * [11] J. Janno, N. Kinash, Reconstruction of an order of derivative and a source term in a fractional diffusion equation from final measurements, Inverse Problems, 34 (2018) 02507.
  * [12] B. Jin, R. Lazarov, Z. Zhou, Numerical methods for time-fractional evolution equations with nonsmooth data: A concise overview, Comput. Methods Appl. Mech. Engrg., 346 (2019) 332–358.
  * [13] A.S. Joujchi, M.H. Derakhshan, H.R. Marasi, An efficient hybrid numerical method for multi-term time fractional partial differential equations in fluid mechanics with convergence and error analysis, Commun. Nonlinear Sci. Numer. Simulations, 114 (2022) 106620.
  * [14] J. Kempainen, K. Ruotsalainen, Boundary integral solution of the time-fractional diffusion equation, Integr. Equa. Oper. Theory, 64 (2009) 239–249.
  * [15] J. Kempainen, J. Siljander, V. Vergara, R. Zacher, Decay estimates for time-fractional and other nonlocal in time subdiffusion equations in $\mathbb{R}^{d}$, Math. Ann., 366 (2016) 941–979.
  * [16] I. Kim, H-H. Kim, S. Lim, On $L_{p}(L_{q})-$theory for the time fractional evolution equations with variable coefficients, Advances Math., 306 (2017) 123–176.
  * [17] M. Kirane, N. Tatar, Absence of local and global solutions to an elliptic system with time-fractional dynamical boundary conditions, Siberian Math. J., 48 (2007) 477–488.
  * [18] A. Kochubei, General fractional calculus, evolution equations, and renewal processes, Integr. Equ. Oper. Theory, 71 (2011) 583–600.
  * [19] M. Krasnoschok, Time-fractional diffusion equation with dynamical boundary condition, Fractional Diff. Calculus, 6 (2016) 151–178.
  * [20] M.V. Krasnoschok, Solvability in Hölder space of an initial boundary value problem for the time-fractional diffusion equation, J. Math. Phys. Anal. Geometry, 12 (2016) 48–77.
  * [21] M. Krasnoschok, V. Pata, N. Vasylyeva, Solvability of linear boundary value problems for subdiffusion equation with memory, J. Integral Equations Appl., 30 (2018) 417–445.
  * [22] M. Krasnoschok, V. Pata, N. Vasylyeva, Semilinear subdiffusion with memory in multidimensional domains, Mathematische Nachrichten, 292 (2019) 1490–1513.
  * [23] M. Krasnoschok, V. Pata, S.V. Siryk, N. Vasylyeva, Equivalent definitions of Caputo derivatives and applications to subdiffusion equations, Dynamics of PDE, 17 (2020) 383–402.
  * [24] M. Krasnoschok, S. Pereverzyev, S.V. Siryk, N. Vasylyeva, Regularized reconstruction of the order in semilinear subdiffusion with memory, (In: Cheng J., Lu S., Yamamoto M. (Eds.) Inverse Problems and Related Topics ICIP2 2018), Springer Proceedings in Mathematics$\&$Statistics, 310 (2020) 205–236, doi:10.1007/978-981-15-1592-7-10.
  * [25] M. Krasnoschok, N. Vasylyeva, On a solvability of a nonlinear fractional reaction-diffusion system in the Hölder spaces, Nonlinear Studies, 20 (2013) 589–619.
  * [26] M. Krasnoschok, N. Vasylyeva, Existence and uniqueness of the solutions for some initial-boundary value problems with the fractional dynamic boundary condition, International J. Part. Diff. Equa., 2013 (2013) ID 796430.
  * [27] A. Krogh, Oxygen diffusive shunts under conditions of heterogeneous oxygen delivery, J. Physiol., 52 (1919) 409–415.
  * [28] O.A. Ladyzhenskaia, V.A. Solonnikov, N.N. Ural’tseva, Linear and quasilinear parabolic equations, Academic Press, New York, 1968.
  * [29] R.M. Leach, D.F. Treacher, ABC of oxygen. Oxygen transport-2. Tissue hypoxia, Clinical Rev. BMJ, 317 (1998) 1370–1373.
  * [30] J. Lecoq, A. Parpaleix, E. Roussakis, M. Ducros, Y.G. Houssen, S.A. Vinogradov, et al., Simultaneous two photon imaging of oxygen and blood flow in deep cerebral vessels, Nature Medicine., 17 (2011) 893–901.
  * [31] Z. Li, Y. Liu, M. Yamamoto, Initial-boundary value problems for multi-term time-fractional diffusion equations with positive constant coefficients, Appl. Math. Comput., 257 (2015) 381–397.
  * [32] Y. Liu, M. Yamamoto, Uniqueness of orders and parameters in multi-term time-fractional diffusion equations by in exact date, arXiv:2206.02108v1, 2022.
  * [33] D. Liu, N. Wood, N. Witt, A. Hughes, S. Thom, X. Xu, Computational analysis of oxygen transport in the retinal arterial network, Current Eye Res., 34 (2009) 945–956.
  * [34] C. Lizama, G.M. Guérékata, Bounded mild solutions for semilinear integro-differential equations, J. Integral Equa. Appl., 5 (1993) 75-78.
  * [35] Y. Luchko, M. Yamamoto, General time-fractional diffusion equation: Some uniqueness and existence results for the iinitial-boundary value problems, Fract. Calc. Appl. Anal., 19 (2016) 676–695.
  * [36] Y. Luchko, A. Suzuki, M. Yamamoto, On the maximum principle for the multi-term fractional transport equation, J. Math. Anal. Appl., 505 (2022) 125579.
  * [37] V.F. Marales-Delgado, J.F. Gómez-Aguilar, K.M. Saad, M.A. Khan, P. Agarwal, Analytical solution for oxygen diffusion from capillary to tissues involving external force effects: A fractional calculus approach, Physica A, 523 (2019) 48–65.
  * [38] A.S. Popel, Oxygen diffusive shunts under conditions of heterogeneous oxygen delivery, J. Theor. Biol., 96 (1982) 533-541.
  * [39] A.S. Popel, Theory of oxygen transport to tissues, Crit. Rev. Biomed. Eng., 17 (1989) 257–321.
  * [40] T. Sandev, I.M. Sokolov, R. Metzler, A. Chechkin, Beyond monofractional kinetics, Chaos, Solitons$\&$Fractals, 102 (2017) 210–217.
  * [41] T. Sandev, R. Metzler, A. Chechkin, From continuous time random walks to the generalized diffusion equation, Fract. Calc. Appl. Anal., 21 (2018) 19–28.
  * [42] T.W. Secomb, R. Hsu, E.Y. Park, M.W. Dewhirst, Green’s function methods for analysis of oxygen delivery to tissue by microvascular networks, Annals Biomed. Engineering, 32 (2004) 1519–1529.
  * [43] C-S. Sin, Well-posedness of general Caputo-type fractional differential equations, Fract. Calc. Appl. Anal., 21 (2018) 819–832.
  * [44] S.V. Siryk, A note on the application of the Guermond-Pasquetti mass lumping correction technique for convection-diffusion problems, J. Comput. Phys., 376 (2019) 1273–1291; DOI: 10.1016/j.jcp.2018.10.016.
  * [45] S.V. Siryk, Analysis of lumped approximations in the finite-element method for convection-diffusion problems, Cybernetics and Systems Analysis, 49 No. 5, (2013) 774-784; DOI: 10.1007/s10559-013-9565-5
  * [46] V. Srivastava, K.N. Rai, A multi-term fractional diffusion equation for oxygen delivery through a capillary to tissues, Math. Comput. Modelling, 51 (2010) 616–624.
  * [47] N. Vasylyeva, L. Vynnytska, On multidimensional moving boundary problem goverened by anomalous diffusion: analytical and numerical study, NoDEA: Nonlinear Differ. Equa. Appl., 22 (2015) 543–577.
  * [48] V. Vergara, R. Zacher, Optimal decay estimates for time-fractional and other nonlocal subdiffusion equations via energy methods, SIAM J. Math. Anal., 47 (2015) 210–239.
  * [49] V. Vergara, R. Zacher, Stability, instability and blowup for time fractional and other nonlocal in time semilinear subdiffusion equations, J. Evol. Equa., 17 (2017) 599–626.
  * [50] C.Y. Wang, J.B. Bassingthwaighte, Capillary supply regions, Math. Biosci., 173 (2001) 103.
  * [51] M. Yamamoto, Fractional derivatives and time-fractional ordinary differential equations in $L_{p}-$space, arXiv:2201.07094v1, 2022.
  * [52] R. Zacher, Maximal regularity of type $L_{p}$ for abstract parabolic Voltera equations, J. Evol. Equa., 5 (2005) 79–103.
  * [53] J. Zhang, F. Liu, Z. Lin, V. Anh, Analytical and numerical solutions of a multi-term time-factional Burgers fluid model, Appl. Math. Comput., 356 (2019) 1–22.