# Time-fractional porous medium equation: Erdélyi-Kober integral equations,
compactly supported solutions, and numerical methods

Belen López, Hanna Okrasińska-Płociniczak,
Łukasz Płociniczak, Juan Rocha∗ Departamento de Matemáticas, Universidad de
Las Palmas de Gran Canaria, Campus de Tafira Baja, $35017$ Las Palmas de Gran
Canaria, Spain.Department of Mathematics, Wroclaw University of Environmental
and Life Sciences, ul. C.K. Norwida 25, 50-275 Wroclaw, PolandFaculty of Pure
and Applied Mathematics, Wrocław University of Science and Technology, Poland,
corresponding author<EMAIL_ADDRESS>
###### Abstract

The time-fractional porous medium equation is an important model of many
hydrological, physical, and chemical flows. We study its self-similar
solutions, which make up the profiles of many important experimentally
measured situations. We prove that there is a unique solution to the general
initial-boundary value problem in the one-dimensional setting. When
supplemented with boundary conditions from the physical models, the problem
exhibits a self-similar solution described with the use of the Erdélyi-Kober
fractional operator. Using a backward shooting method, we show that there
exists a unique solution to our problem.

The shooting method is not only useful in deriving the theoretical results. We
utilize it to devise an efficient numerical scheme to solve the governing
problem along with two ways of discretizing the Erdélyi-Kober fractional
derivative. Since the latter is a nonlocal operator, its numerical realization
has to include some truncation. We find the correct truncation regime and
prove several error estimates. Furthermore, the backward shooting method can
be used to solve the main problem, and we provide a convergence proof. The
main difficulty lies in the degeneracy of the diffusivity. We overcome it with
some regularization. Our findings are supplemented with numerical simulations
that verify the theoretical findings.

Keywords: time-fractional porous medium equation, Erdélyi-Kober fractional
operator, numerical method.

AMS Classification: 34A08, 65M12, 76S05.

##  1 Introduction

Our main interest is the following time-fractional porous medium problem on
the half-line

$\begin{cases}\partial_{t}^{\alpha}u=\left(D(u)u_{x}\right)_{x},&x\in\mathbb{R}_{+},\;t\in(0,T),\;\alpha\in(0,1)\\\
u(x,0)=0,&x\in\mathbb{R}_{+},\\\ u(0,t)=M,&t\in(0,T),\end{cases}$ (1)

where $T>0$ is the final time, $M>0$ the initial value, while the Caputo
fractional derivative is defined with the help of the fractional integral
$I^{\alpha}_{t}$

$\partial_{t}^{\alpha}u(x,t):=I^{1-\alpha}_{t}u_{t}(x,t),\quad
I^{\alpha}_{t}u(x,t)=\frac{1}{\Gamma(\alpha)}\int_{0}^{t}(t-s)^{\alpha-1}u(x,s)ds,\quad\alpha\in(0,1).$
(2)

For the diffusivity $D$ we assume that it is a $C^{1}(\mathbb{R}_{+})$
function with

$D(0)=0,\quad D(u)>0,\quad D^{\prime}(u)>0,\quad u>0.$ (3)

The most important example is the typical power type (or Brookes-Correy model,
as known in hydrology [4]) diffusivity for which $D(u)\propto u^{m}$ with
$m>0$. Note that we assume the _degeneracy_ , that we allow for a situation
when $D(0)=0$. This has a profound effect on the solution and is the main
reason for the finite speed of propagation (for a comprehensive treatment see
[49, 50]). The mathematical treatment of the considered PDE (1) has attracted
some recent attention, especially for the time-fractional case. For example,
bounded weak solutions of the degenerate and non-degenerate cases have been
found in [52] in the case of vanishing Dirichlet data and the problem on a
bounded domain $x\in\Omega\subset\mathbb{R}$. A very general abstract study of
the associated problem has also been given in [1]. Moreover, in [7] a similar
problem has been considered in the full space $x\in\mathbb{R}^{d}$ and the
authors have proved the existence and uniqueness of a complactly supported
solution provided that the initial data have this property. It seems that the
problem on the half-line has not been investigated adequately in the time-
fractional setting. In this paper, we present some further analytical and
numerical results that go beyond our initial works [41, 36, 34]. Our main
observation is based on the fact that the considered initial and boundary data
are self-similar, which allow for a transformation of the governing PDE into
an ordinary integro-differential equation. The main evolution operator then
becomes the Erdélyi-Kober (EK) fractional operator that has previously been
found in complex analysis and special functions [17, 45]. Furthermore, when
analyzing stochastic processes derived from fractional Brownian motion, the EK
operator becomes the main driving force to evolve diffusive dynamics [33].

The problem (1) models moisture imbibition in the essentially one-dimensional
porous medium with the boundary. The initially dry domain is suddenly exposed
to a constant concentration of moisture at the boundary. This is a typical
setting for measuring the properties of the medium [51]. In this setup, the
obvious type of solution to look for is the self-similar profile, which is a
function of the self-similar variable $x/t^{1/2}$. However, as some new
experiments on modern materials show, diffusion can undergo slower (sub-
diffusive) and faster (super-diffusive) dynamics [9, 19, 28, 42, 53, 47, 8,
55]. Further examples come from biology [22, 48], physics [44, 5], and other
fields of science. In [12] it was suggested that the time-fractional Caputo
derivative is adequate for modeling such a phenomenon. This approach was
successful in the sense that the modeling outcome reproduced the experimental
data with good precision. A physical derivation of the model in the
hydrological setting is given in [35].

There is a substantial amount of various numerical methods for the diffusion
with the time-fractional derivative (for a review, see [6, 11, 23]). We would
like to stress the fact that most of the results consider linear and only
space-dependent diffusivity $D=D(x)$. The approaches are based on different
kinds of Caputo discretization method and spatial scheme. For example, the
reader can consult several approaches in [18, 21] for the so-called L1 Caputo
discretization scheme and convolution quadrature [43, 13]. Finite difference
methods were considered in [46], finite elements in [25, 10], and spectral
methods in [26, 24]. Finding the error estimates for the case with nonsmooth
initial data requires some additional care and it is crucial for many
applications. This is especially relevant for time-dependent diffusivity [15,
29, 37]. As for some excursions from linearity, there are some recent methods
concerning the semilinear diffusion in [14, 2, 16]. The quasilinear case is
just beginning to be investigated, and we can refer the reader to our previous
work concerning this important problem [39, 38]. We would like to stress the
fact that most of these results were considered only the non-degenerate case
for which $D\geq D_{0}>0$. We have developed several numerical approaches for
the time-fractional degenerate case in some of our previous work. For example,
in [36] a fast quadrature was devised to solve 1. This result was further
refined to yield a second-order scheme in [31, 32]. Since the nonlocal nature
of the Caputo (or EK) operator increases the computational cost of all of the
numerical methods used to approximate them, we believe that there is a need
for developing fast and accurate schemes for solving equations governed by
these operators, especially in the degenerate case.

This paper has the following structure. In the next section we give a short
proof of the uniqueness of a general time-fractional porous medium problem in
$\Omega\subseteq\mathbb{R}$, where $\Omega$ can be bounded or not. The idea of
the proof is to modify the classical approach to the Caputo derivative
setting. Having the uniqueness, we proceed in Section 3 to the study of self-
similar solutions of the main problem (1). We use the backward shooting method
to prove the existence of such solutions and combine it with the previous
uniqueness result. In Section 4 we use the hints of our analytical reasoning
to devise efficient numerical methods for approximating the Erdélyi-Kober
fractional operator and the exact solution of (1). There we also present
several error estimates and the convergence proof. We illustrate the theory by
several numerical examples.

##  2 Uniqueness of the weak solution

Although we are mostly interested in (1) as the model of moisture imbibition
in porous media, in this section we consider a general problem for which we
present the proof of uniqueness. To this end, consider

$\begin{cases}\partial_{t}^{\alpha}u=\left(D(u)u_{x}\right)_{x},&x\in(a,b),\;t\in(0,T),\;\alpha\in(0,1)\\\
u(x,0)=u_{0}(x),&x\in(a,b),\\\
u(a,t)=\mu(t),\;u(b,t)=\nu(t),&t\in(0,T),\end{cases}$ (4)

in which we allow for a general initial and boundary conditions. Define
$\Omega_{T}:=(a,b)\times(0,T)$ as the domain of the solution. Moreover, let

$K(z):=\int_{0}^{z}D(s)ds,$ (5)

then the PDE can be written as $\partial^{\alpha}_{t}u=(K(u))_{xx}$. As in the
classical case, we cannot expect that the above will enjoy classical solutions
but rather weak ones (although there are other options [50, 52]). Let $\chi\in
C^{1}(\Omega_{T})$ be the arbitrary test function. By multiplication and
integration of the PDE we obtain the following

$\int_{a}^{b}\int_{0}^{T}\partial^{\alpha}_{t}u\;\chi
dtdx-\int_{a}^{b}\int_{0}^{T}(K(u))_{xx}\chi dtdx=0.$ (6)

Now, we integrate by parts to move the derivatives into test functions. First,
for the time derivative by (2) and Fubini’s theorem we have

$\begin{split}\int_{0}^{T}\partial^{\alpha}_{t}u\;\chi
dt&=\int_{0}^{T}(I^{1-\alpha}_{t}u_{t})\chi
dt=\frac{1}{\Gamma(1-\alpha)}\int_{0}^{T}\left(\int_{0}^{t}(t-s)^{-\alpha}u_{s}(x,s)ds\right)\chi(x,t)dt\\\
&=\frac{1}{\Gamma(1-\alpha)}\int_{0}^{T}\left(\int_{s}^{T}(t-s)^{-\alpha}\chi(x,t)dt\right)u_{s}(x,s)ds=\int_{0}^{T}u_{s}(J^{1-\alpha}_{s}\chi)ds,\end{split}$
(7)

where we define the right-sided fractional integral,

$J^{\alpha}_{t}u(x,t)=\frac{1}{\Gamma(\alpha)}\int_{t}^{T}(s-t)^{\alpha-1}u(x,s)ds,\quad\alpha>0.$
(8)

Now, by integrating by parts and renaming the integration variable, we can
obtain

$\int_{0}^{T}\partial^{\alpha}_{t}u\;\chi
dt=u_{0}(x)J^{1-\alpha}_{t}\chi(x,0)-\int_{0}^{T}u\frac{\partial}{\partial
t}(J^{1-\alpha}_{t}\chi)dt.$ (9)

For the space derivatives the integration by parts along with boundary
conditions gives the following

$\int_{a}^{b}(K(u))_{xx}\chi
dx=D(\nu(t))\chi(b,t)-D(\mu(t))\chi(a,t)-\int_{a}^{b}(K(u))_{x}\chi_{x}dx.$
(10)

Therefore, we define the _weak solution_ of (4) as a $H^{1}(\Omega_{T})$
function that satisfies

$\begin{split}\int_{0}^{\infty}\int_{0}^{T}&\left[(K(u))_{x}\chi_{x}-u\frac{\partial}{\partial
t}(J^{1-\alpha}_{t}\chi)\right]dtdx\\\
&=\int_{0}^{T}\left[D(\mu(t))\chi(a,t)-D(\nu(t))\chi(b,t)\right]dt-\int_{a}^{b}u_{0}(x)J^{1-\alpha}_{t}\chi(x,0)dx,\quad\chi\in
C^{1}(\Omega_{T}).\end{split}$ (11)

This approach to defining the weak solution mimics the classical case with
$\alpha=1$. Some relevant existence results are given in [50].

We now turn to the uniqueness proof. It is a generalization of the simple and
elegant argument originally devised by Ladyzenskaya [20] and later frequently
used in many cases [50].

######  Theorem 1.

There can be at most one solution to (11).

###### Proof.

Assume that $u_{1}$ and $u_{2}$ are some solutions to (11). Taking the
difference cancels the initial and boundary conditions, yielding

$\int_{a}^{b}\int_{0}^{T}\left[\left((K(u_{1}))_{x}-(K(u_{2}))_{x}\right)\chi_{x}-(u_{1}-u_{2})\frac{\partial}{\partial
t}(J^{1-\alpha}_{t}\chi)\right]dtdx=0$ (12)

for each $\chi\in C^{1}(\Omega_{T})$. Now, choose a test function of the form

$\chi=J^{\alpha}_{t}(K(u_{1})-K(u_{2})).$ (13)

Of course, the above choice may have to have sufficient regularity in order to
make it a test function; however, by a standard mollification argument, we can
ascertain that it is admissible. The details of such a procedure are described
in detail in [50] and, hence, we omit them here.

The integral with the gradient is non-negative. This can be seen by using
Fubini’s theorem just as in (7) to move the right-sided fractional integral
$J^{\alpha}$ into the usual (left-sided) integral $I^{\alpha}$ and invoking
Lemma 3.1 (ii) from [30]

$\begin{split}\int_{a}^{b}\int_{0}^{T}(K(u_{1})-K(u_{2}))_{x}J^{\alpha}_{t}(K(u_{1})-K(u_{2}))_{x}dtdx=\\\
\int_{a}^{b}\int_{0}^{T}(K(u_{1})-K(u_{2}))_{x}I^{\alpha}_{t}(K(u_{1})-K(u_{2}))_{x}dtdx&=\int_{a}^{b}\int_{0}^{T}|I^{\alpha/2}(K(u_{1})-K(u_{2}))_{x}|^{2}dtdx\geq
0.\end{split}$ (14)

Now, in the other integral we use the semigroup property of the fractional
integral, that is, $J^{\mu}J^{\nu}=J^{\nu+\mu}=J^{\nu}J^{\nu}$ to obtain the
following

$\begin{split}-&\int_{a}^{b}\int_{0}^{T}(u_{1}-u_{2})\frac{\partial}{\partial
t}(J^{1-\alpha}_{t}J^{\alpha}_{t}(K(u_{1})-K(u_{2})))dtdx\\\
&=-\int_{a}^{b}\int_{0}^{T}(u_{1}-u_{2})\frac{\partial}{\partial
t}(J^{1}(K(u_{1})-K(u_{2})))dtdx=\int_{a}^{b}\int_{0}^{T}(u_{1}-u_{2})(K(u_{1})-K(u_{2}))dtdx\geq
0\end{split}$ (15)

since $K$ is increasing, and therefore $K(u_{1})-K(u_{2})$ has the same sign
as $u_{1}-u_{2}$. Therefore, (12) transforms into

$\int_{a}^{b}\int_{0}^{T}\left[|I^{\alpha/2}(K(u_{1})-K(u_{2}))_{x}|^{2}+(u_{1}-u_{2})(K(u_{1})-K(u_{2}))\right]dtdx=0.$
(16)

Because the integrand is non-negative we must have $u_{1}=u_{2}$. This ends
the proof. ∎

##  3 Self-similar solution

Now we can go back to the original problem (1) and look for the solution in a
self-similar form [41, 36]

$u(x,t)=U(\eta),\quad\eta:=xt^{-\frac{\alpha}{2}}.$ (17)

Standard calculations (see, for example, [36]) lead to the ordinary equation
for the unknown profile $U=U(\eta)$

$\left(D(U)U^{\prime}\right)^{\prime}=\left[A-B\eta\frac{d}{d\eta}\right]F_{\alpha}U,\quad
0<\eta<\infty,$ (18)

where $A=1-\alpha$ and $B=\alpha/2$. Note that for different types of boundary
conditions, we obtain different values of the constants $A$ and $B$, however,
the structure of the above equation stays the same. This is why we decided to
leave general constants appearing in the governing equation. The operator
$F_{\alpha}$ is a particular version of the Erdlélyi-Kober fractional operator
[45]

$F_{\alpha}U(\eta)=\frac{1}{\Gamma(1-\alpha)}\int_{0}^{1}(1-s)^{-\alpha}U(s^{-B}\eta)ds,\quad
0<\alpha<1.$ (19)

The boundary conditions are

$U(0)=M,\quad U(\infty)=0.$ (20)

By Theorem 1 we know that the original problem (1) has a unique solution and
it is of the self-similar form. In [41] it has been proved that for the
diffusivity of the porous medium, that is $D(U)\propto U^{m}$ for $m\geq 1$,
there exists a _compactly supported_ solution. We now know that it is
precisely the unique solution of (1). What remains is to consider the general
case $D=D(U)$ satisfying (3). In what follows, we present some initial results
on this topic.

There is no straightforward way to solve (18). In [41] we have adopted a
certain transformation in the case of power-type diffusivity to obtain a
Volterra integral equation for which the theory is known. Since in the general
case no such transformation is available, we follow a different route. The
idea is to use _backward shooting method_ \- the idea that was used in the
classical setting [3]. However, in this non-local version, it gains much more
depth and meaning. To wit, we assume that we are looking for a completely
supported solution with support $[0,\eta^{*}]$ for fixed $\eta^{*}>0$. That
is, we have $U(\eta)=0$ for $\eta\geq\eta^{*}$. We then consider _initial
value problem_ starting from $\eta=\eta^{*}$ going backward to $\eta=0$.
Finally, by varying $\eta^{*}$ we can adjust the value of the solution at
$\eta=0$ in order to have $U(0)=M$. This is precisely the idea of the
numerical method presented in the next section. The crucial observation is to
note that in the definition of the EK operator (19) the function that is
operated is evaluated at $s^{-\alpha/2}\eta$ for fixed $\eta>0$. But since,
for $s\in(0,1)$ we always have $s^{-\alpha/2}>1$ it appears that imposing an
initial value would not yield an amenable problem. That is to say, in order to
compute $F_{\alpha}U(\eta)$ for $\eta$ close to $0$ we would require knowledge
of $U$ over almost the whole half-line, which is certainly not admissible in a
step-by-step calculation. In other words, we cannot advance our solution from
$\eta=0$ onward.

First, let us observe how the solution to our problem can behave. The
following result states many a priori properties of the solution.

######  Proposition 1.

Fix $\eta^{*}>0$. Let $U=U(\eta)$ be a solution of (18) with (3) in the left
neighborhood of $\eta^{*}$ such that $U(\eta)=0$ for $\eta\geq\eta^{*}$ .
Then, it is positive, decreasing, and the following holds

$\int_{0}^{1}\frac{D(s)}{s}ds<\infty,$ (21)

along with

$\lim\limits_{\eta\rightarrow\eta^{*}}-D(U)U^{\prime}(\eta)=0.$ (22)

###### Proof.

Since $U(\eta)=0$ for $\eta\geq\eta^{*}$ we can integrate (18) from $\eta$ to
$\eta^{*}$ to obtain

$-D(U)U^{\prime}(\eta)=A\int_{\eta}^{\eta^{*}}F_{\alpha}U(z)dz-B\int_{\eta}^{\eta^{*}}z\left(F_{\alpha}U\right)^{\prime}(z)dz=(A+B)\int_{\eta}^{\eta^{*}}F_{\alpha}U(z)dz+B\eta
F_{\alpha}U(\eta),$ (23)

where in the second equality, we have calculated by parts and used the
definition of the EK operator (19). From this it immediately follows that
$U^{\prime}(\eta)<0$ and $U(\eta)\geq 0$ which by (19) implies that

$F_{\alpha}U(\eta)\leq\frac{U(\eta)}{\Gamma(2-\alpha)}.$ (24)

Moreover, letting $\eta\rightarrow\eta^{*}$ we obtain the no-flux condition
(22). Hence, using the monotonicity and going back to the integrated equation,
we have the following

$\begin{split}-D(U)U^{\prime}(\eta)&\leq\frac{A}{\Gamma(2-\alpha)}\eta^{*}U(\eta)-B\eta^{*}\int_{\eta}^{\eta^{*}}\left(F_{\alpha}U\right)^{\prime}(z)dz\\\
&=\frac{A}{\Gamma(2-\alpha)}\eta^{*}U(\eta)+B\eta^{*}F_{\alpha}U(\eta)\leq\frac{A+B}{\Gamma(2-\alpha)}\eta^{*}U(\eta).\end{split}$
(25)

If we now divide by $U(\eta)$ and integrate from arbitrary $\eta_{1}$ to
$\eta_{2}$ we obtain

$-\int_{\eta_{1}}^{\eta_{2}}\frac{D(U(z))U^{\prime}(z)}{U(z)}dz\leq\frac{A+B}{\Gamma(2-\alpha)}\eta^{*}\left(\eta_{2}-\eta_{1}\right)\leq\frac{A+B}{\Gamma(2-\alpha)}(\eta^{*})^{2}.$
(26)

A change of the variable $s=U(z)$ along with the monotonicity of $U$ lets us
write

$\int_{U(\eta_{2})}^{U(\eta_{1})}\frac{D(s)}{s}ds<\frac{A+B}{\Gamma(2-\alpha)}(\eta^{*})^{2}.$
(27)

Letting $\eta_{2}\rightarrow\eta^{*}$ implies $U(\eta_{2})\rightarrow 0$ which
concludes the proof. ∎

From the above proof, we see that our compactly supported solution confirms
every physical intuition: it is a positive, bounded function with finite speed
of propagation and vanishing flux at the interface. The condition (21) is
necessary for the existence of the compact support. For example, with power-
type diffusion, we have

$\int_{0}^{1}\frac{s^{m}}{s}ds<\infty\quad\text{iff}\quad m>0,$ (28)

which confirms our previous results.

To proceed further we integrate (18) twice from $\eta$ to $\eta^{*}$, use the
vanishing boundary conditions at $\eta=\eta^{*}$, and integrate by parts to
arrive at the integral equation

$K(U(\eta))=\int_{\eta}^{\eta^{*}}\left((A+B)(z-\eta)+Bz\right)F_{\alpha}U(z)dz=:\int_{\eta}^{\eta^{*}}G(\eta,z)F_{\alpha}U(z)dz,$
(29)

In the following we present the main existence result.

######  Theorem 2.

Fix $\eta^{*}>0$. There exists a solution to (18) with $U(\eta)=0$ for
$\eta\geq\eta^{*}$ and (22).

###### Proof.

We will apply the Leray-Schauder fixed point theorem. First, since $D$ is
increasing, the function $K$ defined in (5) is convex. Therefore, there exists
a unique positive solution $x_{\lambda}$ to the equation $K(x)=\lambda x$ for
any $\lambda>0$. Let $X=C[0,\eta^{*}]$ be the Banach space of continuous
functions on $[0,\eta^{*}]$ with the norm
$\|U\|:=\max_{0\leq\eta\leq\eta^{*}}|U(\eta)|$. From (29) we can obtain the a
priori bound for the solution. First, by the fundamental estimate of the
Erdélyi-Kober operator (24) we have the following

$K(U(\eta))\leq(A+2B)\eta^{*}\int_{0}^{\eta^{*}}F_{\alpha}U(z)dz\leq\frac{A+2B}{\Gamma(2-\alpha)}(\eta^{*})^{2}\|U\|.$
(30)

By taking the maximum on the left-hand side and using the continuity of $U$ we
further have the following

$K(\|U\|)\leq\frac{A+2B}{\Gamma(2-\alpha)}(\eta^{*})^{2}\|U\|.$ (31)

Since $K$ is convex, from simple geometrical considerations, we must have
$\|U\|\leq x_{\lambda}$ with $\lambda:=(A+2B)/\Gamma(2-\alpha)(\eta^{*})^{2}$.
Therefore, we have the a priori upper bound for any solution to (18).

Having the bound for the solution, we define the operator $N:X\mapsto X$ by
the formula

$N(y)(\eta)=\int_{\eta}^{\eta^{*}}G(\eta,z)F_{\alpha}K^{-1}y(z)dz,$ (32)

which is well defined because $K$ is monotone and hence $K^{-1}$ exists. If
$y$ is the fixed point of $N$, then $U=K^{-1}y$ will be the solution of (18).
Since the integrand in the definition of $N$ is a continuous function of
$\eta$, $z$, and $y$ defined in a bounded and closed set ($y$ is bounded), the
operator $N$ is compact. Therefore, by the standard version of the Leray-
Schauder theorem (for ex. Theorem 6.A in [54]), that is, a priori bounded and
compact operator has a fixed point, we conclude that problem (18) has a
solution. ∎

Now we know that for each $\eta^{*}$ there exists a complactly supported
solution of our problem that is a bounded decreasing function. In the
following we show that, at least for small $U(0)=M$, we can determine that
there is a $\eta^{*}$ such that the solution attains $M$ for $\eta=0$.

######  Proposition 2.

Assume (3) and (21). For sufficiently small $M>0$ there exists a unique
solution to equation (18) with $U(0)=M$ and $U(\eta)=0$ for
$\eta\geq\eta^{*}$.

###### Proof.

Define the continuous function $f(\eta^{*}):=U(0;\eta^{*})$. Our goal is to
show that the function $\eta^{*}\mapsto f(\eta^{*})-M$ has exactly one zero.
If we take $\eta_{1}\rightarrow 0$ and $\eta_{2}\rightarrow\eta^{*}$ in (27)
we obtain

$\int_{0}^{f(\eta^{*})}\frac{D(s)}{s}ds\leq\frac{A+B}{\Gamma(2-\alpha)}(\eta^{*})^{2}.$
(33)

Therefore, by assumption we have $f(\eta^{*})\rightarrow 0$ for
$\eta^{*}\rightarrow 0$. If $M$ is small enough, then there exists $\eta^{*}$
such that $f(\eta*)=M$ and the existence is proved.

In order to ascertain uniqueness, we will show that $f$ is an increasing
function. To this end, assume on the contrary that there are
$\eta_{1}<\eta_{2}$ such that $f(\eta_{1})\geq f(\eta_{2})$. From the
monotonicity of the solution $U$ it follows that there exist $\eta_{0}$ such
that $U_{1}(\eta_{0})=U_{2}(\eta_{0})$ and $U_{1}(\eta)<U_{2}(\eta)$ for
$\eta_{0}<\eta\leq\eta_{1}$. Here, we have denoted $U_{i}$ as the solution to
the problem with the support $[0,\eta_{i}]$. We find that each $U_{i}$
satisfies (29). Subtracting the respective equations from $\eta=\eta_{0}$ we
obtain the following

$0=\int_{\eta_{0}}^{\eta_{1}}G(\eta_{0},z)\left(F_{\alpha}U_{2}(z)-F_{\alpha}U_{1}(z)\right)dz+\int_{\eta_{1}}^{\eta_{2}}G(\eta_{0},z)F_{\alpha}U_{2}(z)dz.$
(34)

However, both above terms are strictly positive and, hence, we arrive at a
contradiction. The function $f$ is strictly increasing. ∎

Therefore, assuming (3) and (27) we have proved that there exists a self-
similar solution to the problem (1) and by Theorem 1 we know that it is
unique.

##  4 Numerical methods

In this section, we develop an efficient numerical method for both
approximating the Erdélyi-Kober fractional operator (19) and solving the main
time-fractional porous medium equation (1).

###  4.1 Erdelyi-Kober operator

We start by numerically approximating the EK operator (19) acting on any
bounded function $U:[0,\infty)\mapsto\mathbb{R}_{+}$. Note that for a moment
we do not assume that $U$ satisfies (29) nor has a compact support. Some
general quadratures for the EK operator have been analyzed in [40] where a
thorough error analysis has also been given. Here, we focus only on the main
difference between the paper cited: here, our operator involves the solution
evaluated at $s^{-B}\eta$. In [40] only the positive exponent case was
considered. This furnishes a radical change in the numerical analysis. To
provide a concrete examples, we devise two schemes: of first (rectangle) and
second (trapezoid) order. First, introduce a uniform grid of with a step $h>0$

$\eta_{n}=nh,\quad n\geq 0.$ (35)

If $n=0$, we immediately see from (19) that

$F_{\alpha}U(0)=\frac{U(0)}{\Gamma(2-\alpha)},$ (36)

therefore, we must focus on the case $n\geq 1$. To this end, change the
integration variable in the definition (19) according to $z=s^{-B}\eta$ to
obtain

$F_{\alpha}U(\eta_{n})=\frac{\eta^{\frac{1-\alpha}{B}}}{B\Gamma(1-\alpha)}\int_{\eta_{n}}^{\infty}(\eta^{-\frac{1}{B}}-z^{-\frac{1}{B}})^{-\alpha}z^{-\frac{1}{B}-1}U(z)dz,$
(37)

where now the forward-nonlocal property of the EK operator is evident. That is
to say, the value $F_{\alpha}U(\eta_{n})$ depends on $U(z)$ for
$z\geq\eta_{n}$. Since the integral is improper, in order to evaluate it
numerically, we have to truncate it at some point, say $\eta_{N}$. We will
choose the optimal value for this truncation later. Therefore,

$F_{\alpha}U(\eta_{n})=\frac{\eta^{\frac{1-\alpha}{B}}}{B\Gamma(1-\alpha)}\int_{\eta_{n}}^{\eta_{N}}(\eta^{-\frac{1}{B}}-z^{-\frac{1}{B}})^{-\alpha}z^{-\frac{1}{B}-1}U(z)dz+R_{N}(\eta_{n}),$
(38)

with the remainder

$R_{N}(\eta_{n})=\frac{\eta^{\frac{1-\alpha}{B}}}{B\Gamma(1-\alpha)}\int_{\eta_{N}}^{\infty}(\eta^{-\frac{1}{B}}-z^{-\frac{1}{B}})^{-\alpha}z^{-\frac{1}{B}-1}U(z)dz.$
(39)

Now, we can write

$F_{\alpha}U(\eta_{n})=\frac{\eta^{\frac{1-\alpha}{B}}}{B\Gamma(1-\alpha)}\sum_{i=n}^{N-1}\int_{\eta_{i}}^{\eta_{i+1}}(\eta^{-\frac{1}{B}}-z^{-\frac{1}{B}})^{-\alpha}z^{-\frac{1}{B}-1}U(z)dz+R_{N}(\eta_{n}),$
(40)

and approximate the function $U$ in the small interval
$[\eta_{i},\eta_{i+1})$. The two simplest choices are the rectangle and
trapezoid approximation for which

$U(z)\approx U(\eta_{i+1}),\quad z\in[\eta_{i},\eta_{i+1}),$ (41)

and

$U(z)\approx U(\eta_{i})+\frac{U(\eta_{i+1})-U(\eta_{i})}{h}(z-\eta_{i}),\quad
z\in[\eta_{i},\eta_{i+1}),$ (42)

respectively. Plugging the above into the EK integral reveals that

$F_{\alpha}U(\eta_{n})\approx\widehat{F}_{\alpha,N}U(\eta_{n}):=\sum_{i=n}^{N}a^{(r,t)}_{in}U(\eta_{i}),$
(43)

with the following _positive_ weights that can be computed by a
straightforward calculation

$a^{(r)}_{in}=\begin{cases}\dfrac{1}{\Gamma(2-\alpha)},&i=n=0,\\\
0,&n=0,\;0<i\leq N,\\\ 0,&i=n>0,\\\
\dfrac{1}{\Gamma(2-\alpha)}\left[\left(1-\left(\dfrac{i}{n}\right)^{-\frac{1}{B}}\right)^{1-\alpha}-\left(1-\left(\dfrac{i-1}{n}\right)^{-\frac{1}{B}}\right)^{1-\alpha}\right],&n<i\leq
N,\;n>0.\end{cases}$ (44)

and

$\begin{split}a^{(t)}_{in}&=\begin{cases}\dfrac{1}{\Gamma(2-\alpha)},&i=n=0,\\\
0,&n=0,\;0<i\leq N,\\\ a^{(r)}_{(n+1)n}-d_{nn},&i=n>0,\\\
d_{(i-1)n}-d_{in}+a^{(r)}_{(i+1)n},&n<i\leq N-1,\;n>0,\\\
d_{(M-1)n},&i=N,\;n>0,\end{cases}\\\
d_{in}&=\frac{n}{\Gamma(1-\alpha)}\left[\beta\left(\left(\frac{i}{n}\right)^{-\frac{1}{B}};1-B,1-\alpha\right)-\beta\left(\left(\frac{i+1}{n}\right)^{-\frac{1}{B}};1-B,1-\alpha\right)\right]-ia^{(r)}_{(i+1)n},\end{split}$
(45)

where the superscripts $(r)$ and $(t)$ denote the rectangle and trapezoid
rules, respectively. Here, $\beta(z;a,b)$ is the Euler incomplete beta
function. The following result gives the error bounds for the discretization
operator $\widehat{F}_{\alpha,N}$.

######  Theorem 3.

Let $U:[0,\infty)\mapsto\mathbb{R}_{+}$ be a sufficiently bounded and smooth
function. Moreover, set $N=\gamma n$, where

$\gamma=\begin{cases}[h^{-B}]+1,&\text{rectangle quadrature},\\\
[h^{-2B}]+1,&\text{trapezoid quadrature}.\end{cases}$ (46)

Then, we have the following error bounds

$\|F_{\alpha}U(\eta_{n})-\widehat{F}_{\alpha,N}U(\eta_{n})\|_{\infty}\leq\begin{cases}\left(\max_{z\geq\eta_{N}}|U(z)|+\max_{0\leq
z\leq\eta_{N}}|U^{\prime}(z)|\right)\dfrac{h}{\Gamma(2-\alpha)},&\text{rectangle
quadrature},\\\ \left(\max_{z\geq\eta_{N}}|U(z)|+\frac{1}{2}\max_{0\leq
z\leq\eta_{N}}|U^{\prime\prime}(z)|\right)\dfrac{h^{2}}{\Gamma(2-\alpha)},&\text{trapezoid
quadrature}.\end{cases}$ (47)

###### Proof.

We will prove only the rectangle case, the other is completely analogous. By
Taylor series, we immediately have

$U(z)=U(\eta_{i+1})+U^{\prime}(\zeta_{i})(\eta_{i+1}-z),\quad
z\in[\eta_{i},\eta_{i+1}),$ (48)

and plugging it into (40) yields

$F_{\alpha}U(\eta_{n})=\widehat{F}_{\alpha,N}U(\eta_{n})+P_{N,h}(\eta_{n})+R_{N}(\eta_{n}),$
(49)

where

$|P_{N,h}(\eta_{n})|\leq\max_{0\leq
z\leq\eta_{N}}|U^{\prime}(z)|\frac{\eta^{\frac{1-\alpha}{B}}}{B\Gamma(1-\alpha)}\sum_{i=n}^{N-1}\int_{\eta_{i}}^{\eta_{i+1}}(\eta^{-\frac{1}{B}}-z^{-\frac{1}{B}})^{-\alpha}z^{-\frac{1}{B}-1}(\eta_{i+1}-z)dz.$
(50)

But since $|\eta_{i+1}-z|\leq h$ we can compute the sum explicitly and change
back the integration variable $s^{-B}\eta_{n}=z$ to obtain the following

$\begin{split}|P_{N,h}(\eta_{n})|&\leq h\max_{0\leq
z\leq\eta_{N}}|U^{\prime}(z)|\frac{\eta^{\frac{1-\alpha}{B}}}{B\Gamma(1-\alpha)}\int_{\eta_{n}}^{\eta_{N}}(\eta^{-\frac{1}{B}}-z^{-\frac{1}{B}})^{-\alpha}z^{-\frac{1}{B}-1}(\eta_{i+1}-z)dz\\\
&=h\max_{0\leq
z\leq\eta_{N}}|U^{\prime}(z)|\frac{1}{\Gamma(1-\alpha)}\int_{(\frac{n}{N})^{1/B}}^{1}(1-s)^{-\alpha}ds=\frac{h}{\Gamma(2-\alpha)}\max_{0\leq
z\leq\eta_{N}}|U^{\prime}(z)|\left(1-\gamma^{-\frac{1}{B}}\right)^{1-\alpha}.\end{split}$
(51)

As for the truncation remainder (39) we can simply estimate

$\begin{split}|R_{N}(\eta_{n})|&\leq\max_{z\geq\eta_{N}}|U(z)|\frac{\eta^{\frac{1-\alpha}{B}}}{B\Gamma(1-\alpha)}\int_{\eta_{N}}^{\infty}(\eta^{-\frac{1}{B}}-z^{-\frac{1}{B}})^{-\alpha}z^{-\frac{1}{B}-1}dz\\\
&=\max_{z\geq\eta_{N}}|U(z)|\frac{1}{\Gamma(1-\alpha)}\int_{0}^{\left(\frac{n}{N}\right)^{1/B}}(1-s)^{-\alpha}ds=\frac{\max_{z\geq\eta_{N}}|U(z)|}{\Gamma(2-\alpha)}\left(1-\left(1-\gamma^{-\frac{1}{B}}\right)^{1-\alpha}\right).\end{split}$
(52)

Since

$\left(1-\gamma^{-\frac{1}{B}}\right)^{1-\alpha}\leq
1,\quad\left(1-\left(1-\gamma^{-\frac{1}{B}}\right)^{1-\alpha}\right)\leq\gamma^{-\frac{1}{B}},$
(53)

where the second inequality follows from convexity, by our assumption (46) on
$\gamma$ we have

$|P_{N,h}(\eta_{n})|+|R_{N}(\eta_{n})|\leq\left(\max_{z\geq\eta_{N}}|U(z)|+\max_{0\leq
z\leq\eta_{N}}|U^{\prime}(z)|\right)\frac{h}{\Gamma(2-\alpha)},$ (54)

which concludes the proof. ∎

As we can see, to obtain an optimal error, the truncation has to be chosen
according to the grid spacing $h$. The optimality in this sense is associated
with the same order of both remainders for $h\rightarrow 0^{+}$. Note also
that the higher the order of the quadrature, the larger the interval over
which we have to integrate. A numerical illustration of the above theorem can
be presented by choosing a function with an explicitly known EK operator. Let
$B=\alpha/2$

$U(\eta)=\min\left\\{1,\eta^{\mu}\right\\},\quad\mu>0,$ (55)

for which

$F_{\alpha}U(\eta)=\begin{cases}\frac{1-(1-\eta^{\frac{2}{\alpha}})^{1-\alpha}}{\Gamma(2-\alpha)}+\eta^{\mu}\left(\frac{\Gamma\left(1-\frac{\alpha\mu}{2}\right)}{\Gamma\left(2-\frac{\alpha(2+\mu)}{2}\right)}-\frac{\beta\left(\eta^{\frac{2}{\alpha}};1-\frac{\alpha\mu}{2},1-\alpha\right)}{\Gamma(1-\alpha)}\right),&0\leq\eta<1,\\\
\frac{1}{\Gamma(2-\alpha)},&\eta\geq 1.\end{cases}$ (56)

We can now easily compute the discretization error. In Fig. 1 we depict the
maximum error of approximating the EK operator with $\widehat{F}_{\alpha,N}$
with $N$ chosen according to the optimal choice (46). The error is plotted
with respect to the grid spacing $h$, and the respective orders of
approximation are clearly seen. As can be inferred, the graphs increase with a
slope corresponding to the quadrature order. However, we note that due to the
higher computational complexity of the trapezoid scheme in both function
evaluations and the larger $\gamma$, this method is more expensive for the
same $h$ compared to the simple rectangle quadrature.

Figure 1: The maximal discretization error with respect to the grid spacing
$h$ and $\alpha=0.5$. Reference lines are added for comparison.

###  4.2 Integro-differential equation

We can now proceed to discretization of the main equation (18). The strategy
is to consider the integral form (29) rather than the original and solve the
problem backwards. If some initial value is prescribed, say $U(0)=M$, we can
use the shooting method and look for the zero of a function $\eta^{*}\mapsto
M-U(0,\eta^{*})$.

For what follows, we fix the value of the wetting front $\eta^{*}$. Since we
already know how to discretize the EK operator, it is just a matter of
choosing the correct quadrature for the integral in (29). Similarly, as
before, we choose either the rectangle or trapezoid methods. To this end, we
naturally choose the integration horizon to $\eta^{*}$, that is, we choose $N$
in $\widehat{F}_{\alpha,N}$ according to

$\eta^{*}=\eta_{N}=Nh,$ (57)

by possibly adjusting $h$ so that $N$ is an integer. Thus a discretization of
the integral equation (29) can be found by splitting the integral into parts

$K(U(\eta_{n}))=\sum_{j=n}^{N-1}\int_{\eta_{j}}^{\eta_{j+1}}G(\eta_{n},z)F_{\alpha}U(z)dz.$
(58)

Now, approximating the EK operator on each subinterval by a constant or linear
function yields the following discretization

$K(U_{n})=\sum_{j=n}^{N}b^{(r,t)}_{jn}\sum_{i=n}^{N}a^{(r,t)}_{ij}U_{i}=\sum_{i=n}^{N}\left(\sum_{j=n}^{i}a^{(r,t)}_{ij}b^{(r,t)}_{jn}\right)U_{j}=:\sum_{i=n}^{N}c^{(r,t)}_{in}U_{i},$
(59)

where by $U_{i}$ we have denoted the numerical approximation to the exact
solution $U(\eta_{i})$. Weights $b^{(r,t)}_{jn}$ correspond to the usual
rectangle $(r)$ and trapezoid $(t)$ product quadratures for the integral. Note
that we do not approximate the kernel $G$ since it can be computed directly.
More specifically,

$b^{(r)}_{nn}=0,\quad
b^{(r)}_{jn}=\frac{h^{2}}{2}\left((A+2B)(2j-1)-2(A+B)n\right),\quad n<j\leq
N,$ (60)

and

$\begin{split}b^{(t)}_{jn}&=\begin{cases}b^{(r)}_{(n+1)n}-\widehat{d}_{nn},&j=n,\\\
\widehat{d}_{(j-1)n}-\widehat{d}_{jn}+b^{(r)}_{(j+1)n},&n<i\leq N-1,\\\
\widehat{d}_{(M-1)n},&j=N,\\\ \end{cases}\\\
\widehat{d}_{jn}&=\frac{h^{2}}{6}\left((A+2B)(3j-1)-3(A+B)n\right),\end{split}$
(61)

for rectangle and trapezoid quadratures, respectively. The scheme’s
coefficients $c_{in}$ can then be computed as the convolution of $a_{ij}$ and
$b_{jn}$, which can be done very efficiently with the use of Fast Fourier
Transform. Note from (44) and (60) that both rectangle weights vanish for
$i=j=n$. This means that the numerical scheme (59) is _explicit_ which makes
the method simple and fast. That is, to compute $U_{n}$ we have to know the
values of $U_{i}$ for $n<i\leq N$. Since $D$ is positive, its integral $K$ is
increasing and, hence, has a well-defined inverse. Therefore, we can write

$U_{n}=K^{-1}\left(\sum_{i=n+1}^{N}c^{(r)}_{in}U_{i}\right).$ (62)

This is especially relevant for the important case of power-law diffusivity
$K(u)=u^{m+1}/(m+1)$. Of course, when the analytical form of $K^{-1}$ is not
available, it is probably better to apply a root-finding algorithm to (59),
especially in the implicit trapezoid case. The overall procedure is now
$\eta$-stepping and computing the values of $U_{n}$ backwards for
$n=N-1,N-2,N-3,\cdots,1,0$.

Note that since the solution has a compact support terminating at
$\eta=\eta^{*}=\eta_{N}$, we always have

$U_{N}=0.$ (63)

The question arises how to start the scheme (59) in order _not to_ obtain the
trivial solution $U_{i}\equiv 0$. The answer is straightforward for the
trapezoid method, for we have

$K(U_{N-1})=c^{(t)}_{(N-1)(N-1)}U_{N-1}=a^{(t)}_{(N-1)(N-1)}b^{(t)}_{(N-1)(N-1)}U_{N-1},$
(64)

which is a nonlinear algebraic equation to be solved for $U_{N-1}$. The
explicit rectangle method does not have this form, and we have to prescribe
the starting value in a different way. To this end, let us return to the
integral equation (29) evaluated for $n=N-1$. If we use the constant function
approximation at the _left_ endpoint of the interval $(\eta_{N-1},\eta_{N})$
we obtain

$K(U(\eta_{N-1}))=\int_{\eta_{N-1}}^{\eta_{N}}G(\eta_{N-1},z)F_{\alpha}U(z)dz\approx
F_{\alpha}U(\eta_{N-1})\int_{\eta_{N-1}}^{\eta_{N}}G(\eta_{N-1},z)dz=F_{\alpha}U(\eta_{N-1})b^{(r)}_{N(N-1)}$
(65)

Doing the same left approximation in the rectangle quadrature for the EK
operator (43) yields the equation to be solved for the starting value
$U_{N-1}$

$K(U_{N-1})=a^{(r)}_{N(N-1)}b^{(r)}_{N(N-1)}U(\eta_{N-1}).$ (66)

Note the similarity with (64). For ease of implementation, the whole
$\eta$-stepping scheme is summarized in Algorithm 1.

Algorithm 1 Main $\eta$-stepping scheme for solving (18) with (LABEL:eqn:IC)
by the rectangle (r) or trapezoid (t) quadrature. The absolute tolerance is
given by a fixed $0<\epsilon<1$.

$N\in\mathbb{N}$

Define weights $a^{(r,t)}$, $b^{(r,t)}$, and $c^{(r,t)}$ (use FFT) (44), (45),
(59), (60), (61)

Define the function $K=K(z)$ by (5)

$U_{0:N}\leftarrow 0$

while $|U_{0}-1|\geq\epsilon$ do

Fix $\eta^{*}>0$ $\triangleright$ A root finding solver does that.

$h\leftarrow\eta^{*}/N$

Solve $K(U_{N-1})=a^{(r,t)}_{N(N-1)}b^{(r,t)}_{N(N-1)}U_{N-1}$ for $U_{N-1}$

for $n=N-2:-1:0$ do

Solve $K(U_{n})=\sum_{i=n}^{N-1}c^{(r,t)}_{in}U_{i}$ for $U_{n}$

end for

end while

Having described the numerical scheme, we can proceed to proving that it is
convergent. The main difficulty is the degeneracy of our equation, that is,
the fact that $D(0)=0$. In order to overcome it, we consider a family of
regularizations $D_{h}$ of the diffusivity that converge to $D$ as
$h\rightarrow 0$. For each $D_{h}$ we obtain a solution $U_{h,n}$ that
converges to $U$ as we refine the grid. Let $D_{h}$ be a family of functions
satisfying

$0<\epsilon(h)\leq D_{h}(z),\quad\epsilon(h)\rightarrow 0\quad\text{as}\quad
h\rightarrow 0,$ (67)

and

$\|D_{h}-D\|_{\infty}\leq\epsilon(h).$ (68)

The following result states the convergence proof.

######  Theorem 4.

Let $D_{h}$ be the family of regularizations of $D$ satisfying (67). Suppose
that the weights of the quadrature (59) satisfy $0<c_{in}\leq Ch$ and that the
and the quadrature of the integral in (5) have order $p>0$. Then, there exists
an $\epsilon(h)$ such that when $U_{h,n}$ is a solution of (59) with $D_{h}$,
and $U$ is an exact solution of (18) we have

$|U_{h,n}-U(\eta)|=O(h^{p-\delta}|\ln h|)\quad\text{when}\quad h\rightarrow
0,\quad nh\rightarrow\eta\in[0,\eta^{*}],$ (69)

where $0<\delta<1$ is arbitrary.

###### Proof.

Denote the error by $e_{n}:=U(\eta_{n})-U_{h,n}$ and let the quadrature error
for the integral (29) be denoted by $\rho_{n}(h)$, that is,

$\int_{\eta_{n}}^{\eta^{*}}G(\eta_{n},z)F_{\alpha}U(z)dz=\sum_{i=n}^{N-1}c_{in}U_{i}+\rho_{n}(h)\quad\text{where}\quad|\rho_{n}(h)|\leq\rho(h)\rightarrow
0,\quad h\rightarrow 0.$ (70)

Then, from this and (59) we have

$K_{h}(U(\eta_{n}))-K_{h}(U_{h,n})=\sum_{i=n}^{N-1}c_{in}e_{i}+\rho_{n}(h),$
(71)

where $K_{h}$ is the integral corresponding to $D_{h}$. Now, by the mean value
theorem we can write
$K_{h}(U(\eta_{n}))-K_{h}(U_{n,h})=K^{\prime}(V_{n})e_{n}=D_{h}(V_{n})e_{n}$
for some intermediate value $V_{n}$, and hence

$D_{h}(V_{n})|e_{n}|\leq\sum_{i=n}^{N-1}c_{in}|e_{i}|+\rho(h).$ (72)

Now, by the construction of $D_{h}$ we can write

$\epsilon(h)|e_{n}|\leq\sum_{i=n}^{N-1}c_{in}|e_{i}|+\rho(h),$ (73)

which is owing to the fact that $c_{in}\leq Ch$

$|e_{n}|\leq
Ch\epsilon(h)^{-1}\sum_{i=n}^{N-1}|e_{i}|+\epsilon(h)^{-1}\rho(h),$ (74)

which is a form amenable for the discrete version of the Grönwall inequality
(for ex. heorem 7.1 from [27] applied for $f_{i}:=e_{N-i}$). Therefore,

$|e_{n}|\leq\epsilon(h)^{-1}\rho(h)\left(1+Ch\epsilon(h)^{-1}\right)^{n}\leq\epsilon(h)^{-1}\rho(h)\left(1+\frac{Cnh\epsilon(h)^{-1}}{n}\right)^{n}\leq\epsilon(h)^{-1}\rho(h)e^{C\eta^{*}\epsilon(h)^{-1}},$
(75)

since $nh\rightarrow\eta\leq\eta^{*}$. Now, if we choose

$\epsilon(h)=C\eta^{*}\left(\delta\ln\frac{1}{h}\right)^{-1}\rightarrow
0\quad\text{as}\quad h\rightarrow 0,$ (76)

for some arbitrary $0<\delta<1$, we obtain

$|e_{n}|\leq\frac{\delta}{C\eta^{*}}\rho(h)h^{-\delta}|\ln h|.$ (77)

If now $\rho(h)=O(h^{p})$ with $p>0$ and $h\rightarrow 0$, we have

$|e_{n}|=O(h^{p-\delta}|\ln h|),\quad h\rightarrow 0,$ (78)

what finishes the proof. ∎

From the above proof we thus see that the order of the scheme for the
regularized solution is almost $p$, that is, the order of the quadrature for
(29). The actual order is less by an arbitrary small number $\delta$ and a
logarithmic factor.

We illustrate our theory by some numerical experiments. In what follows, we
always choose the rectangle scheme in approximating the solution. Our
simulations indicated that although the trapezoidal method is superior when
discretizing the pure EK operator (19) it is very expensive when applied to
the nonlinear equation (18). This computational cost comes from a large number
of special functions needed to calculate the weights of the trapezoid method
(45) - especially the incomplete beta function. In effect, the temporal and
spatial complexity of the algorithm can be prohibitively large. Moreover, the
method is implicit without a significant stability gain, and thus requires
solving a nonlinear equation in each iteration step. As a benchmark, we have
calculated the time ratio of computations needed to obtain the wetting front
position with the diffusivity $D(u)=u^{2}$ for number of subdivisions
$N=2^{8}$ for different values of $\alpha$. In Tab. 1 we present the quantity

$\tau=\frac{\text{time of computations for trapezoid method}}{\text{time of
computations for the rectangle method}}.$ (79)

Immediately we see that computations with the trapezoidal method are at least
_one hundred_ times slower than with the rectangle method. We can conclude
that the increase in accuracy for the second order method does not compensate
the high increase in computational cost. We have thus decided that a less
accurate but much faster explicit rectangle method will be the scheme of
choice. An efficient second-order explicit scheme for the power-law case, i.e.
$D(u)\propto u^{m}$ has been devised in [31] by different means that cannot be
generalized to the arbitrary diffusivity.

$\alpha$ | 0.1 | 0.25 | 0.5 | 0.75 | 0.9
---|---|---|---|---|---
$\tau$ | 91 | 230 | 168 | 197 | 380
Table 1: Time ratio $\tau$ defined in (79) for computing the wetting front
position with $D(u)=u^{2}$, $N=2^{8}$ with the trapezoid and rectangle method.
Figure 2: An exemplary plot of solutions to (18) with $U(0)=0$ for
$\alpha=0.5$ calculate with the Algorithm 1. Different diffusivities are
indicated in the legend.

In what follows we always solve (18) with the initial condition $U(0)=1$ with
the Algorithm 1 where we use Newton’s iteration for finding $\eta^{*}$. Some
exemplary plots of typical solutions of the porous medium equation are
depicted in Fig. 2 for several choices of the diffusivity

$D_{BC}(u)=u^{m},\quad D_{exp}(u)=1-e^{-u},\quad m\geq 1,$ (80)

where the first choice is the typical porous medium power-type diffusivity (in
hydrology known as the Brooks-Correy model) and the second is the exponential
formula. By a simple limit test, we can verify that the necessary condition
for the compact support (21) is satisfied for each of the diffusivities.

To illustrate convergence, we present two tests. The first one is an empirical
calculation of the convergence order using extrapolation (Aitken’s method)

$\text{order
}\approx\log_{2}\frac{|\eta^{*}_{(2N)}-\eta^{*}_{(N)}|}{|\eta^{*}_{(4N)}-\eta^{*}_{(2N)}|},$
(81)

in which we compare the wetting front (the worst-case scenario) calculated for
different numbers of subdivisions $N$ (hence, twice or quadruple smaller grid
spacing $h$). The results for two different diffusivities are presented in the
Tab. 2. As can be seen, the numerical estimates are consistent with the
theoretical predictions that the rectangle quadrature should retain its order.
This is not obvious since it is widely known that even for the _linear
diffusion_ , the discretizations of fractional derivatives may converge with
lower order than $1$ depending on the regularity of the solution (for a
comprehensive account, see [46]. We can see that for $\alpha=0.75$ the method
converged even faster.

$\alpha$ | $0.1$ | $0.25$ | $0.5$ | $0.75$ | $0.9$
---|---|---|---|---|---
order for $D_{BC}$ | 0.95 | 0.98 | 0.99 | 1.255 | 0.97
order for $D_{exp}$ | 0.96 | 0.98 | 0.99 | 1.08 | 0.98
Table 2: Estimated order of quadrature for the scheme (59) based on
extrapolation (81) applied for calculating the wetting front $\eta^{*}$. Two
diffusivities has been chosen from (80): $D_{BC}$ with $m=1$ and $D_{exp}$.
The base number of iteration is $N=300$.

The second test we conduct is once again based on finding the wetting front,
but not in the classical case, since then the exact values are available (see
[32]). The results are gathered in the Tab. 3. The results are decent based on
the fact that a small number of steps have been taken. The error decays at a
rate $O(N^{-1})$ as $N\rightarrow\infty$. This verifies the fact that our
method is convergent even in the classical case. However, if one wants to
conduct precision calculations on the wetting front for power type
diffusivity, it is recommended to use our second order method [31]. The
present one has the advantage of being fast and robust when it comes to
choosing a general form of diffusivity.

N | 10 | 50 | 100 | 200 | 500 | 1000
---|---|---|---|---|---|---
error | $5.5\times 10^{-2}$ | $3.3\times 10^{-2}$ | $2.2\times 10^{-2}$ | $1.3\times 10^{-2}$ | $7.0\times 10^{-3}$ | $4.0\times 10^{-3}$
Table 3: Absolute errors of calculating the wetting front $\eta^{*}$ for
different values of subdivisions of the interval $[0,\eta^{*}]$ for
$D_{B}C(u)=u$ and $\alpha=1$. The reference exact values were taken from [32].

##  5 Conclusion

The time-fractional porous medium equation models several important
experimental settings in material science, hydrology, and construction
engineering. We have proved that the problem with general diffusivity has a
unique solution that has a self-similar form. The main role was played here by
the Erdélyi-Kober fractional operator and its careful analysis. On the
practical side, we have devised a robust numerical method that can be easily
used by practitioners.

In our future work, we plan to resign from the small initial value requirement
and to consider a generalized version of (1) where we will allow for a
nonlocal in space operator. This will enlarge the number of possible modeling
situations and include the superdiffusive case, which has also been found in
many experiments.

## Acknowledgement

Ł.P. has been supported by the National Science Centre, Poland (NCN) under the
grant Sonata Bis with a number NCN 2020/38/E/ST1/00153.

## References

  * [1] Goro Akagi.  Fractional flows driven by subdifferentials in Hilbert spaces.  Israel Journal of Mathematics, 234(2):809–862, 2019.
  * [2] Mariam Al-Maskari and Samir Karaa.  Numerical approximation of semilinear subdiffusion equations with nonsmooth initial data.  SIAM Journal on Numerical Analysis, 57(3):1524–1544, 2019.
  * [3] FV Atkinson and LA Peletier.  Similarity profiles of flows through porous media.  Archive for Rational Mechanics and Analysis, 42(5):369–379, 1971\.
  * [4] Royal Harvard Brooks and Arthur Thomas Corey.  Hydraulic properties of porous media and their relationship to drainage design.  PhD thesis, Colorado State University. Libraries, 1963.
  * [5] Diego del Castillo-Negrete, BA Carreras, and VE Lynch.  Nondiffusive transport in plasma turbulence: a fractional diffusion approach.  Physical Review Letters, 94(6):065003, 2005.
  * [6] Kai Diethelm and Neville J Ford.  Analysis of fractional differential equations.  Journal of Mathematical Analysis and Applications, 265(2):229–248, 2002.
  * [7] Jean-Daniel Djida, Juan J Nieto, and Iván Area.  Nonlocal time-porous medium equation: weak solutions and finite speed of propagation.  Discrete Continuous Dyn. Syst. Ser. B, 2018.
  * [8] A El Abd, SE Kichanov, M Taman, KM Nazarov, DP Kozlenko, and Wael M Badawy.  Determination of moisture distributions in porous building bricks by neutron radiography.  Applied Radiation and Isotopes, 156:108970, 2020.
  * [9] Abd El-Ghany El Abd and Jacek J Milczarek.  Neutron radiography study of water absorption in porous building materials: anomalous diffusion analysis.  Journal of Physics D: Applied Physics, 37(16):2305, 2004.
  * [10] Neville J Ford, Jingyu Xiao, and Yubin Yan.  A finite element method for time fractional partial differential equations.  Fractional Calculus and Applied Analysis, 14(3):454–474, 2011.
  * [11] Roberto Garrappa.  Numerical solution of fractional differential equations: A survey and a software tutorial.  Mathematics, 6(2):16, 2018.
  * [12] DN Gerasimov, VA Kondratieva, and OA Sinkevich.  An anomalous non-self-similar infiltration and fractional diffusion equation.  Physica D: Nonlinear Phenomena, 239(16):1593–1597, 2010.
  * [13] Bangti Jin, Buyang Li, and Zhi Zhou.  Correction of high-order BDF convolution quadrature for fractional evolution equations.  SIAM Journal on Scientific Computing, 39(6):A3129–A3152, 2017.
  * [14] Bangti Jin, Buyang Li, and Zhi Zhou.  Numerical analysis of nonlinear subdiffusion equations.  SIAM Journal on Numerical Analysis, 56(1):1–23, 2018.
  * [15] Bangti Jin, Buyang Li, and Zhi Zhou.  Subdiffusion with a time-dependent coefficient: analysis and numerical solution.  Mathematics of Computation, 88(319):2157–2186, 2019.
  * [16] Samir Karaa.  Galerkin type methods for semilinear time-fractional diffusion problems.  Journal of Scientific Computing, 83(3):1–22, 2020.
  * [17] Virginia S Kiryakova and Bader N Al-Saqabi.  Transmutation method for solving Erdélyi–Kober fractional differintegral equations.  Journal of Mathematical Analysis and Applications, 211(1):347–364, 1997.
  * [18] Natalia Kopteva.  Error analysis of the L1 method on graded and uniform meshes for a fractional-derivative problem in two and three dimensions.  Mathematics of Computation, 88(319):2135–2155, 2019.
  * [19] Michel Küntz and Paul Lavallée.  Experimental evidence and theoretical analysis of anomalous diffusion during water infiltration in porous building materials.  Journal of Physics D: Applied Physics, 34(16):2547, 2001.
  * [20] Olga Aleksandrovna Ladyzenskaya.  On the linear and quasilinear parabolic equations.  Differential Equations and Their Applications, pages 273–279, 1967\.
  * [21] TAM Langlands and Bruce I Henry.  The accuracy and stability of an implicit solution method for the fractional diffusion equation.  Journal of Computational Physics, 205(2):719–736, 2005.
  * [22] M Levandowsky, BS White, and FL Schuster.  Random movements of soil amebas.  Acta Protozoologica, 36:237–248, 1997.
  * [23] Changpin Li and Fanhai Zeng.  Numerical methods for fractional calculus, volume 24.  CRC Press, 2015.
  * [24] Xianjuan Li and Chuanju Xu.  A space-time spectral method for the time fractional diffusion equation.  SIAM Journal on Numerical Analysis, 47(3):2108–2131, 2009.
  * [25] Hong-lin Liao, William McLean, and Jiwei Zhang.  A discrete Gronwall inequality with applications to numerical schemes for subdiffusion problems.  SIAM Journal on Numerical Analysis, 57(1):218–237, 2019.
  * [26] Yumin Lin and Chuanju Xu.  Finite difference/spectral approximations for the time-fractional diffusion equation.  Journal of computational physics, 225(2):1533–1552, 2007.
  * [27] Peter Linz.  Analytical and numerical methods for Volterra equations.  SIAM, 1985.
  * [28] DA Lockington and JY Parlange.  Anomalous water absorption in porous materials.  Journal of Physics D: Applied Physics, 36(6):760, 2003.
  * [29] Kassem Mustapha.  FEM for time-fractional diffusion equations, novel optimal error analyses.  Mathematics of Computation, 87(313):2259–2272, 2018.
  * [30] Kassem Mustapha and Dominik Schötzau.  Well-posedness of hp-version discontinuous Galerkin methods for fractional diffusion wave equations.  IMA Journal of Numerical Analysis, 34(4):1426–1446, 2014.
  * [31] Hanna Okrasińska-Płociniczak and Łukasz Płociniczak.  Second order scheme for self-similar solutions of a time-fractional porous medium equation on the half-line.  Applied Mathematics and Computation, 424:127033, 2022.
  * [32] W Okrasiński and S Vila.  Power series solutions to some nonlinear diffusion problems.  Zeitschrift für angewandte Mathematik und Physik ZAMP, 44(6):988–997, 1993.
  * [33] Gianni Pagnini.  Erdélyi-Kober fractional diffusion.  Fractional calculus and applied analysis, 15(1):117–127, 2012.
  * [34] Łukasz Płociniczak.  Approximation of the Erdélyi–Kober operator with application to the time-fractional porous medium equation.  SIAM journal on applied mathematics, 74(4):1219–1237, 2014.
  * [35] Łukasz Płociniczak.  Analytical studies of a time-fractional porous medium equation. derivation, approximation and applications.  Communications in Nonlinear Science and Numerical Simulation, 24(1-3):169–183, 2015.
  * [36] Łukasz Płociniczak.  Numerical method for the time-fractional porous medium equation.  SIAM journal on numerical analysis, 57(2):638–656, 2019.
  * [37] Łukasz Płociniczak.  Error of the Galerkin scheme for a semilinear subdiffusion equation with time-dependent coefficients and nonsmooth data.  Computers & Mathematics with Applications, 27:181–191, 2022.
  * [38] Łukasz Płociniczak.  Linear Galerkin-Legendre spectral scheme for a degenerate nonlinear and nonlocal parabolic equation arising in climatology.  Applied Numerical Mathematics, 179:105–124, 2022.
  * [39] Łukasz Płociniczak.  A linear Galerkin numerical method for a quasilinear subdiffusion equation.  Applied Numerical Mathematics, 185:203–220, 2022.
  * [40] Łukasz Płociniczak and Szymon Sobieszek.  Numerical schemes for integro-differential equations with Erdélyi-Kober fractional operator.  Numerical Algorithms, 76(1):125–150, 2017.
  * [41] Łukasz Płociniczak and Mateusz Świtała.  Existence and uniqueness results for a time-fractional nonlinear diffusion equation.  Journal of mathematical analysis and applications, 462(2):1425–1434, 2018.
  * [42] Fangzhou Ren, Chunsheng Zhou, Qiang Zeng, Zhidong Zhang, Ueli Angst, and Wei Wang.  Quantifying the anomalous water absorption behavior of cement mortar in view of its physical sensitivity to water.  Cement and Concrete Research, 143:106395, 2021.
  * [43] Achim Schädle, María López-Fernández, and Christian Lubich.  Fast and oblivious convolution quadrature.  SIAM Journal on Scientific Computing, 28(2):421–438, 2006.
  * [44] Stefan Schaufler, WP Schleich, and VP Yakovlev.  Scaling and asymptotic laws in subrecoil laser cooling.  EPL (Europhysics Letters), 39(4):383, 1997.
  * [45] Ian Naismith Sneddon.  The use in mathematical physics of Erdélyi-Kober operators and of some of their generalizations.  Fractional Calculus and its applications, pages 37–79, 1975.
  * [46] Martin Stynes, Eugene O’Riordan, and José Luis Gracia.  Error analysis of a finite difference method on graded meshes for a time-fractional diffusion equation.  SIAM Journal on Numerical Analysis, 55(2):1057–1079, 2017.
  * [47] HongGuang Sun, Mark M Meerschaert, Yong Zhang, Jianting Zhu, and Wen Chen.  A fractal Richards equation to capture the non-Boltzmann scaling of water transport in unsaturated media.  Advances in Water Resources, 52:292–295, 2013.
  * [48] Titiwat Sungkaworn, Marie-Lise Jobin, Krzysztof Burnecki, Aleksander Weron, Martin J Lohse, and Davide Calebiro.  Single-molecule imaging reveals receptor–G protein interactions at cell surface hot spots.  Nature, 550(7677):543, 2017.
  * [49] Juan Luis Vázquez.  An introduction to the mathematical theory of the porous medium equation.  In Shape optimization and free boundaries, pages 347–389. Springer, 1992.
  * [50] Juan Luis Vázquez.  The porous medium equation: mathematical theory.  Oxford University Press on Demand, 2007.
  * [51] Licheng Wang and Tamon Ueda.  Mesoscale modeling of water penetration into concrete by capillary absorption.  Ocean Engineering, 38(4):519–528, 2011.
  * [52] Petra Wittbold, Patryk Wolejko, and Rico Zacher.  Bounded weak solutions of time-fractional porous medium type and more general nonlinear and degenerate evolutionary integro-differential equations.  Journal of Mathematical Analysis and Applications, 499(1):125007, 2021.
  * [53] Z Wu, HS Wong, C Chen, and NR Buenfeld.  Anomalous water absorption in cement-based materials caused by drying shrinkage induced microcracks.  Cement and Concrete Research, 115:90–104, 2019.
  * [54] E Zeidler.  Nonlinear functional analysis and its applications I.(Fixed Point Theorems) 1986.
  * [55] Yixin Zhao, Yang Wu, Chuanlong Dong, Songbai Han, Derek Elsworth, and Linfeng He.  Hydraulic characterization and modeling of water diffusivity through direct neutron radiography measurement on unsaturated cracked sandstone.  International Journal of Heat and Mass Transfer, 196:123256, 2022\.