\documentclass[reqno]{amsart} \usepackage{hyperref} \usepackage{graphicx} \AtBeginDocument{{\noindent\small \emph{Electronic Journal of Differential Equations}, Vol. 2018 (2018), No. 173, pp. 1--27.\newline ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu} \thanks{\copyright 2018 Texas State University.} \vspace{8mm}} \begin{document} \title[\hfilneg EJDE-2018/173\hfil Fractional differential equations with time dependent order] {Well-posedness, regularity, and asymptotic behavior of continuous and discrete solutions of linear fractional integro-differential equations with time-dependent order} \author[E. Cuesta, R. Ponce \hfil EJDE-2018/173\hfilneg] {Eduardo Cuesta, Rodrigo Ponce} \address{Eduardo Cuesta \newline Department of Applied Mathematics, E.T.S.I. of Telecomunication, Campus Miguel Delibes, University of Valladolid 47011, Spain} \email{eduardo@mat.uva.es} \address{Rodrigo Ponce \newline Universidad de Talca, Instituto de Matem\'atica y F\'isica, Casilla 747, Talca, Chile} \email{rponce@inst-mat.utalca.cl} \thanks{Submitted April 27, 2018. Published October 17, 2018.} \subjclass[2010]{45A05, 45E10, 45N05, 65R20, 65J08, 65J10} \keywords{Fractional integrals; Banach spaces; variable order; \hfill\break\indent convolution quadratures} \begin{abstract} We study the well-posedness of abstract time evolution fractional integro-differential equations of variable order $u(t)=u_0 + \partial^{-\alpha(t)}Au(t)+f(t)$. Also we study the asymptotic behavior as $t\to +\infty$, and the regularity of solutions. Moreover, we present the asymptotic behavior of the discrete solution provided by a numerical method based on convolution quadratures, inherited from the behavior of the continuous solution. In this equation $A$ plays the role of a linear operator of sectorial type. Several definitions proposed in the literature for the fractional integral of variable order are discussed, and the differences between the solutions provided for each of them are illustrated numerically. The definition we chose for this work is based on the Laplace transform, and we discuss the reasons for this choice. \end{abstract} \maketitle \numberwithin{equation}{section} \newtheorem{theorem}{Theorem}[section] \newtheorem{corollary}[theorem]{Corollary} \newtheorem{remark}[theorem]{Remark} \allowdisplaybreaks \section{Introduction}\label{Int} In the previous decades, fractional calculus has become a very active field of research in the framework of evolution phenomena because many of these phenomena, classically described by means of evolution equations involving integer order derivatives and/or integrals, have turned out to be better suited if non integer integrals/derivatives are introduced (see \cite{Dalir-10,Hilfer-00,Trujillo-06,Liu-16,Sheng-12,Samko-93} and references therein). The prototype linear fractional equation with constant order for evolution phenomena can be written in integral form as \begin{equation}\label{In-0} u(x,t) = u_0(x) + \int_0^t \frac{(t-s)^{\alpha-1}}{\Gamma(\alpha)} (A u)(x,s)\,\mathrm{d} s + f(x,t),\quad t>0,\; x\in\Omega, \end{equation} where $\Omega\subset\mathbb{R}^m$ is the spatial domain, $A$ is a closed linear operator acting on $u(x,t)$ in a convenient functional set, and the integral term stands for the fractional integral in the sense of Riemann-Liouville with order of integration (or viscosity parameter) $\alpha\in\mathbb{R}$, $1<\alpha<2$. More recently, fractional models involving non constant order of integration (viscosity function, $\alpha=\alpha(t)$ or $\alpha=\alpha(x,t)$) have received special attention as a natural extension of those with constant order \cite{Ak-In-Ba-17,At-Pi-11,Ch-Zh-Zh-16,Ch-Liu-Liu-Bo-16, Ji-Xu-Lin-17,Li-Li-Wu-17,Lorenzo-02,Zh-Li-An-Tu-09,Ra-Di-Ma-12,Ross-95, Samko-95,Scarpi-72}. The main reason is that practical applications modeling has demonstrated that the freedom to choose a variable order of integration/derivation instead of the constant one allows us a finer tuning in mathematical modeling. In other words, with variable order of integration more accurate modeling can be achieved. This fact applies over a large variety of applications: control theory \cite{Diaz-09,Ingman-04,Malesza-15}, mechanical engineering \cite{Coimbra-03}, image processing \cite{Cuesta-15}, remote sensing \cite{Zhou-16}, physical applications \cite{Bhrawy-17,Bhrawy-17-2}, and some others \cite{Bhrawy-15,Pedro-08,Ramirez-10,Ramirez-11}. But also theoretical properties of such kind of equations attracted the interest of researchers \cite{Almeida-09,Liu-16}, in this way it is interesting to show how these properties differ from the ones satisfied by the fractional equations with constant order. One can accept that for linear scalar equations of type \eqref{In-0}, i.e.\ for $A$ standing for a constant, some properties, as for instance the existence and uniqueness of solution, can be easily proven. However, these properties are in general not straightforwardly extended to abstract formulations of these equations as the one considered in the present work. Nonetheless there exists a extensive literature related to non local equations in time, in particular related to Volterra equations and the well-posedness, but to our knowledge, there are not general results establishing appropriate conditions on the viscosity function $\alpha(t)$ to ensure the well-posedness of that equations within a framework of general operators setting. That is why we devote this work to the study of some relevant properties of fractional equations of variable order. To be more precise, the first contribution of this work is to establish conditions, as weak as possible, for the viscosity functions $\alpha(t)$ in order to guarantee the well-posedness of linear abstract evolution equations of fractional type with time dependent order in a very general functional setting as it is the framework of complex Banach spaces. Once these conditions are stated we show the regularity exhibited by the solution as well as the asymptotic behavior as the time goes to infinity. These results are accompanied by the study of the asymptotic behavior of the numerical solution provided by backward Euler based convolution quadrature method, which is inherited from the asymptotic behavior of the continuous solution, and it is accompanied as well by several numerical experiments illustrating the theoretical results. Notice that the viscosity functions $\alpha(t)$ can be assumed to be depending on spatial variables, however in the present framework this dependence is meaningless. This paper is organized as follows. Section \ref{Back} is devoted to present a preliminary discussion on some of definitions existing in the literature for the fractional integrals with time dependent order, and to motivate our choice from all those mentioned. In Section \ref{Wp} we establish the conditions under which we can ensure the well-posedness of initial value problems of fractional type with time dependent order, and we prove the well-posedness under such assumptions. In Sections \ref{Rg} and \ref{Ab} we show the regularity of the analytic solution at time level $t=0$, and the asymptotic behavior as time goes to infinity respectively. In Section \ref{Dis} we set a time discretization based on the backward Euler convolution quadrature and we show how the asymptotic behavior as the number of steps goes to infinity is inherited from the asymptotic behavior of the analytic solution. Finally, in Section \ref{Ne} we illustrate numerically how the results of the fractional integration significantly depends on the definition we choose, and moreover we illustrate the behavior of the solutions of initial boundary value problems of fractional type with time dependent order depending on the choice of $\alpha(t)$. \section{Background on fractional integrals with time dependent order}\label{Back} The study of fractional integration with order varying in time $\alpha(t)>0$, $t\ge 0$, can be addressed in different manners depending on the definition one adopts. Assume that the order of integration/derivation ranges between $1$ and $2$, i.e.\ $1<\alpha<2$. A natural generalization of the Riemann-Liouville definition \begin{equation}\label{Back-0} \partial^{-\alpha}g(t) := \int_0^t \frac{(t-s)^{\alpha-1}}{\Gamma(\alpha)} g(s)\,\mathrm{d} s, \quad t>0, \end{equation} seems to be the one directly obtained by replacing in \eqref{Back-0} the constant order $\alpha$ by a function $\alpha(t)$. In fact, for $\alpha{:}(0,+\infty)\to (1,2)$, this definition reads \begin{equation}\label{Back-1} \partial^{-\alpha(t)}g(t) := \int_0^t \frac{(t-s)^{\alpha(t-s)-1}} {\Gamma(\alpha(t-s))} g(s)\,\mathrm{d} s,\quad t>0. \end{equation} Definition \eqref{Back-1} stands for a convolution integral $(k*g)(t)$, where the convolution kernel is given by \begin{equation}\label{Back-7} k(t)=\frac{t^{\alpha(t)-1}}{\Gamma(\alpha(t))},\quad t>0. \end{equation} However, in classical operational calculus, the numerical solution of equations involving convolution terms, e.g. the equation \eqref{Back-1}, requires the evaluation of the Laplace transform $K$ of the convolution kernel $k$, namely $K=\mathcal{L}k$, instead of $k$. In that case, the analysis of \eqref{Back-1} in terms of the Laplace transform seems to be, in general, difficult if not unaffordable, since an explicit expression of the Laplace transform of $k$, is in general not explicitly reachable. Other definitions can be found in the literature, see for example \cite{Ingman-04,Ross-95}, \begin{equation}\label{Back-2} \partial^{-\alpha(t)}g(t) := \frac{1}{\Gamma(\alpha(t))}\int_0^t (t-s)^{\alpha(s)-1} g(s)\,\mathrm{d} s,\quad t>0, \end{equation} or, with a more general formulation \cite{Malinowska-15}, \begin{equation}\label{Back-2-0} \partial^{-\alpha(\cdot,\cdot)}g(t) := \int_0^t \frac{(t-s)^{\alpha(t,s)-1}}{\Gamma(\alpha(t,s))} g(s)\,\mathrm{d} s,\quad t>0. \end{equation} Other definitions can be proposed from the definitions above. Note that the analysis of fractional integrations of variable order according definitions \eqref{Back-2} and \eqref{Back-2-0} neither looks like a straightforward matter since these equations have not convolution structure, excepting might be \eqref{Back-2-0} with a convenient choice of $\alpha(\cdot,\cdot)$. The definition of fractional integral with time dependent order we will adopt in the present work was basically proposed in \cite{Scarpi-72}, and the main feature is that it is given in terms of the Laplace transform of the convolution kernel. In fact, let $\tilde\alpha(z)$ be the Laplace transform of $\alpha(t)$, then the fractional integral of order $\alpha(t)$ is defined as the convolution integral \begin{equation}\label{Back-4} \partial^{-\alpha(t)}g(t) := \int_0^t k(t-s) g(s) \,\mathrm{d} s, \quad t>0, \end{equation} where \begin{equation}\label{Back-5} k(t) := (\mathcal{L}^{-1} K)(t),\quad t>0, \quad \text{with}\quad K(z):=\frac{1}{z^{z\tilde\alpha(z)}}, \quad z\in \mathcal{D}(K)\subset \mathbb{C}. \end{equation} We will discuss the domain of definition for $K$, $\mathcal{D}(K)$, or, equivalently, the domain of definition of the function $\tilde\alpha(z)$. Note that the definition of fractional integral of variable order according \eqref{Back-4}-\eqref{Back-5} has been already used in practical instances, in models related to image processing \cite{Cuesta-15} where only $K$ (the Laplace transform of $k$) was necessary, i.e.\ $k$ was not explicitly required. Note also that definition \eqref{Back-4}-\eqref{Back-5} is consistent with \eqref{Back-0} in the sense that, if the viscosity function $\alpha(t)$ is constant, i.e.\ $\alpha(t)=$const., then the convolution kernel turns out to be the one in \eqref{Back-0} $$ k(t) = \frac{t^{\alpha-1}}{\Gamma(\alpha)},\quad t>0. $$ Certainly the main advantage of definition \eqref{Back-4}-\eqref{Back-5} versus definitions of type \eqref{Back-1}, \eqref{Back-2}, or \eqref{Back-2-0} (the last one depending on the choice of $\alpha$), comes out when one addresses the solvability of an abstract integral equation of type \begin{equation}\label{Back-6} u(t) = u_0 + \int_0^t k(t-s) (A u)(s) \,\mathrm{d} s, \quad t>0, \end{equation} where $A$ stands for a linear operator in an abstract functional setting $X$, $A:D(A)\subset X\to X$, e.g. $A$ standing for the Laplacian operator $\Delta$ in $\mathbb{R}^n$. In the following sections, we address the solvability of fractional integral equations \eqref{Back-6} of time dependent order in the sense of definition \eqref{Back-4}-\eqref{Back-5}. Also we prove some properties of the the solution. \section{Well-posedness of fractional initial value problems} \label{Wp} Let $X$ be a complex Banach space, and let $\alpha(t)$ be a function $\alpha{:}(0,+\infty)\to (1,2)$. Consider the abstract integral equation of fractional type of variable order \begin{equation}\label{Wp-1} u(t) = u_0 + \partial^{-\alpha(t)}(Au)(t) = u_0 + \int_0^t k(t-s) (Au)(s) \,\mathrm{d} s,\quad t>0, \end{equation} where $A{:}D(A)\subset X\to X$ is a linear and closed operator of sectorial type in $X$, $u_0\in X$ is the initial data, and $k$ is a convolution kernel defined as in \eqref{Back-4}--\eqref{Back-5}, for a given complex function $\tilde\alpha(z)$ defined in certain domain $D\subset\mathbb{C}$ to be discussed below. Recall that a linear and closed operator $A$ is of \emph{sectorial type}, or $\theta$-sectorial \cite{Henry-81,Pazy-83} in $X$, if there exist $0<\theta<\pi/2$, $\omega\in\mathbb{R}$, and $M>0$ such that \begin{equation}\label{Wp-10} \|( zI - A)^{-1}\|_{X\to X}\le \frac{M}{|z-\omega|},\quad \text{for } z\notin \omega_+ + S_\theta, \end{equation} where $I$ stands for the identity operator in $X$, $\omega_+=\max\{\omega,0\}$, and \begin{equation}\label{Wp-11} \omega_+ + S_\theta := \{\omega_+ + \xi: \xi\in\mathbb{C}, |\arg (-\xi)|<\theta\}\,. \end{equation} Henceforth, for simplicity of notation, and if not confusing, for $B\in \mathcal{L}(X)$, we will write $\|B\|$ instead of $\|B\|_{X\to X}$. For the simplicity of the notation as well, hereafter we will denote $$ g(z):=z\tilde\alpha(z),\quad g_R(z):=\operatorname{Re}g(z), \quad g_I(z):=\operatorname{Imag}g(z),\quad h(z)=z^{g(z)}, $$ for $z\in D\subset \mathbb{C}$. Therefore the Laplace transform of $k$ can be written as $$ K(z)= \frac{1}{h(z)}. $$ Assume that there exist constants $D_\alpha,m_1,m_2,\varepsilon,R>0$, $0<\varepsilon^*<1$, and $0<\theta<\pi/2$, satisfying the following assumptions: \begin{itemize} \item[(A1)] The function $\alpha(t)$ admits Laplace transform $\tilde\alpha(z)$ in the complex domain $\operatorname{Re}z\ge D_\alpha$. \item[(A2)] The real part of $g(z)$, $g_R(z)$, is bounded by $$ 10$), or in finite dimensional operators such as matrices $\Delta_h\in {\mathcal M}_{M\times M}(\mathbb{R})$, in particular the ones coming out from most classical discretizations of $\Delta$. Before stating the existence and uniqueness result, we look for a convenient representation of the time evolution operator associated with the solution of \eqref{Wp-1}. In fact, since by Assumption (A1) the function $\alpha(t)$ admits Laplace transform, we can write \eqref{Wp-1} in the frequency domain as $$ U(z) = \frac{u_0}{z} + K(z)A U(z),\quad z\notin \omega_+ + S_\theta, $$ where $U(z)$ stands for the Laplace transform of the solution $u(t)$ of \eqref{Wp-1}, and $K$ is the Laplace transform of $k$. Therefore, according the notation stated above, $U(z)$ reads in terms of the resolvent of $A$ as \begin{equation}\label{Wp-8} U(z) =\frac{1}{z} (I-K(z)A)^{-1} u_0 = \frac{h(z)}{z} \,(h(z)I - A)^{-1} u_0\,. \end{equation} The well-posedness of \eqref{Wp-1} requires the existence of a bounded evolution operator $E(t)$ such that the generalized solution may be written as \begin{equation}\label{Wp-12} u(t) = E(t) u_0,\quad t>0. \end{equation} Let us show that $E(t)$ exists and, by the inversion formula of the Laplace transform, and (A1)--(A4), the evolution operator admits the expression \begin{equation}\label{Wp-9} E(t) = \frac{1}{2\pi \mathrm{i}}\int_\Gamma \mathrm{e}^{t z} \frac{h(z)}{z}\,(h(z)I - A)^{-1} \,\mathrm{d} z,\quad t>0, \end{equation} where $\Gamma$ is a convenient complex path connecting $-\infty\mathrm{i}$ and $+\infty\mathrm{i}$ with increasing imaginary part (i.e.\ positively oriented), and surrounding the sector $\omega_+ + S_\theta$. Therefore, the key point is to find one of such complex paths providing a convergent integral in \eqref{Wp-9}. In this regard define the complex paths (see Figure \ref{fig1}) \begin{gather*} \Gamma_1 : \gamma_1(\rho):= \omega_+ + a + \rho \mathrm{e}^{\pm\mathrm{i} (\pi-\theta)}, \quad 0\le \rho \le \rho_0,\\ \Gamma_2 : \gamma_2(\rho):= \rho \mathrm{e}^{\pm\mathrm{i} \varphi},\quad \rho \ge \rho_1, \end{gather*} where \begin{equation}\label{Wp-13} \rho_0=\frac{\omega_++a}{\cos(\theta)(1-\frac{\tan(\theta)}{\tan(\pi-\varphi)})}, \quad \rho_1=|z_0|=\frac{\omega_++a}{\cos(\pi-\varphi) (\frac{\tan(\pi-\varphi)}{\tan(\theta)}-1)}\,, \end{equation} $z_0$ and $\bar z_0$ are the intersection points of $\Gamma_1$ and $\Gamma_2$, i.e.\ $z_0=\omega_+ + a+ \rho_0\mathrm{e}^{\mathrm{i}(\pi-\theta)} = \rho_1\mathrm{e}^{\mathrm{i}\varphi}$, $a$ is positive constant (we will see below the role $a$ plays), and $\pm$ means that we are considering both parts of the paths, the one located in the upper half-complex plane joint with the symmetric one. The angle $\varphi$ satisfies $$ m_2\frac{\pi}{2}<\varphi\le \varepsilon^*(\pi-\theta)<\pi-\theta, $$ and without loss of generality we can set $\varphi = \varepsilon^*(\pi-\theta)$. Notice that $\varphi$ exists thanks to Assumption (A4). \begin{figure}[htb] \begin{center} \includegraphics[width=0.7\textwidth]{fig1} % path.eps \end{center} \caption{Complex paths $\Gamma_1$ and $\Gamma_2$} \label{fig1} \end{figure} Now define the continuous complex path positively oriented $\Gamma$ consisting of the union of $\Gamma_1^{1/m_2}$, and $\Gamma_2^{1/m_2}$, where $\Gamma_j^{1/m_2}$ is defined as $(\gamma_j(\rho))^{1/m_2}$, for $j=1,2$, and show that the integral \eqref{Wp-9} along $\Gamma$ is convergent. This will prove that the evolution operator $E(t)$ is well-defined. Consider first $\Gamma_1^{1/m_2}$. The following discussion will focus only on the part of $\Gamma_1^{1/m_2}$ lying in the upper half-complex plane $\operatorname{Im}z\ge0$, and the same applies for $\Gamma_2^{1/m_2}$. In the lower half-complex plane $\operatorname{Im}z\le0$ the proof follows easily. We have to prove that if $z\in\Gamma_1^{1/m_2}$, then $h(z)\notin \omega_+ +S_\theta$. To this end, we set $z\in\Gamma_1^{1/m_2}$, and denote $$ \xi=|\xi|\mathrm{e}^{\mathrm{i}\eta}=\omega_++a+\rho\mathrm{e}^{\mathrm{i}(\pi-\theta)},\quad \text{for certain } \rho>0,\text{ and } 0\le\eta\le \varphi, $$ such that $z=\xi^{1/m_2}$. Then, there is $\tilde\eta$ such that $$ h(z)=|h(z)|\mathrm{e}^{\mathrm{i} \tilde\eta}, $$ and there is also $\tau\ge 0$ such that $$ \tilde z = |\tilde z|\mathrm{e}^{\mathrm{i}\tilde \eta}=\omega_++\tau\mathrm{e}^{\mathrm{i}(\pi-\theta)}. $$ belongs to the boundary of $\omega_++S_\theta$, and has the same argument as $h(z)$. It is straightforward to prove that \begin{gather*} |\tilde z|=\frac{\omega_+\sin(\pi-\theta)}{\sin(\pi-\theta-\tilde\eta)}, \\ |h(z)|=\Big(\frac{(\omega_++a)\sin(\pi-\theta)}{\sin(\pi-\theta-\eta)} \Big)^{g_R(z)/m_2}\mathrm{e}^{-\eta g_I(z)/m_2}, \\ \tilde\eta = \frac{1}{m_2}\Big(g_I(z)\log\Big(\frac{(\omega_++a) \sin(\pi-\theta)}{\sin(\pi-\theta-\eta)}\Big) + \eta g_R(z)\Big). \end{gather*} Hence, the proof reduces to show that \begin{equation}\label{Wp-7} |h(z)|>|\tilde z|,\quad z\in \Gamma_1^{1/m_2}. \end{equation} Here we discuss several cases. First, if $\omega_+=0$, then inequality \eqref{Wp-7} obviously holds. On the other hand, if $\omega_+>0$ and $g_I=0$, then $\tilde\eta =\eta g_R(z)/m_2$ and therefore inequality \eqref{Wp-7} reduces to $$ \frac{\omega_+\sin(\pi-\theta)} {\sin\big(\pi-\theta-\eta\frac{g_R(z)}{m_2}\big)} <\Big(\frac{(\omega_++a)\sin(\pi-\theta)}{\sin(\pi-\theta-\eta)}\Big)^{g_R(z)/m_2}. $$ In that case, $g_R(z)=m_1=m_2$ leads to $\alpha(t)$ constant, and the last inequality straightforwardly holds. However, without additional assumptions, it is easy to find a naive function $g_R(z)$ (non constant, e.g. with $\eta=0$) for which this equality does not satisfy. Therefore, in the general case (i.e. with $g_I(z)$ not necessarily $0$) with $\omega_+>0$ additional assumptions are required, e.g. assumptions of type $g_I(z)/m_2$ being small enough joint with $\omega_+< 1$. Since additional assumptions for $\omega_+>0$ look like very unrealistic, now and hereafter we will assume that $\omega_+=0$, i.e. $\omega\le 0$. Consider now $\Gamma_2^{1/m_2}$, and $z\in\Gamma_2^{1/m_2}$. Therefore $z=(\rho \mathrm{e}^{\mathrm{i}\varphi})^{1/m_2}$, for certain $\rho\ge\rho_1$. In that case we have to prove that $\operatorname{Arg}(h(z))<\pi-\theta$, but this is straightforward in view of inequality $$ |\operatorname{Arg}(h(z))|\le |\log|z|\,g_I(z)| + \frac{\varphi}{m_2}g_R(z), $$ and Assumption (A3). In fact $$ |\log|z|\,g_I(z)| + \frac{\varphi}{m_2}g_R(z) < (1-\varepsilon^*)(\pi-\theta) + \varepsilon^*(\pi-\theta) = \pi-\theta, $$ and therefore $h(z)\notin S_\theta$. On the other hand, by Assumption (A2), $$ \operatorname{Arg} z =\frac{\varphi}{m_2} =\frac{\varepsilon^*(\pi-\theta)}{m_2} >\frac{\pi}{2}, $$ and therefore $\Gamma_2^{1/m_2}$ falls inside the half-complex plane $\operatorname{Re}z<0$. In conclusion, we are now in a position to prove the following theorem where, under Assumptions (A1)--(A4), we prove that the integral \eqref{Wp-9} is convergent, and equivalently Problem \eqref{Wp-1} is well-posed. \begin{theorem}\label{Wp-Th1} Let $\alpha(t)$ be a function belonging to $L^1(0,+\infty)$ satisfying (A1)--(A3). Assume also that the operator $A$ satisfies (A4). Then problem \eqref{Wp-1} is well-posed. \end{theorem} \begin{proof} The proof reduces to state that the integral in \eqref{Wp-9} is convergent for which we consider the positively oriented complex path $\Gamma$ consisting of the union of $\Gamma_1^{1/m_2}$ and $\Gamma_2^{1/m_2}$ defined above. Therefore we can write $$ E(t) = \sum_{j=1}^2 I_j,\quad\text{where}\quad I_j=\frac{1}{2\pi \mathrm{i}}\int_{\Gamma_j^{1/m_2}} \mathrm{e}^{t z} \frac{h(z)}{z}\,(h(z)I - A)^{-1} \,\mathrm{d} z, \quad j=1,2. $$ Now we prove that $I_1$ and $I_2$ are bounded. First of all, for $00$, and \eqref{Wp-12} is a solution of \eqref{Wp-1}. \end{remark} \begin{remark} \rm We have shown that the evolution operator $E(t)$, $t>0$, associated with \eqref{Wp-1} admits an integral representation \eqref{Wp-9} defined along a suitable complex path. However the path we set $\Gamma$ cannot be strictly placed in the half-complex plane $\operatorname{Re}z<0$, since the integrand does not admit holomorphic extension to the real line $\operatorname{Re}z\le 0$ ($\operatorname{Imag}z=0$). Therefore, no exponential decay will be obtained as one could expect having in mind the theory of classical $C_0$-semigroups. \end{remark} \section{Continuous solution: Regularity at $t=0$}\label{Rg} The regularity at $t=0^+$ has been already studied for $\alpha(t)=$const. in \cite{Cu-07}. In this section we extend the study to the case of $\alpha$ depending on time, i.e. $\alpha(t)$ instead of $\alpha$. Let $\delta$ be a positive non-zero constant, and define the set of functions \begin{equation}\label{Rg-1} \mathcal{F}_\delta := \big\{f{:}(0,T]\to X : \text{ measurable, } \sup_{0< t\le T} \|t^\delta f(t)\| <+\infty\big\}, \end{equation} equipped with the norm \begin{equation}\label{Rg-2} |||f|||_\delta := \sup_{0< t\le T} \|t^\delta f(t)\|,\quad f\in \mathcal{F}_\delta. \end{equation} \begin{theorem}\label{Rg-Th1} Let $u(t)$ be the solution of \eqref{Wp-1} under Assumptions {\rm (A1)--(A4)}. If $u_0\in D(A)$, then the derivative $u'(t)$, is bounded for $0< t \le T$, and $u'(t)=O(t^{m_1-1})$ as $t\to 0^+$. \end{theorem} \begin{proof} We first notice that from the resolvent identity \eqref{Wp-8} it follows that the operator $A$ commutes with the evolution operator, that is, $E(t)A\xi=AE(t)\xi$, for all $\xi\in D(A)$. Let $\Gamma$ be again the complex path defined in Section \ref{Back}. So, given $\xi\in D(A)$ we can write \begin{align*} u(t) &= E(t)\xi = \frac{1}{2\pi \mathrm{i}} \int_\Gamma \mathrm{e}^{zt} \frac{h(z)}{z}(h(z)I - A)^{-1}\xi \,\mathrm{d} z\\ &= \frac{1}{2\pi \mathrm{i}} \int_\Gamma \mathrm{e}^{zt} \Big( \frac{\xi}{z} + \frac{1}{z}(h(z)I - A)^{-1} A \xi \Big)\,\mathrm{d} z\\ &= \xi + \frac{1}{2\pi \mathrm{i}} \int_\Gamma \mathrm{e}^{zt} \frac{1}{z}(h(z)I - A)^{-1} A \xi\,\mathrm{d} z, \end{align*} for $t>0$. Note that above equality has been reached by using that $$ \frac{h(z)}{z}(h(z)I - A)^{-1}\xi = \Big(\frac{h(z)}{z}+A-A\Big)(h(z)I - A)^{-1} \xi = \xi + \frac{1}{z}\,(h(z) - A)^{-1} A \xi. $$ Therefore, denoting the derivative of $E(t)$ with respect to $t$ as $E'(t)$, we have $$ u'(t) = E'(t)\xi = \frac{1}{2\pi \mathrm{i}} \int_\Gamma \mathrm{e}^{zt} (h(z) - A)^{-1} A \xi\,\mathrm{d} z. $$ For the convenience, in the definition of $\Gamma$, we consider the parameter $a$ as a function of time; in fact $a$ takes the value $1/t^{m_2}$. So we have $$ u'(t) = \sum_{j=1}^2 I_j,\quad \text{where}\quad I_j=\frac{1}{2\pi \mathrm{i}} \int_{\Gamma_j^{1/m_2}} \mathrm{e}^{zt} (h(z)I - A)^{-1} A \xi\,\mathrm{d} z, $$ and where each individual integral will be studied separately. Note that, by the definition of $\Gamma$, we have that $h(\Gamma) \subset\mathbb{C}\setminus\{\omega_++S_\theta\}$, and therefore, for $j=1,2$, the choice of $\Gamma$ allows us to apply the sectorial property \eqref{Wp-10}--\eqref{Wp-11} of $A$. We consider two separated cases, with the same notation as in Theorem \ref{Wp-Th1}: $\omega=0$, and $\omega<0$. Hereafter in this proof we assume that $C>0$ is a generic positive constant independent of $t$. We first set $\omega=0$. Note that, by Assumption (A3) there exist $c_m,C_M>0$ such that \begin{equation}\label{Rg-2-0} c_m \le \mathrm{e}^{-\operatorname{Arg}(z)g_I(z)}\le C_M. \end{equation} If $z\in\Gamma_1^{1/m_2}$, then $z=(1/t^{m_2} + \rho\mathrm{e}^{\mathrm{i}(\pi-\theta)})^{1/m_2}$, $0\le \rho\le \rho_0$. In that case, and having in mind that $t\to 0^+$, we have \begin{equation}\label{Rg-3} |\mathrm{e}^{zt}| \le \mathrm{e}^{(1/t^{m_2})^{1/m_2} t} = \mathrm{e}, \end{equation} and there exists $C>0$ such that \begin{equation} \begin{aligned} |h(z)| &= |z|^{g_R(z)}\mathrm{e}^{-\operatorname{Arg}(z)g_I(z)} \\ &= \Big(\frac{1}{t^{m_2}} + \rho \mathrm{e}^{\mathrm{i}(\pi -\theta)}\Big)^{g_R(z)/m_2} \mathrm{e}^{-\operatorname{Arg}(z)g_I(z)} \\ &\geq \Big(\frac{1}{t^{m_2}\sin(\theta)}\Big)^{m_1/m_2}c_m \geq \frac{c_m C }{t^{m_1}}. \end{aligned} \label{Rg-4} \end{equation} On the other hand, \begin{equation}\label{Rg-5} \operatorname{length}(\Gamma_1^{1/m_2}) \le C\Big(\frac{1/t^{m_2}}{\cos(\theta) \big(1-\frac{\tan(\theta)}{\tan(\pi-\varphi)}\big)}\Big)^{1/m_2}. \end{equation} Therefore, regarding \eqref{Rg-3}--\eqref{Rg-5} we have \begin{align*} \|I_1\| &\leq \frac{M}{2\pi} \int_{\Gamma_1^{1/m_2}} \big|\frac{\mathrm{e}^{zt}}{h(z)}\big||\,\mathrm{d} z|\\ &\leq \frac{M\mathrm{e} t^{m_1}}{2\pi c_m\sin(\theta)^{m_1/m_2}} \int_{\Gamma_1^{1/m_2}} |\,\mathrm{d} z|\\ &= \frac{M\mathrm{e} t^{m_1}}{2\pi c_m\sin(\theta)^{m_1/m_2}} operatorname{length}(\Gamma_1^{1/m_2})\\ &\leq \frac{M\mathrm{e} C t^{m_1-1}}{\pi c_m\sin(\theta)^{m_1/m_2}} . \end{align*} If $z\in\Gamma_2^{1/m_2}$, then $z=(\rho\mathrm{e}^{\mathrm{i}\varphi})^{1/m_2}$, $\rho\ge \rho_1$, where $\rho_1$ was computed in \eqref{Wp-13}. In that case, there exists $C>0$ such that \begin{equation} \begin{aligned} |h(z)| &= |z|^{g_R(z)}\mathrm{e}^{-\operatorname{Arg}(z)g_I(z)} \ge c_m \rho_1^{g_R(z)/m_2} \\ &= \frac{c_m (1/t)^{g_R(z)}} {\big(\cos(\pi-\varphi)\big(\frac{\tan(\pi-\varphi)}{\tan(\theta)}-1 \big)\big)^{g_R(z)/m_2}} \\ &\ge \frac{C}{t^{m_1}}. \end{aligned}\label{Rg-7} \end{equation} On the other hand, \begin{equation}\label{Rg-6} |\,\mathrm{d} z| =\frac{\rho^{1/m_2-1}}{m_2}\,\mathrm{d}\rho. \end{equation} Therefore, by \eqref{Rg-7}, and \eqref{Rg-6} we have \begin{align*} \|I_2\| &\leq \frac{M}{2\pi} \int_{\Gamma_2^{1/m_2}}\big|\frac{\mathrm{e}^{zt}}{h(z)}\big||\,\mathrm{d} z|\\ &\leq \frac{CM}{2\pi}\int_{\rho_1}^{+\infty} \mathrm{e}^{\rho^{1/m_2}\cos(\varphi/m_2)} \frac{t^{m_1}\rho^{1/m_2-1}}{m_2}\,\mathrm{d} \rho\quad (\nu=\rho^{1/m_2})\\ &\leq \frac{CM t^{m_1}}{2\pi} \int_{\rho_1^{1/m_2}}^{+\infty} \mathrm{e}^{\nu\cos(\varphi/m_2)}\,\mathrm{d}\nu\\ &\leq \frac{CM t^{m_1}}{2\pi} \frac{\mathrm{e}^{\cos(\varphi/m_2)/t}}{|\cos(\varphi/m_2)|}\\ &\leq Ct^{m_1-1}. \end{align*} The statement for $\omega=0$ is now straightforward. If $\omega<0$, then we have $$ \|u'(t)\| \le \frac{M}{2\pi}\int_\Gamma \big|\frac{\mathrm{e}^{zt}}{h(z)-\omega} \big||\,\mathrm{d} z| \|A\xi\|, $$ In this case, the bound \begin{align*} |h(z)-\omega|\ge |h(z)|-|\omega| &= \Big| \frac{1}{t^{m_2}} + \rho\mathrm{e}^{\mathrm{i}(\pi-\theta)}\Big|^{g_R(z)/m_2} \mathrm{e}^{-\operatorname{Arg}(z)g_I(z)} -|\omega|\\ &\geq \frac{c_m\sin(\theta)^{m_1/m_2}}{t^{m_1}} - |\omega|, \end{align*} and the same process as for $\omega=0$, leads to the bound $$ \|u'(t)\|\le \frac{C t^{m_1-1}}{c_m\sin(\theta)^{m_1/m_2}- t^{m_1}|\omega|}, $$ and the proof is complete. \end{proof} As occurs in classical abstract parabolic ordinary differential equations, no regularity in time is expected for no regular initial data $u_0$. In fact, for $u_0$ merely belonging to $X$, we have next the Corollaries. \begin{corollary}\label{Rg-Co1} Let $u(t)$ be the solution of \eqref{Wp-1} satisfying {\rm(A1)--(A4)}. If $u_0$ belongs to $X$, then $u'$ belongs to $\mathcal{F}_1$. \end{corollary} \begin{proof} Since the proof follows the one of Theorem \ref{Rg-Th1}, we just sketch it. Once again, let $\Gamma$ be the path defined in Section \ref{Wp}. Then, we can write $$ u'(t) = E'(t)\xi = \frac{1}{2\pi \mathrm{i}} \int_\Gamma \mathrm{e}^{zt} h(z)(h(z) - A)^{-1}\xi \,\mathrm{d} z, \quad t>0. $$ On the other hand, sectorial property (A4) of $A$, now with $\omega=0$, and the change of variable $\nu=tz$ lead to $$ \|u'(t)\| \le \frac{M}{2\pi} \int_\Gamma |\mathrm{e}^{t z}|\,\mathrm{d} z =\frac{M}{2\pi t} \int_{\Gamma^*} |\mathrm{e}^\nu| \,\mathrm{d} \nu,\quad t>0, $$ where the complex path $\Gamma^*$ results from the change of variable $\nu=tz$. The boundedness of the last integral is straightforward. Proceeding similarly for $\omega<0$ we complete the proof. \end{proof} The proof of the next corollary follows similar steps as that of Theorem \ref{Rg-Th1}, thus we omit it. \begin{corollary}\label{Rg-Co2} Let $u(t)$ be the solution of \eqref{Wp-1} satisfying {\rm (A1)--(A4)}. If $u_0\in D(A)$, then $u''\in \mathcal{F}_{2-m_1}$. \end{corollary} \section{Continuous solution: Asymptotic behavior}\label{Ab} In this section we study the asymptotic behavior in norm of the solution of \eqref{Wp-1} as $t$ approaches $+\infty$, through the associated evolution operator $E(t)$. We remark that the study when $\alpha(t)$ is constant can be found in \cite{Cu-07}. \begin{theorem}\label{Ab-Th1} Let $E(t)$ be the evolution operator \eqref{Wp-9} associated with problem \eqref{Wp-1} satisfying {\rm (A1)--(A4)}. Then, there exists a constant $C>0$ independent of $t$ such that \begin{equation} \label{Ab-1} \|E(t)\|\leq \frac{CM}{1+|\omega| t^{m_1}},\quad\text{as } t\to+\infty. \end{equation} \end{theorem} Notice that, in view of Theorem \ref{Ab-Th1}, the asymptotic behavior of the solution $u(t)$ of \eqref{Wp-1} is independent of the initial value, i.e.\ this holds for $u_0$ merely belonging to $X$. \begin{proof} Let $E(t)$ be the evolution operator \eqref{Wp-9}, \[ E(t)=\frac{1}{2\pi \mathrm{i}}\int_{\Gamma}\mathrm{e}^{z t}\frac{h(z)}{z} \left(h(z)I-A\right)^{-1}\,\mathrm{d} z,\quad t>0, \] where $\Gamma$ is a suitable path connecting $-\mathrm{i}\infty$ and $+\mathrm{i}\infty$ positively oriented. In fact for the convenience of the proof we choose again $\Gamma$ as the union of $\Gamma_1^{1/m_2}$ and $\Gamma_2^{1/m_2}$ defined in Section \ref{Wp}, and again with $a$ depending on time as $1/t^{m_2}$ . In Section \ref{Wp} we proved that $E(t)$ is bounded, for $t>0$, and now we get a finer bound. To this end, we write $$ E(t)=\sum_{j=1}^2 I_j(t), \quad t>0,\quad\text{where}\quad I_j(t)=\frac{1}{2\pi \mathrm{i}}\int_{\Gamma_j^{1/{m_2}}}\mathrm{e}^{z t} \frac{h(z)}{z}\left(h(z)I-A\right)^{-1}\,\mathrm{d} z, $$ for $j=1,2$. Consider $t>0$ large enough, and study separately $I_1(t)$ and $I_2(t)$. Along this proof $C>0$ will denote a generic constant independent of $t$. \smallskip \noindent\textbf{First part.} Consider $I_1(t)$, and $z\in\Gamma_1^{1/m_2}$. By the definition of $\Gamma_1^{1/m_2}$, $$ z=\Big( \frac{1}{t^{m_2}} + \rho \mathrm{e}^{\mathrm{i} (\pi-\theta)}\Big)^{1/m_2},\quad \text{for some } 0\le \rho\le \rho_0, $$ where $\rho_0$ was computed in \eqref{Wp-13}. Note that $|\mathrm{e}^{zt}|\le \mathrm{e}$, for $z\in\Gamma_1^{1/m_2}$. Therefore, if $\omega=0$, then there exists $C>0$ such that parametrizing $\Gamma_1^{1/m_2}$ we have \begin{equation}\label{Ab-1-0} \begin{aligned} \|I_1(t)\| &\le \frac{CM\mathrm{e}}{2\pi}\int_{\Gamma_1^{1/m_2}} \frac{1}{|z|} |\,\mathrm{d} z| \le \frac{CM\mathrm{e}\, t^{m_2}}{\pi\sin(\theta)}\int_0^{\rho_0} \,\mathrm{d} \rho \\ &= \frac{CM\mathrm{e}}{\pi\sin(\theta)\cos(\theta) \big(1-\frac{\tan(\theta)}{\tan(\pi-\varphi)}\big)}, \end{aligned} \end{equation} and therefore, $I_1(t)$ is bounded, for $t>0$. If $\omega<0$, then $$ \|I_1(t)\|\le \frac{M}{2\pi} \int_{\Gamma_1^{1/m_2}} \frac{|h(z)|\, |\mathrm{e}^{zt}|}{|z|\,|h(z)-\omega|} |\,\mathrm{d} z|. $$ For $z\in\Gamma_1^{1/m_2}$ there exists $C>0$ such that $|z|\leq C/t$, and by \eqref{Rg-2-0}, \begin{equation}\label{Ab-2} |h(z)| = |z|^{g_R(z)} \mathrm{e}^{-\operatorname{Arg}(z)g_I(z)} \le C_M \big(\frac{C}{t}\big)^{g_R(z)} \le \frac{CC_M}{t^{m_1}}. \end{equation} Moreover, there exists $C>0$ such that \begin{equation}\label{Ab-3} |z|\ge \Big(\frac{1}{t^{m_2}}\sin(\theta)\Big)^{1/m_2}\quad \Longrightarrow\quad \frac{1}{|z|}\le Ct. \end{equation} Also, we have \begin{equation}\label{Ab-4} \frac{1}{|h(z)-\omega|}\leq \frac{1}{\sin(\theta)} \frac{t^{m_2}}{1+|\omega|t^{m_2}},\quad \text{for all } z\in\Gamma_1^{1/m_2}. \end{equation} Finally, since \begin{equation}\label{Ab-5} \int_{\Gamma_1^{1/m_2}}|\,\mathrm{d} z| = \operatorname{length}(\Gamma_1^{1/m_2}) \le C \Big(\frac{1/t^{m_2}}{\cos(\theta) \big(1-\frac{\tan(\theta)}{\tan(\pi-\varphi)}\big)}\Big)^{1/m_2} \le \frac{C}{t}, \end{equation} for some $C>0$, bounds \eqref{Ab-2}--\eqref{Ab-5} lead us to $$ \|I_1(t)\|\le \frac{CM}{|\omega|t^{m_1}},\quad \text{as } t\to +\infty, $$ and by the boundedness of $E(t)$ we have \begin{equation}\label{Ab-6} \|I_1(t)\|\le \frac{CM}{1+|\omega|t^{m_1}},\quad \text{as } t\to +\infty. \end{equation} \noindent\textbf{Second part.} Consider now $I_2(t)$, and $z\in\Gamma_2^{1/m_2}$. By the definition of $\Gamma_2^{1/m_2}$, $z=(\rho\mathrm{e}^{\mathrm{i}\varphi})^{1/m_2}$ for some $\rho\ge \rho_1$. First of all observe that the choice of $\Gamma$, and more precisely the choice of $\varphi$, implies an exponential decay of $|\mathrm{e}^{zt}|$ over $\Gamma_2^{1/m_2}$. In this case we have to split the analysis into two parts: $|z|\le R$, and $|z|\ge R$, for $R>0$ large enough. Denote $\Gamma_{2,l}^{1/m_2}=\Gamma_{2}^{1/m_2}\cap \{z\in\mathbb{C}:|z|\le R\}$, and $\Gamma_{2,u}^{1/m_2}=\Gamma_{2}^{1/m_2}\cap \{z\in\mathbb{C}:|z|\ge R\}$. Consider first $\omega=0$, then \begin{align*} \|I_2(t)\| &\leq \frac{M}{2\pi}\int_{\Gamma_2^{1/m_2}} \big|\frac{\mathrm{e}^{zt}}{z} \big||\,\mathrm{d} z|\\ & = \frac{M}{2\pi}\Big(\int_{\Gamma_{2,l}^{1/m_2}} \big|\frac{\mathrm{e}^{zt}}{z}\big| |\,\mathrm{d} z| + \int_{\Gamma_{2,u}^{1/m_2}} \big|\frac{\mathrm{e}^{zt}}{z} \big| |\,\mathrm{d} z|\Big)\\ &= \frac{M}{2\pi} (I_{2,l}(t) + I_{2,u}(t)). \end{align*} On the one hand, parametrizing $\Gamma_{2,l}^{1/m_2}$ along the interval $\rho_1\le \rho\le R^{m_2}$, we have \begin{equation} \begin{aligned} I_{2,l}(t) &= \frac{1}{m_2}\int_{\rho_1}^{R^{m_2}} \frac{\exp (\cos(\varphi/m_2)t\rho^{1/m_2}) |\frac{1}{m_2}\rho^{1/m_2-1}\mathrm{e}^{\mathrm{i}\varphi/m_2}|} {|(\rho\mathrm{e}^{\mathrm{i}\varphi})^{1/m_2}|} \,\mathrm{d} \rho \\ &= \frac{1}{m_2}\int_{\rho_1}^{R^{m_2}} \frac{ \exp (\cos(\varphi/m_2)t\rho^{1/m_2})}{\rho} \,\mathrm{d} \rho \\ &\leq \int_{t\rho_1^{1/m_2}}^{tR} \frac{\exp (\cos(\varphi/m_2)\mu)}{\mu} \, \,\mathrm{d} \mu \quad (t\rho^{1/m_2}=\mu) \\ &\leq \int_{t\rho_1^{1/m_2}}^{+\infty} \frac{\exp (\cos(\varphi/m_2)\mu)}{\mu} \,\mathrm{d}\mu \le C. \end{aligned} \label{Ab-7} \end{equation} Notice that bound $C$ in \eqref{Ab-7} does not depend on $t$ because the integral lower limit $t\rho_1^{1/m_2}$ does not. Moreover, parametrizing $\Gamma_{2,u}^{1/m_2}$ along the interval $\rho\ge R^{m_2}$ we have \begin{equation} \begin{aligned} I_{2,u}(t) &\leq \frac{1}{m_2} \int_{R^{m_2}}^{+\infty} \frac{\exp (\cos(\varphi/m_2)t\rho^{1/m_2})}{\rho}\,\mathrm{d}\rho \\ &= \int_{tR}^{+\infty} \frac{\exp (\cos(\varphi/m_2)\mu)}{\mu}\,\mathrm{d}\rho\quad (t\rho^{1/m_2}=\mu) \\ &\leq \frac{\exp (\cos(\varphi/m_2)tR)}{-\cos(\varphi/m_2)tR}. \end{aligned}\label{Ab-8} \end{equation} Observe that, for $\omega=0$, \eqref{Ab-8} shows an exponential decay, more than needed to state the estimate \eqref{Ab-1} for $\omega=0$. Therefore the bounds \eqref{Ab-7} and \eqref{Ab-8} lead to the statement of theorem. Now, we consider $\omega<0$. As in the case $\omega=0$, we can write \begin{align*} \|I_2(t)\| &\leq \frac{M}{2\pi}\Big(\int_{\Gamma_{2,l}^{1/m_2}}\frac{|h(z)| \,|\mathrm{e}^{zt}|}{|z|\,|h(z)-\omega|}|\,\mathrm{d} z| + \int_{\Gamma_{2,u}^{1/m_2}}\frac{|h(z)|\,|\mathrm{e}^{zt}|}{|z|\,|h(z)-\omega|}|\,\mathrm{d} z| \Big)\\ &= \frac{M}{2\pi} (I_{2,l}(t) + I_{2,u}(t)). \end{align*} At this point, it is straightforward to show that for $z\in\Gamma_2^{1/m_2}$, \begin{equation}\label{Ab-9} |h(z) - \omega| \ge |\omega|\sin(\theta), \end{equation} and again \begin{equation}\label{Ab-10} |h(z)| = |z|^{g_R(z)} \mathrm{e}^{-\operatorname{Arg}(z)g_I(z)} \le C_M \rho^{g_R(z)/m_2}. \end{equation} Therefore, with some abuse of the notation, we have \begin{align*} I_{2,l}(t) &\leq \frac{C}{m_2 |\omega| \sin(\theta)} \int_{\rho_1}^{R^{m_2}} \rho^{g_R(z)/m_2 - 1} \exp (\cos(\varphi/m_2)t\rho^{1/m_2}) \,\mathrm{d} \rho\\ &\leq \frac{C}{|\omega|\sin(\theta)} \int_0^{+\infty} \mu^{g_R(z)-m_2} \exp (\cos(\varphi/m_2)t\mu) \mu^{m_2-1} \,\mathrm{d} \mu\quad (\rho^{1/m_2}=\mu)\\ &\leq \frac{C}{|\omega|\sin(\theta)} \int_0^{1} \mu^{m_1-1} \exp (\cos(\varphi/m_2)t\mu) \,\mathrm{d} \mu\\ &\quad + \frac{1}{|\omega|\sin(\theta)} \int_1^{+\infty} \mu^{m_2-1} \exp (\cos(\varphi/m_2)t\mu) \,\mathrm{d} \mu\\ &\leq \frac{C}{|\omega|\sin(\theta)} \Big(\frac{1}{(\cos(\varphi/m_2))^{m_1} t^{m_1}} + \frac{1}{(\cos(\varphi/m_2))^{m_2} t^{m_2}}\Big). \end{align*} Since we are assuming that $t\gg0$, we conclude that there exists $C>0$ such that \begin{equation}\label{Ab-11} I_{2,l}(t) \le \frac{C}{|\omega|t^{m_1}}. \end{equation} Finally, admitting again some abuse of notation, the parametrization of $\Gamma_2^{1/m_2}$ and \eqref{Ab-9} leads to \begin{equation} \begin{aligned} I_{2,u}(t) &\leq \int_{R^{m_2}}^{+\infty} \exp (\cos(\varphi/m_2)t\rho^{1/m_2}) \frac{\rho^{g_R(z)/m_2}}{\rho^{1/m_2}} \frac{1}{|h(z)-\omega|} \frac{\rho^{1/m_2-1}}{m_2}\,\mathrm{d}\rho \\ &\leq \frac{1}{|\omega|\sin(\theta)m_2} \int_{R^{m_2}}^{+\infty} \exp (\cos(\varphi/m_2)t\rho^{1/m_2}) \,\mathrm{d}\rho \\ &= \frac{1}{|\omega|\sin(\theta)t^{m_2}} \int_{t R}^{+\infty} \exp (\cos(\varphi/m_2)\mu)\mu^{m_2-1} \,\mathrm{d} \mu\quad (t\rho^{1/m_2}=\mu) \\ &= \frac{\Gamma(m_2-1)}{|\omega|\sin(\theta)t^{m_2}\cos(\varphi/m_2)^{m_2}}. \end{aligned} \label{Ab-12} \end{equation} Note that a different bound, finer than \eqref{Ab-12}, could be achieved, however the asymptotic behavior for $\omega<0$ is restricted by \eqref{Ab-6}; therefore the statement of theorem is satisfied and the proof is complete. \end{proof} \begin{remark} \rm By the uniqueness of the Laplace transform, the evolution operator $E(t)$ defined in \eqref{Wp-9} satisfies the equation \begin{equation}\label{eqResolent} E(t)\xi=\xi+\int_0^t k(t-s)AE(s)\xi\,\mathrm{d} s, \end{equation} for all $\xi\in X$. This means that the family of bounded operators $\{E(t)\}_{t\geq 0}$ is a {\em resolvent family}. These families of operators were introduced by Da Prato and Ianelli in \cite[Definition 1]{Da-Ia-80}, as an extension of the notion of $C_0$-semigroups to solve integro--differential equations. For general kernels $k$ in \eqref{eqResolent} and under the $1$-regularity of $k$ (see Definition in \cite[Chapter 1, Section 3]{Pruss-93}) it was proved that $\lim_{t\to +\infty} \|E(t)\|=0$, see \cite{Li-Ve-04}. However, the results in \cite{Li-Ve-04} show that if $k$ is $1$-regular, then $\|E(t)\|\leq \frac{C}{t}$ whereas the Theorem \ref{Ab-Th1} above provides a better description of the behavior of $\|E(t)\|$ as $t\to +\infty$. \end{remark} \section{Discrete solution: Definition and asymptotic behavior}\label{Dis} The time discretization of \eqref{Wp-1} has been addressed in the literature by several means, numerical inverse Laplace transform \cite{Palencia-04}, collocation methods \cite{Brunner-04}, Adomian decomposition methods \cite{Duan-13,Duan-12}, among others. In this work we focus on the convolution quadrature based methods whose convergence and stability have been deeply studied \cite{Cuesta-06,Cuesta-03}, even within the wide framework described in this work for \eqref{Wp-1}, and not only but also in the more general context of Volterra equations \cite{Cuesta-06-2,Cuesta-07}. \subsection{Convolution quadratures}\label{Dis-S1} Let $g:(0,+\infty)\to X$ be a function belonging to $L^1((0,+\infty),X)$, $N\in\mathbb{N}$, $T>0$, $\tau=T/N$, and denote $t_n=n\tau$, for $0\le n\le N$. In this section we briefly recall how convolution quadratures formulate \cite{Lubich-88}, in fact $$ \int_0^{t_n} k(t_n-s)g(s)\,\mathrm{d} s \approx \sum_{j=0}^n k_{n-j}g_j, \quad 0\le n\le N, $$ for certain weights $k_n$ defined below, and where $g_n=g(t_n)$, for $0\le n \le N$. First of all if $k:(0,+\infty)\to \mathbb{R}$ is a convolution kernel admitting Laplace transform $K$ (e.g. the ones we are considering in this paper) the inversion formula for the Laplace transform ensures the existence of a complex path $\Gamma$ connecting $-\mathrm{i}\infty$ and $+\mathrm{i}\infty$ positively oriented, such that \begin{align*} \int_0^{t_n} k(t_n-s)g(s)\,\mathrm{d} s &= \int_0^{t_n} \Big(\frac{1}{2\pi \mathrm{i}}\int_\Gamma \mathrm{e}^{\lambda (t_n-s)} K(\lambda)\,\mathrm{d} \lambda\Big)g(s)\,\mathrm{d} s\\ &= \frac{1}{2\pi \mathrm{i}}\int_\Gamma K(\lambda) \Big(\int_0^{t_n} \mathrm{e}^{\lambda (t_n-s)}g(s)\,\mathrm{d} s \Big)\,\mathrm{d} \lambda\\ &= \frac{1}{2\pi \mathrm{i}}\int_\Gamma K(\lambda) y_\lambda(t_n)\,\mathrm{d} \lambda, \end{align*} where $y_\lambda(t)$ is the solution of the initial value problem \begin{equation}\label{Dis-1} y'(t)=\lambda y(t) + g(t),\quad t>0,\quad \text{with } y(0)=0. \end{equation} Now, we set the numerical solution $\{y_n\}_{n\ge 0}$ of \eqref{Dis-1}, provided by a multistep linear method of $m$ steps \begin{equation}\label{Dis-2} \sum_{j=0}^{m} \alpha_j y_{n+j-m} = \tau \sum_{j=0}^{m} \beta_j (\lambda y_{n+j-m}+g_{n+j-m}), \end{equation} where $y_n$ stands for the approximation to $y(t_n)$, for $n\ge 0$. These methods admit a compact formulation in terms of generating functions. In fact, if \begin{gather*} P(\xi)=\alpha_0\xi^m+\alpha_1\xi^{m-1}+\ldots+\alpha_m\xi^0,\quad Q(\xi)=\beta_0\xi^m+\beta_1\xi^{m-1}+\ldots+\beta_m\xi^0, \\ Y(\xi)=\sum_{j=0}^{+\infty}y_j\xi^j\quad\text{and}\quad G(\xi)= \sum_{j=0}^{+\infty}g_j\xi^j, \end{gather*} then the numerical method \eqref{Dis-2} can be compactly written $$ P(\xi)Y(\xi)=\tau Q(\xi)(\lambda Y(\xi)+G(\xi)), $$ or equivalently \begin{equation}\label{Dis-3} Y(\xi)=\Big(\frac{\sigma(\xi)}{\tau}-\lambda\Big)^{-1}G(\xi), \end{equation} where $\sigma(\xi)$ stands for the characteristic function corresponding to the multistep linear method, i.e.\ $\sigma(\xi)=P(\xi)/Q(\xi)$. Note that all formal series are valid for $|\xi|\le r$, $00$, in that case for $\tau>0$ small enough). The Cauchy formula along the complex path $S(\nu)=r\mathrm{e}^{\nu\mathrm{i}}$, for $-\pi<\nu\le \pi$, allows us to write \begin{equation}\label{Dis-8-1} u_n=D_n u_0, \end{equation} where \begin{equation}\label{Dis-8} D_n:=\frac{1}{2\pi \mathrm{i}}\int_S \frac{1}{(1-\xi)\xi^n} h\big(\frac{\sigma(\xi)}{\tau}\big)\Big(h\big(\frac{\sigma(\xi)}{\tau}\big)I - A\Big)^{-1} \,\mathrm{d} z. \end{equation} Now, applying these ideas to \eqref{Wp-1}, we choose as multistep linear method the backward Euler method whose characteristic function reads $\sigma(\xi)=1-\xi$. In that case, the change of variable $z=(1-\xi)/\tau$, and a convenient path deformation lead to the following expression of the discrete evolution operator $D_n$ \begin{equation}\label{Dis-6} D_n = \frac{1}{2\pi \mathrm{i}}\int_{\Gamma} \frac{h(z)}{z}\,r_n(\tau z)\,(h(z)I- A)^{-1} \,\mathrm{d} z, \end{equation} where $r_n(z)=1/(1-z)^n$, $n\ge 1$. \subsection{Asymptotic behavior} The convergence, and stability of the method \eqref{Dis-6}, and results related to the representation of the numerical solution have been already studied in \cite{Cuesta-06-2,Cuesta-07,Cuesta-03}. In this section we extend the asymptotic behavior of the discrete solution studied for $\alpha(t)$ constant \cite{Cu-07,Li-17} to $\alpha(t)$ varying in time. \begin{theorem} Let $\alpha(t){:}(0,+\infty)\to(1,2)$ a function belonging to $L^1(0,+\infty)$, satisfying {\rm (A1)--(A3)}, and let $A$ be an operator satisfying {\rm (A4)}. Assume that $\omega\leq 0$. Then there exists a constant $C>0$ independent of $n$ and $\tau$ such that the numerical solution \eqref{Dis-8-1} and \eqref{Dis-6} satisfies $$ \|u_n\|\le \frac{CM}{1+|\omega|t_n^{m_1}},\quad\text{as } n\to +\infty. $$ \end{theorem} \begin{proof} First of all, we notice that, by the representation given in \cite{Cuesta-06-2,Cuesta-03} where the backward Euler based convolution quadrature is also considered, and under Assumption (A4), the discrete evolution operator \eqref{Dis-6} admits the following representation in terms of the continuous evolution operator \eqref{Wp-9}, \begin{equation}\label{Dis-7} D_n=\int_0^{+\infty} E(s)\rho_n(s)\,\mathrm{d} s,\quad n\ge 1, \end{equation} where $\rho_n(t)$ is the measure given by \[ \rho_n(t):= \frac{\mathrm{e}^{-t/\tau}}{\tau (n-1)} \big(\frac{t}{\tau}\big)^{n-1},\quad n\ge 0\,. \] Note that $\int_0^{+\infty} \rho_n(s)\,\mathrm{d} s = 1$. Therefore, the numerical solution inherits some properties of the continuous one through this representation. In fact, since $E(t)$ is bounded as we stated in Section \ref{Wp}, the discrete evolution operator $D_n$ is bounded as well, and therefore the numerical solution is bounded independently of the regularity of $u_0$. On the other hand, consider the representation \eqref{Dis-6} of $D_n$ with the complex path $\Gamma$ defined in Section \ref{Wp}. Therefore we can write $$ D_n=\sum_{j=1}^2 I_j^n,\quad\text{where}\quad I_j^n:=\frac{1}{2\pi\mathrm{i}}\int_{\Gamma_j^{1/m_2}} \frac{h(z)}{z}\,r_n(\tau z) \,(h(z)I- A)^{-1} \,\mathrm{d} z,\quad j=1,2. $$ Now we prove that both integrals, $I_1^n$ and $I_2^n$, satisfy the statement of the Theorem. Now and hereafter we assume that $\tau$ is small enough in each instance of the proof. Since the proof follows that of Theorem \ref{Ab-Th1}, we only present the key points. \smallskip \noindent\textbf{First part.} Consider first $I_1^n$, and $z\in\Gamma_1^{1/m_2}$. Therefore $$ z=\Big( \frac{1}{t_n^{m_2}} + \rho \mathrm{e}^{\mathrm{i} (\pi-\theta)}\Big)^{1/m_2},\quad\text{for some } 0\le \rho\le \rho_0. $$ One straightforwardly has that $|1-\tau z|\ge 1-1/n$, for $z\in\Gamma_1^{1/m_2}$ and $n$ large enough, and therefore \begin{equation}\label{Dis-8-0} |r_n(\tau z)|=\frac{1}{|1-\tau z|^n} \le C, \quad \text{for } z\in\Gamma_1^{1/m_2}. \end{equation} If $\omega=0$, then there exists $C>0$ such that, as in \eqref{Ab-1-0} $$ \|I_1^n\|\le \frac{CM}{2\pi}\int_{\Gamma_1^{1/m_2}} \frac{1}{|z|} |\,\mathrm{d} z| \le \frac{CM t_n^{m_2}}{\pi\sin(\theta)}\int_0^{\rho_0} \,\mathrm{d} \rho = \frac{CM}{\pi\sin(\theta)\cos(\theta) \big(1-\frac{\tan(\theta)}{\tan(\pi-\varphi)}\big)}, $$ and the boundedness of $I_1^n$ is proven. If $\omega<0$, then applying \eqref{Ab-2}--\eqref{Ab-5} there exists $C>0$ such that $$ \|I_1^n\|\le \frac{M}{2\pi} \int_{\Gamma_1^{1/m_2}} \frac{|h(z)|\, |r_n(\tau z)|}{|z|\,|h(z)-\omega|} |\,\mathrm{d} z| \le \frac{CM}{1+|\omega|t_n^{m_1}},\quad n\ge 1,\quad (n\text{ large enough}). $$ \noindent\textbf{Second part.} Consider now $I_2^n$, and $z\in\Gamma_2^{1/m_2}$. Denote as in Theorem \ref{Ab-Th1}, $\Gamma_{2,l}^{1/m_2}$, and $\Gamma_{2,u}^{1/m_2}$. First, we observe that, by the choice of $\Gamma$, and more precisely of $\varphi$, there exists $\eta>0$ such that, for $\xi=|\xi|\mathrm{e}^{\mathrm{i}\varphi}$, it holds $1/|1-\xi|\le \mathrm{e}^{\eta\cos(\varphi)|\xi|}$, and hence we easily have \begin{equation}\label{Dis-13} |r_n(\tau z)| =\big|\frac{1}{(1-\tau z)^n}\big| \le \exp \big(\eta\cos(\varphi/m_2)\rho^{1/m_2} t_n\big), \end{equation} for $n\ge 1$ and $\rho\ge \rho_1$. At this point it is important to notice that by the choice of $\varphi$, $|\exp (\eta\cos(\varphi/m_2)\rho^{1/m_2} t_n)|$ decays exponentially over $\Gamma_2^{1/m_2}$, as $n\to +\infty$. If $\omega=0$, then \begin{align*} \|I_2^n\| &\leq \frac{M}{2\pi}\int_{\Gamma_2^{1/m_2}} \big|\frac{r_n(\tau z)}{z}\big| |\,\mathrm{d} z|\\ &= \frac{M}{2\pi}\Big(\int_{\Gamma_{2,l}^{1/m_2}} \big|\frac{r_n(\tau z)}{z}\big| |\,\mathrm{d} z| + \int_{\Gamma_{2,u}^{1/m_2}} \big|\frac{r_n(\tau z)}{z}\big| |\,\mathrm{d} z|\Big) \\ &= \frac{M}{2\pi} (I_{2,l}^n + I_{2,u}^n). \end{align*} On the one hand, parametrizing $\Gamma_{2,l}^{1/m_2}$ along the interval $\rho_1\le \rho\le R^{m_2}$, regarding \eqref{Dis-13}, and following the steps of the estimation in \eqref{Ab-7}, we obtain the boundedness of $I_{2,l}^n$, i.e.\ there exists $C>0$ such that \begin{equation}\label{Dis-14} I_{2,l}^n \le C. \end{equation} Moreover, parametrizing $\Gamma_{2,u}^{1/m_2}$ along the interval $\rho\ge R^{m_2}$ we have \begin{equation} \begin{aligned} I_{2,u}^n &\leq \int_{R^{m_2}}^{+\infty} \frac{\big|\frac{1}{m_2}\rho^{1/m_2-1}\mathrm{e}^{\mathrm{i}\varphi/m_2}\big|} {|1-\tau(\rho\mathrm{e}^{\mathrm{i}\varphi})^{1/m_2}|^n |(\rho\mathrm{e}^{\mathrm{i}\varphi})^{1/m_2}|}\,\mathrm{d}\rho \\ &\leq \frac{1}{m_2} \int_{R^{m_2}}^{+\infty} \frac{1}{\rho |\tau(\rho\mathrm{e}^{\mathrm{i}\varphi})^{1/m_2}|^n }\,\mathrm{d}\rho = \frac{1}{(\tau R)^n}, \end{aligned} \label{Dis-15} \end{equation} which is bounded if $R$ is large enough, in fact if $\tau R\ge C$ for some $C>1$. The bounds \eqref{Dis-14} and \eqref{Dis-15} lead to the statement of theorem for $\omega=0$. Now, we set $\omega<0$. As in the case $\omega=0$, we can estimate $\|I_2^n\|$ as \begin{align*} \|I_2^n\| &\leq \frac{M}{2\pi}\Big(\int_{\Gamma_{2,l}^{1/m_2}}\frac{|h(z)|\, |r_n(\tau z)|}{|z|\,|h(z)-\omega|}|\,\mathrm{d} z| + \int_{\Gamma_{2,u}^{1/m_2}}\frac{|h(z)|\,|r_n(\tau z)|}{|z|\,|h(z)-\omega|}|\,\mathrm{d} z| \Big)\\ &= \frac{M}{2\pi} (I_{2,l}^n + I_{2,u}^n). \end{align*} In view of \eqref{Ab-9} and \eqref{Ab-10} it is easy to show that there exists $C>0$, for $z\in\Gamma_2^{1/m_2}$, \begin{equation}\label{Dis-18} I_{2,l} \le \frac{C}{|\omega|t_n^{m_1}}\quad\text{as }\quad n\to+\infty. \end{equation} Finally, and admitting again some abuse of the notation, we have \begin{align*} I_{2,u}^n &\leq \int_{R^{m_2}}^{+\infty} \frac{\rho^{g_R(z)/m_2}}{\rho^{1/m_2}} \frac{1}{\big|1-\tau(\rho\mathrm{e}^{\mathrm{i} \varphi})^{1/m_2}\big|^n} \frac{1}{|h(z)-\omega|}\frac{\rho^{1/m_2-1}}{m_2}\,\mathrm{d}\rho\\ &\leq \frac{C}{|\omega|\sin(\theta)m_2} \int_{R^{m_2}}^{+\infty} \frac{1}{(\tau\rho^{1/m_2})^n} \,\mathrm{d}\rho\\ &= \frac{C}{|\omega|\sin(\theta)\tau^{m_2}} \int_{\tau R}^{+\infty} \frac{1}{\mu^{n-m_2+1}}\,\mathrm{d} \mu\quad (\mu=\tau\rho^{1/m_2})\\ &= \frac{CR^{m_2}}{|\omega|\sin(\theta)(n-m_2)(\tau R)^{n}}. \end{align*} Assuming again that $\tau R\ge C>1$, the bound \eqref{Dis-18} is valid for $I_{2,u}^n$ as well, and the proof is complete. \end{proof} \begin{remark} \rm Note that the numerical solution inherits the asymptotic behavior of the continuous one, in fact in case $\omega=0$ the numerical solution is merely bounded, and if $\omega<0$ the numerical solution decays as $1/t_n^{m_1}$, for any $u_0\in X$ (i.e. not necessarily belonging to $D(A)$), \end{remark} \section{Numerical experiments}\label{Ne} We devote this section to illustrate numerically some of theoretical aspects discussed in previous sections. \subsection{On the different definitions of fractional integrals of variable order} In this sub-section we consider the three definitions we \ discussed in Section \ref{Back} for the fractional integral equations of time dependent order, recall \begin{itemize} \item[Def.1:] $\partial^{-\alpha(t)}g(t) := \int_0^t \frac{(t-s)^{\alpha(t-s)-1}}{\Gamma(\alpha(t-s))} g(s)\,\mathrm{d} s$, $t>0$. \item[Def.2:] $ \partial^{-\alpha(t)}g(t) := \frac{1}{\Gamma(\alpha(t))}\int_0^t (t-s)^{\alpha(s)-1} g(s)\,\mathrm{d} s$, $t>0$. \item[Def.3:] $\partial^{-\alpha(t)}g(t) := \int_0^t k(t-s) g(s) \,\mathrm{d} s$, $t>0$, where $k(t) := (\mathcal{L}^{-1}K)(t)$, and $K(z):= \frac{1}{z^{z\tilde\alpha(z)}}$. \end{itemize} We also consider several choices for the viscosity function $\alpha(t)$: \begin{gather} \alpha_1(t) = 0.4\sin(t) + 1.5,\quad (\tilde\alpha(z)=0.4/(z^2+1) + 1.5/z).\label{Ne-Aa1}\\ \alpha_2(t) = 0.4\cos(t) + 1.5\quad (\tilde\alpha(z)=0.4z/(z^2+1) + 1.5/z).\label{Ne-Aa2}\\ \alpha_3(t) = 0.1\mathrm{e}^{-t} + 1\quad (\tilde\alpha(z)=0.1/(z+1) + 1/z).\label{Ne-Aa3} \end{gather} In Figures \ref{fig1}--\ref{fig4} we show the results of the fractional integration with variable order according the definitions Def.1--Def.3, with the viscosity functions \eqref{Ne-Aa1}--\eqref{Ne-Aa3}, and applied to the function $g(t)=t^2$. The numerical results in Figures \ref{fig1}--\ref{fig4} has been obtained by means of the trapezoidal rule for Def.1 and Def.2, and the backward Euler based convolution quadrature \cite{Cuesta-03} for Def.3. For all of them the final time is $T=10$, the number of time steps is $N=1000$, and therefore the time step size is $h=0.01$. Figures \ref{fig1}--\ref{fig4} are organized as follows. For each figure, the first row shows the results of integrating $g(t)$ according the three definitions Def.1 -- Def.3, and the second row shows the differences between the results reached with these definitions. Finally, Figure 2 shows the results of integrating with integer and constant order $\alpha(t)=1$, and Figures \ref{fig2}--\ref{fig4} show the result of using the viscosity functions $\alpha_1(t)$, $\alpha_2(t)$, and $\alpha_3(t)$ respectively. Regarding Figures \ref{fig2}--\ref{fig5} we must point out several facts: \begin{itemize} \item The integer integration ($\alpha(t)=1$, Figure \ref{fig1}) matches perfectly with the expected results. To be more precise, since trapezoidal rule is exact for polynomials up to degree $2$ ($g(t)=t^2$) Def. 1 and Def. 2 coincide, in other words the difference between them vanishes (first column second row, Figure \ref{fig1}). On the other hand, the differences between the results provided by Def.1 and Def.2, and the results provided by Def.3 are within the order of the backward Euler convolution quadrature (see Thm.\ 2.2 \cite{Lubich-04} with $p=1$, $\beta=3$, and $\eta=1$), i.e.\ differences are $O(1)$ (second row, first and second column Figure \ref{fig2}). \item The qualitative behavior of the integration differs depending on the definition but also on the viscosity function $\alpha(t)$, it keeps clearer if Def.2 and Def.3 are compared. In fact, observe that profile provided by Def.1 seems to be less affected by the choice of $\alpha(t)$ (first column, first row Figures \ref{fig2}--\ref{fig4}), however the profiles provided by Def.2 and Def.3 are more oscillatory for $\alpha_1(t)$, that for $\alpha_2(t)$ and $\alpha_3(t)$ (second and third column, first row Figures \ref{fig2}--\ref{fig3}). \end{itemize} \begin{figure}[htb] \begin{center} \includegraphics[width=0.9\textwidth]{fig2} % ScalarInteger.eps \end{center} \caption{Integer integration of order $1$ of $g(t)=t^2$.} \label{fig2} \end{figure} \begin{figure}[htb] \begin{center} \includegraphics[width=0.9\textwidth]{fig3} % Exp1.eps \end{center} \caption{Fractional integration of $g(t)=t^2$ with $\alpha_1(t)=0.4\sin(t)+3/2$.} \label{fig3} \end{figure} \begin{figure}[htb] \begin{center} \includegraphics[width=0.9\textwidth]{fig4} % Exp2.eps \end{center} \caption{Fractional integration of $g(t)=t^2$ with $\alpha_2(t)=0.4\cos(t)+3/2$.} \label{fig4} \end{figure} \begin{figure}[htb] \begin{center} \includegraphics[width=0.9\textwidth]{fig5} % Exp3.eps \end{center} \caption{Fractional integration of $g(t)=t^2$ with $\alpha_3(t)=0.1\mathrm{e}^{-t} + 1$.} \label{fig5} \end{figure} \subsection{Initial and boundary value problem} We devote this sub-section to illustrate numerically the behavior of the solution of \eqref{Back-4}--\eqref{Back-5} when a non scalar equation is considered. Consider the integral initial-boundary value problem \begin{equation}\label{Ne-IBVP} u(x,t) = u_0(x) +\int_0^t k(t-s)(\Delta u)(x,s) \,\mathrm{d} s, \quad 0\le t \le T,\; x\in[a,b], \end{equation} where $\Delta$ stands for the one dimensional Laplacian in an interval $[a,b]$ with homogeneous Dirichlet/Newmann boundary conditions, $u_0$ stands for the initial data, and $k$ is a given convolution kernel defined as in \eqref{Back-5}. To illustrate numerically our results we consider here the viscosity functions $\alpha_1(t)$, $\alpha_3(t)$, \eqref{Ne-Aa1} and \eqref{Ne-Aa3} respectively, and to compare to the constant order fractional equations we also consider $\alpha_0(t)=\alpha$ constant. Notice that \eqref{Ne-IBVP} fits in the abstract framework stated for \eqref{Wp-1}. Since differences in the behavior of 2-D models cannot be accurately observed in a printed version i.e.\ merely by showing some selected frames, we restrict our attention to 1-D examples. In the context of 2-D models we refer the readers e.g.\ to \cite{Cuesta-15} where image processing problems have been addressed by considering piecewise constant viscosity functions, always within the framework of \eqref{Wp-1}. First of all we discretize the operator $\Delta$ over an uniform spatial mesh of size $h>0$, $x_m=a+m h$ with $h=(b-a)/M$. To this end we set a second-order difference scheme whose formulation leads to the system of integral equations \begin{equation}\label{Ne-IBVP2} \mathbf{u}(t) = \mathbf{u}_0 +\int_0^t k(t-s)(\Delta_h \mathbf{u})(s) \,\mathrm{d} s, \quad 0\le t\le T, \end{equation} where the vector $\mathbf{u}_0$ stands for the restriction of $u_0$ to the spatial mesh, and $\Delta_h$ stands for the three--diagonal matrix corresponding to the mentioned difference scheme for the Laplacian, including homogeneous Dirichlet/Newmann discrete boundary conditions. Notice that the sectorial property of $\Delta$ is inherited by $\Delta_h$, therefore the semi-discrete problem \eqref{Ne-IBVP2} fits in sectorial framework of \eqref{Wp-1} as well. The time discretization is carried out by means of the backward Euler based convolution quadrature method described in Section \ref{Dis-S1}. The fully discrete problem now reads \begin{equation}\label{Ne-IBVP3} \mathbf{u}_n = \mathbf{u}_0 + \sum_{j=1}^n q_{n-j} \Delta_h \mathbf{u}_{j}, \quad 1\le n\le N, \end{equation} where the vector $\mathbf{u}_n$ of size $(M+1)^2\times 1$ represents the approximation to the analytical solution $u(\cdot,t_n)$ restricted to the spatial mesh, with $t_n=n\tau$, for $0\le n \le N$ and the step size $\tau=T/N >0$. Let us highlight that each convolution weight $q_j$, for $j\ge0$, is the j-$th$ coefficient of the expansion provided when evaluating the Laplace transform of $k$ in $\sigma(\xi)/\tau$ according \eqref{Dis-4}, and all of them have been computed by means of Fast Fourier Transform techniques \cite{Lubich-88-2}. Since the issues related to the convergence and stability have been precisely studied in \cite{Cuesta-06-2,Cuesta-07}, in this section we will focus on showing the different behaviors of the solutions depending on the choice of $\alpha(t)$. In fact in Figure \ref{fig5} we consider three initial data (first column) on the spatial intervals $[0,1], [0,\pi]$ and $[0,\sqrt[]{\pi}]$ respectively, and homogeneous Dirichlet boundary conditions. The time discretization is carried out with the final time $T=5$, and $N=500$, and the spatial mesh is $x_m = mh$ with $h=(b-a)/M$ and $M=200$, on the respective intervals, and for each initial data. The results shown in columns 2--4 of Figure \ref{fig5} stand for the evolution of $u(x,t)$ at time levels $n=100,200$, and $500$, for each initial data. All results are compared with the results achieved with constant order of integration, in fact we have considered $\alpha_0(t)=\alpha_0$ with $\alpha_0=1.9$. In Figure \ref{fig6} we repeat the first experiment of Figures \ref{fig5}, but changing the boundary conditions now with homogeneous Newmann boundary conditions. In view of experiments in Figures \ref{fig5} and \ref{fig6}, we highlight the following: \begin{itemize} \item In all cases the solutions decay, as $t$ tends to $+\infty$, even for $\alpha_3(t)$ (last row, Figures \ref{fig5}) whose decay turns out to be slower and cannot be clearly observed for $T=5$. However in this case, for longer time e.g. $T=20$, it can be numerically observed that the solution decays to $0$ as well. \item The behavior of solutions strongly depend on the viscosity functions $\alpha(t)$, in fact the oscillatory behavior, or the decay as $t\to+\infty$. \end{itemize} \begin{figure}[htb] \begin{center} \includegraphics[width=0.95\textwidth]{fig6a} % Exp4.eps \includegraphics[width=0.95\textwidth]{fig6b} % Exp5.eps \includegraphics[width=0.95\textwidth]{fig6c} % Exp6.eps \end{center} \caption{$\alpha_0(t)=1.9$, $\alpha_1(t)=0.4\sin(t) + 1.5$, and $\alpha_3(t) = 0.9\mathrm{e}^{-t} + 1$. Boundary conditions: Dirichlet homogeneous.} \label{fig6} \end{figure} \begin{figure}[htb] \begin{center} \includegraphics[width=0.95\textwidth]{fig7} % Exp7.eps \end{center} \caption{$\alpha_0(t)=1.9$, $\alpha_1(t)=0.4\sin(t) + 1.5$, and $\alpha_3(t)=0.9\mathrm{e}^{-t} + 1$. Boundary conditions: Newmann homogeneous.}\label{fig7} \end{figure} \subsection*{Conclusions} In this paper we state conditions rather undemanding for the well-posedness of fractional integral equations of variable order in time, with special regard to the conditions for viscosity functions beyond the conditions stated in the literature for general Volterra equations. The numerical results confirmed that different choices of the viscosity function and/or the choice of the definition of fractional integral lead to solutions with very different profiles. Finally, if a source term $f(t)$ is introduced in the equation \eqref{Wp-1}, and since most parts of the proofs are carried out in terms of the continuous and discrete evolution operators, then no additional and relevant difficulties are expected to extend most of results to the non homogeneous problem. \subsection*{Acknowledgments} E. Cuesta was supported by Spanish Ministerio de Econom\'ia y Competitividad under the Research Grant MTM2014-54710-P. R. Ponce was supported by Fondecyt Grant 11130619. \begin{thebibliography}{99} \bibitem{Ak-In-Ba-17} A.~Akgul, M.~Inc, D.~Baleanu; On solutions of variable-order fractional differential equations. \emph{Int. J. Optim. Control. Theor. Appl.}, \textbf{7} (2017), 112--116. \bibitem{Almeida-09} A.~Almeida, S.~Samko; Fractional and hypersingular operators in variable exponent spaces on metrics measure spaces. \emph{Mediterr. J. Math.}, \textbf{6} (2009), 215--232. \bibitem{At-Pi-11} T.~Atanackovic, S.~Pilipovic; Hamilton's principle with variable order fractional derivatives, \emph{Fract. Cal. Appl. Anal.} \textbf{14} (2011), 94--109. \bibitem{Bhrawy-15} A.~H. Bhrawy, M.~A. Zaky; Numerical simulation for two-dimensional variable-order fractional nonlinear cable equation. \emph{Nonlinear Dyn.}, \textbf{80} (2015), 101--116. \bibitem{Bhrawy-17} A.~H. Bhrawy, M.~A. Zaky; Highly accurate numerical schemes for multi-dimensional space variable-order fractional {S}chr\"odinger equations, \emph{Comp. Math. Appl.}, \textbf{73} (2017), 1100--1117. \bibitem{Bhrawy-17-2} A.~H. Bhrawy, M.~A. Zaky; An improved collocation method for multi-dimensional space-time variable-order fractional {S}chr\"odinger equation. \emph{Appl. Numer. Math.} \textbf{111} (2017), 197--218. \bibitem{Brunner-04} H.~Brunner; \emph{ Collocation {M}ethods for {V}olterra {I}ntegral and {R}elated {F}unctional {E}quations,} vol.~15, Cambridge University Press, New York. P. G. Ciarlet, A. Iserles, R. V. Kohn, M. H. Wright (eds.), 2004. \bibitem{Cuesta-06-2} M.~P. Calvo, E.~Cuesta, C.~Palencia; \emph{Backward {E}uler method as a positivity preserving method for abstract integral equations of convolution type.} Proceedings in 2nd IFAC Workshop on Fractional Differentiation and its Applications--FDA '06, 2006, pp.~530--534. \bibitem{Cuesta-07} M.~P. Calvo, E.~Cuesta, C.~Palencia; Runge-{K}utta convolution quadrature methods for well-posed equations with memory. \emph{ Numer. Math.}, \textbf{107} (2007), no.~4, 589--614. \bibitem{Ch-Zh-Zh-16} W.~Chen, J.~Zhang, J.~Zhang; A variable-order time-fractional derivative model for chloride ions sub-diffusion in concrete structures. \emph{Fract. Cal. Appl. Anal.}, \textbf{16} (2013), 76--92. \bibitem{Ch-Liu-Liu-Bo-16} Y.~Chen, L.~Liu, D.~Liu, D.~Boutat; Numerical study of a class of variable order nonlinear fractional differential equation in terms of {B}ernstein polynomials. \emph{Eng. Physics and Math.} (2016), (in Press). \bibitem{Coimbra-03} C.~F.~M. Coimbra; Mechanics with variable--order differential operators. \emph{Ann. Phys.}, \textbf{12} (2003), no.~11-12, 692--703. \bibitem{Cu-07} E.~Cuesta; Asymptotic behaviour of the solutions of fractional integro--differential equations and some time discretizations. \emph{ Discrete Contin. Dyn. Syst., Proceedings of the 6th AIMS International Conference, suppl.}, (2007), 277--285. \bibitem{Cuesta-15} E.~Cuesta, A.~Dur\'an, M.~Kirane; \emph{On evolutionary integral models for image restoration.} Developments in Medical Image Processing and Computational Vision, Springer International Publishing, 2015, pp.~241--260. \bibitem{Cuesta-06} E.~Cuesta, Ch. Lubich, C.~Palencia; Convolution quadrature time discretizations of fractional diffusion-wave equation. \emph{Math. Comput.}, \textbf{75} (2006), no.~254, 673--696. \bibitem{Cuesta-03} E.~Cuesta, C.~Palencia; A numerical method for an integro-differential equation with memory in {B}anach spaces: {Q}ualitative properties. \emph{SIAM J. Numer. Anal.}, \textbf{41} (2003), no.~4, 1232--1241. \bibitem{Dalir-10} M.~Dalir; Applications of fractional calculus. \emph{Appl. Math. Sci.}, \textbf{4} (2010), no.~21, 1021--1032. \bibitem{Diaz-09} G.~Diaz, C.~F.~M. Coimbra; Nonlinear dynamics and control of a variable order oscillator with applications to the {V}an der {P}ol equation. \emph{Nonlinear Dyn.}, \textbf{56} (2009), no.~1-2, 145--157. \bibitem{Duan-13} J.-S. Duan, T.~Chaolu, R.~Rach, L.~Lu; The adomian decomposition method with convergence acceleration techniques for nonlinear fractional differential equations. \emph{Comp. Math. App.}, \textbf{66} (2013), 728--736. \bibitem{Duan-12} J.~S. Duan, R.~Rach, D.~Baleanu, A.~M. Wazwaz; A review of the {A}domian decomposition method and its applications to fractional differential equations. \emph{Commun. Frac. Calc.}, \textbf{3} (2012), 73--99. \bibitem{Henry-81} D.~Henry; \emph{Geometric {T}heory of {S}emilinear {P}arabolic {E}quations.}, LNM 840, Springer, Berlin, 1981. \bibitem{Hilfer-00} R.~Hilfer (ed.); \emph{Applications of {F}ractional {C}alculus in {P}hysics.} World Scientific Publishing Co. Pte. Ltd., 2000. \bibitem{Ingman-04} D.~Ingman, J.~Suzdalnitsky; Control of damping oscillations by fractional differential operator with time-dependent order. \emph{Comput. Methods Appl. Mech. Eng.}, \textbf{193} (2004), 5585--5595. \bibitem{Ji-Xu-Lin-17} Y.~Jia, M.~Xu, Y.~Lin; A numerical solution for variable order fractional functional differential equation. \emph{Applied Math. Letters}, \textbf{64} (2017), 125--130. \bibitem{Trujillo-06} A.~A. Kilbas, H.~M. Srivastava, J.~J. Trujillo; \emph{Theory and {A}pplications of {F}ractional {D}ifferential {E}quations.} North-{H}oland Mathematics Studies 204, Elsevier B. V., 2006. \bibitem{Li-Li-Wu-17} X.~Li, H.~Li, B.~Wu; A new numerical method for variable order fractional functional differential equations. \emph{Applied Math. Letters}, \textbf{68} (2017), 80--86. \bibitem{Liu-16} Z.~Liu, S.-Da Zeng; Maximum principles for multi--term space-time variable-order fractional diffusion equations and their applications. \emph{Fract. Calc. Appl. Anal.}, \textbf{19} (2016), no.~1, 188--211. \bibitem{Li-17} C.~Lizama; The Poisson distribution, abstract fractional difference equations, and stability. \emph{Proc. Amer. Math. Soc.}, \textbf{145} (2017), no.~9, 3809--3827. \bibitem{Li-Ve-04} C.~Lizama, V.~Vergara; Uniform stability of resolvent families. \emph{Proc. Amer. Math. Soc.}, \textbf{132} (2004), no.~1, 175--181. \bibitem{Palencia-04} M.~L\'opez-Fern\'andez, C.~Palencia; On the numerical inversion of the {L}aplace transform of certain holomorphic mappings. \emph{App. Numer. Math.}, \textbf{51} (2005), 289--303. \bibitem{Lorenzo-02} C.~F. Lorenzo, T.~T. Hartley; Variable order and distributed order fractional operators. \emph{Nonlinear Dyn.}, \textbf{29} (2002), no.~1-4, 57--98. \bibitem{Lubich-88} Ch. Lubich; Convolution quadrature and discretized operational calculus, {I}. \emph{Numer Math.}, \textbf{52} (1988), 129--145. \bibitem{Lubich-88-2} Ch. Lubich; Convolution quadrature and discretized operational calculus, {II}. \emph{Numer. Math.} \textbf{52} (1988), 413--425. \bibitem{Lubich-04} Ch. Lubich; Convolution quadrature revisited. \emph{BIT Numer. Math.} \textbf{44} (2004), 503--514. \bibitem{Sheng-12} R.~L. Magin (ed.); \emph{Fractional {P}rocesses and {F}ractional--{O}rder Signal Processing Techniques and A}pplications. Springer--Verlag London, 2012. \bibitem{Malesza-15} W.~Malesza, D.~Sierociuk, M.~Macias; Solution of fractional variable order differential equation. \emph{2015 American Control Conference Proceedings} (2015). \bibitem{Malinowska-15} A.~B. Malinowska, D.~F.~M. Torres, T.~Odzijewicz; \emph{Advanced {M}ethods in the {F}ractional {C}alculus of {V}ariations.} Springerbriefs in Applied Sciences and Technology, Springer International Publishing AG Switzerland, 2015. \bibitem{Zh-Li-An-Tu-09} V.~Anh P.~Zhuang, F.~Liu, I.~Turner; Numerical methods for the for the variable-order fractional advection-diffusion equation with a nonlinear source term. \emph{SIAM J. Numer. Anal.} \textbf{47} (2009), 1760--1781. \bibitem{Pazy-83} A.~Pazy; \emph{Semigroups of {L}inear {O}perators and {A}pplications to {P}artial {D}ifferential {E}quations.} Springer-Verlag, New York, 1983. \bibitem{Pedro-08} H.~T.~C. Pedro, M.~H. Kobayashi, J.~M.~C. Pereira, C.~F.~M. Coimbra; Variable order modeling of diffusive--convectie effects on the oscillatory flow past a sphere. \emph{J. Vib. Control} \textbf{14} (2008), no.~9-10, 1569--1672. \bibitem{Da-Ia-80} G.~Da Prato, M.~Iannelli; Linear integrodifferential equations in {B}anach spaces. \emph{Rend. Sem. Mat. Univ. Padova} \textbf{62} (1980), 207--219. \bibitem{Pruss-93} J.~Pr\"uss; \emph{Evolutionary {I}ntegral {E}quations and {A}pplications.} Monographs in Mathematics, vol.~87, Birkh{\"a}user Verlag, 1993. \bibitem{Ramirez-10} L.~E.~S. Ramirez, C.~F.~M. Coimbra; On the selection and meaning of variable order operators for dynamics modeling. \emph{Int. J. Differ. Equ.} (2010), 16 pages. \bibitem{Ramirez-11} L.~E.~S. Ramirez, C.~F.~M. Coimbra; On the variable order dynamics of the nonlinear wake caused by a sediment particle. \emph{Phys. D.} \textbf{240} (2011), no.~13, 1111--1118. \bibitem{Ra-Di-Ma-12} A.~Raziminia, A.~F. Dizaji, V.~J. Majd; Solution existence for non-autonomous variable-order fractional differential equations. \emph{Math. and Comput. Model.}, \textbf{55} (2012), 1106--1117. \bibitem{Ross-95} B.~Ross; Fractional integration operator of variable order in the {H}\"older spaces {H}$^{\lambda(x)}$. \emph{Internat. J. Math. \& Math. Sci.}, \textbf{18} (1995), no.~4, 777--788. \bibitem{Samko-95} S.~G. Samko; Fractional integration and differentiation with variable order. \emph{Anal. Mathematica} \textbf{21} (1995), 213--236. \bibitem{Samko-93} S.~G. Samko, A.~A. Kilbas, O.~I. Marichev; \emph{Fractional {I}ntegrals and {D}erivatives.} Gordon and Breach Science Publishers, 1993. \bibitem{Scarpi-72} G.~Scarpi; Sulla posibilit\'a di un modello reologico intermedio di tipo evolutivo. \emph{Rend. Sc. Fis. Mat. e Nat.} \textbf{LII} (1972), 570--575. \bibitem{Zhou-16} Q.~Zhou, J.~Gao, Z.~Wang, Kexue Li; Adaptive variable time fractional anisotropic diffusion filtering for seismic data noise attenuation. \emph{IEEE Trans. Geosci. Remote Sens.} \textbf{54} (2016), no.~4, 1905--1917. \end{thebibliography} \end{document}