# Second order scheme for self-similar solutions of a time-fractional porous
medium equation on the half-line111This is an accepted version of the
manuscript published in Applied Mathematics and Computation 424C (2022),
127033 with DOI: https://doi.org/10.1016/j.amc.2022.127033

Hanna Okrasińska-Płociniczak , Łukasz Płociniczak$\;{}^{,}$222Corresponding
author. Mail<EMAIL_ADDRESS>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, Wyb. Wyspiańskiego 27, 50-370 Wrocław, Poland

###### Abstract

Many physical, biological, and economical systems exhibit various memory
effects due to which their present state depends on the history of the whole
evolution. Combined with the nonlinearity of the process these phenomena pose
serious difficulties in both analytical and numerical treatment. We
investigate a time-fractional porous medium equation that has proved to be
important in many applications, notably in hydrology and material sciences. We
show that solutions of the free boundary Dirichlet, Neumann, and Robin
problems on the half-line satisfy a Volterra integral equation with a non-
Lipschitz nonlinearity. Based on this result we prove existence, uniqueness,
and construct a family of numerical methods that solve these equations
outperforming the usual finite difference approach. Moreover, we prove the
convergence of these methods and support the theory with several numerical
examples.

Keywords: numerical method, porous medium equation, anomalous diffusion,
fractional derivative, Volterra equation, non-Lipschitz

AMS Classification: 35R11, 45G10, 65R20

##  1 Introduction

During the last few decades experimentalists, modellers, theoreticians, and
numerical analysts have been investigating many aspects of nonlocality in
mathematical objects. These inquiries have certainly been motivated by an
increasingly rich experimental evidence for anomalous behaviour,
mathematically nontrivial models, and professional curiosity. Apart from many
applications, probably the one that is most researched is the so-called
anomalous diffusion that arises when the transport is slower or faster than in
the classical case [1, 2]. Probabilistically, this can be characterised with a
deviation from the linear mean-square displacement for a randomly walking
particle [3]. There is a broad experimental evidence of such dynamics in
quantum optics [4], physics of plasma [5], movement of bacteria and amoeba
[6], G-protein on cell surface [7], hydrology [8, 9, 10, 11, 12], and motion
of dislocations in crystal lattice [13].

The nonlocality in the anomalous diffusion in usually modelled with a use of
the fractional derivative [1]. There are many definitions of such operators
and it depends on the situation which one is the most suited for a specific
model. For example, temporal nonlocality (memory) can be described with
Riemann-Liouville or Caputo derivatives while nonlocality in space requires
,,symmetric” operator such as the fractional Laplacian [14, 15, 16]. Physical
derivations of the fractional porous media equation have been given for
example in [17, 18]. We have argued there that the emergence of the time-
fractional derivative is a manifestation of the waiting-time phenomenon
happening in the medium, i.e. the fluid can be trapped in certain regions for
prolonged periods of time which leads to the so-called subdiffusion. Recently,
some experimental evidence has been given in [19] that supports this
hypothesis. On the other hand, the spatial nonlocality in the porous medium
equation is a consequence of the long jump phenomenon (for a physical
derivation see [18]). That is to say, due to medium’s complex geometry and
heterogeneity fluid parcels can traverse relatively long distances in short
amounts of time what in the literature is called superdiffusion.
Mathematically this leads to the fractional gradient operator which is closely
related to the fractional Laplacian [15, 20, 21]. A thorough summary of the
theory of the fractional porous medium equation can be found [16].

Our previous studies investigated the time-fractional porous medium equation
applied as a model for moisture transport in construction materials

$\partial^{\alpha}_{t}u=\left(u^{m}u_{x}\right)_{x},\quad 0<\alpha<1,\quad
m>1,$ (1)

where the fractional derivative is of the Riemann-Liouville type

$\partial^{\alpha}_{t}u(x,t)=\frac{1}{\Gamma(1-\alpha)}\frac{\partial}{\partial
t}\int_{0}^{t}(t-s)^{-\alpha}u(x,s)ds.$ (2)

with the vanishing initial condition and dry medium at infinity

$u(x,0)=0,\quad u(\infty,t)=0,\quad x>0,\quad t>0.$ (3)

The natural function space in which the weak solution of a fractional
diffusion equation should be sought is

$C\left([0,1];H^{1}(\mathbb{R}_{+})\right),$ (4)

where $H^{1}$ is the usual Sobolev space. Below, by the use of integral
equation techniques, we show that the self-similar solution of the above
problem is indeed continuous and differentiable almost everywhere with respect
to space and time. Moreover, in physical derivation of the time-fractional
porous medium equation (1) one arrives at the Caputo derivative on the left-
hand side (see for ex. [17]). However, due to the zero initial condition this
derivative coincides with the Riemann-Liouville one [14] and hence the present
form of (1). Additionally, motivated by the modelling using continuous time
random walk it is interesting to consider some general fractional operators
that govern the evolution in time (for some interesting examples and analyses
see [22, 23]). We leave this topic for future studies.

The parameters $\alpha$ and $m$ are associated with a specific experiment and
are subject to fitting with the data. Based on [12] we can think that usually
$\alpha\in[0.6,1]$ while $m\approx 7-8$. The parameter $\alpha$ denotes the
level of nonlocality of the model with $\alpha=1$ being the classical case.
Since, as will be seen below, our model (1) describes a slower evolution than
in the classical case we speak about subdiffusion.

In previous works we have proved existence and uniqueness of a self-similar
solution to the Dirichlet problem [24], provided some approximations [10],
solved inverse problem of diffusivity identification [25], and constructed a
numerical scheme [26]. This numerical method, although only first order
accurate, provided a superior performance over the finite difference method
applied to the governing PDE. The idea was to use a series of transformations
leading to a nonlinear Volterra equation and then to discretize it. In this
work we follow a similar route but now we introduce several improvements.
First of all, we sketch how to construct an arbitrarily high order scheme and
explicitly devise the second order accurate linear method. Then, we prove its
convergence under some mild assumptions. The difference from our previous work
is that we use the general theory of Volterra equations in order to anticipate
the behaviour of its solution. This lets us to conduct yet another
transformation which facilitates the crucial step in the convergence proof
yielding a stronger result. Moreover, our method is constructed in a way that
it gives exact results for power functions which are the leading-order
approximations to the solution of the time-fractional porous medium equation.
This further improves the accuracy within the second order. For example, as
will be shown below, computing the wetting front to several decimal places
requires using only moderate step sizes.

There is a vast literature concerning numerical methods for fractional
differential equations. For general surveys one can refer to [27, 28]. The
linear fractional diffusion either with temporal or spatial nonlocality has
been solved for example in [29, 30, 31]. The usual approach to constructing a
numerical scheme for the time-fractional diffusion is based on discretizing
nonlocal derivative by, for example, L1 scheme [32, 33], Grünwald-Letnikov
[31], or Convolution quadrature [34]. This provides the quasi-discrete scheme
that can be fully discretized by introducing a finite difference [35], finite
element [36] or other spatial discretization [30].

On the other hand, the case of nonlinear anomalous diffusion has been
rigorously treated only in relatively few papers (especially when the
diffusivity is a function of the solution). Some recent examples include time-
fractional reaction-diffusion equations [37]. One of the very useful tools in
proving convergence for a nonlinear problem is the discrete fractional
Grönwall inequality that was developed and applied to nonlinear time-
fractional diffusion problems in [38, 39] (note that probably the first
account of fractional Grönwall inequality has been given in [40] but see also
[41] and references therein). Numerical methods for equations with non-smooth
data have been also analysed in [42]. On the other hand, the space-fractional
porous medium has been discretized in [43] with the use of finite difference
method applied to the extension problem [44]. Furthermore, in [45] some very
robust numerical methods have been devised, thoroughly analysed, and applied
to a very broad classes of nonlocal equations. Moreover, authors of [46] used
a multi-grid waveform method for the time-fractional case. The nonlinear
source terms have been considered in [47]. The main difficulty in the analysis
of the porous medium equation (even in the classical case) is the degeneracy
of diffusivity, i.e. the coefficient of $u_{x}$ can vanish. This fact makes
both theoretical and numerical studies interesting and non-trivial. For
several years many approaches of numerical treatment of the classical porous
medium equation have been proposed. To paint the background we mention only a
few papers. For example, in [48] authors used a finite difference method with
an interface tracking algorithm. This feature is crucial since due to
degeneracy in many circumstances the solution is characterized with a finite
speed of propagation what translates into a free-boundary problem [49]. Also,
a finite element method for the classical porous medium equation has been
analysed in [50, 51]. This topic is still vigorously investigated and
mathematicians constantly develop some new approaches (for ex. [52, 53, 54]).

The root of our method lies in appropriate discretization of a nonlinear
Volterra equation (for a thorough treatment see [55]). A survey of numerical
methods for such problems can be found in [56, 57]. Notice that the majority
of the theory concerns Lipschitz nonlinearity while our model fails to satisfy
this condition. There are only a few papers that developed rigorous methods
for nonlipschitzian case [58]. Notably, in [59] an iterative scheme has been
proposed. We also mention our previous work [60] where an idea for a proof
similar to the one used in the present paper originated. Finally, in [61, 62]
a highly readable theoretical summaries about non-Lipschitz problems have been
published.

As we mentioned before, this paper concerns a high order numerical method for
finding self-similar solutions of (1). In Section 2 we perform a series of
transformations that reduce the main nonlocal nonlinear PDE into an ordinary
Volterra equation with non-Lipschitz nonlinearity. As a side result we prove a
existence and uniqueness result about self-similar solutions on the half-line
with Dirichlet, Neumann, and Robin conditions. Next, in Section 3 we present a
general scheme of numerical method construction along with explicit formulas
for a second order linear interpolation scheme (trapezoidal). We also prove
that it is convergent provided that the parameter $m$ is sufficiently large.
In Section 4 we illustrate our theory with numerical examples concerning some
exemplary equations and the time-fractional porous medium equation itself. We
close the paper with some concluding remarks.

##  2 Self-similar solutions and Volterra equations

In applications one is frequently interested in finding self-similar solutions
of (1) in the form

$u(x,t)=t^{a}U(\eta),\quad\eta=xt^{-b},$ (5)

where $a,b\geq 0$ are constants to be found and the similarity variable is
denoted by $\eta$. Our main motivation for looking for such solutions lies in
the fact that in the real-world experiments of moisture percolation in porous
media one observes the self-similar profile (see [11, 19, 9, 12]). On the
other hand, in a large generality Lie group methods can be used to look for
many families of self-similar solutions in a systematic way [63, 64, 65]. When
we plug the above into the original PDE (1) we transform the spatial

$\left(u^{m}(x,t)u_{x}(x,t)\right)_{x}=t^{a(m+1)-2b}\frac{d}{d\eta}\left(U^{m}(\eta)\frac{d}{d\eta}U(\eta)\right).$
(6)

and the temporal derivatives

$\begin{split}\partial^{\alpha}_{t}u(x,t)&=\frac{\partial}{\partial
t}\left(I^{1-\alpha}_{t}\left(t^{a}U(xt^{-b})\right)\right)=\frac{1}{\Gamma(1-\alpha)}\frac{\partial}{\partial
t}\int_{0}^{t}(t-z)^{-\alpha}z^{a}U(xz^{-b})dz\\\
&=\frac{1}{\Gamma(1-\alpha)}\frac{\partial}{\partial
t}\left(t^{a-\alpha+1}\int_{0}^{1}(1-s)^{-\alpha}s^{a}U(\eta
s^{-b})ds\right),\end{split}$ (7)

If we further introduce the Erdélyi-Kober fractional integral operator (see
[66])

$I^{\beta,\gamma}_{\delta}U(\eta)=\frac{1}{\Gamma(\gamma)}\int_{0}^{1}(1-z)^{\gamma-1}z^{\beta}U(\eta
z^{\frac{1}{\delta}})dz,$ (8)

we can write (1) as

$t^{\alpha+am-2b}\frac{d}{d\eta}\left(U^{m}\frac{dU}{d\eta}\right)=\left[(1-\alpha+a)-b\eta\frac{d}{d\eta}\right]I^{a,1-\alpha}_{-\frac{1}{b}}U,\quad
0<\alpha<1,$ (9)

Now, the self-similar solution is possible only if

$2b-ma=\alpha,$ (10)

which is the consistency condition. Since there are two unknowns:
$a=a(\alpha,m)$ and $b=b(\alpha,m)$ we need yet another equation. This comes
from the initial-boundary conditions and we assume that initially the medium
is dry (3) which expressed in the self-similar variables (5) yields

$U(\infty)=0.$ (11)

Furthermore, for the boundary conditions we have several self-similar
physically relevant choices. For example, we can consider the Dirichlet

$u(0,t)=1,\quad t>0,$ (12)

which describes constant concentration, Neumann

$u^{m}(0,t)u_{x}(0,t)=-1,\quad t>0,$ (13)

describing constant flux, or Robin problem

$u^{m}(0,t)u_{x}(0,t)=-u(0,t),\quad t>0,$ (14)

which translates into requirement that the flux is proportional to the
concentration at the interface $x=0$. For a physical account of these models
see [17]. Boundary conditions provide us with the second equation for $a$ and
$b$ and eventually we are able to write (1) as

$\frac{d}{d\eta}\left(U^{m}\frac{dU}{d\eta}\right)=AI^{a,1-\alpha}_{-\frac{1}{b}}U-B\eta\frac{d}{d\eta}I^{a,1-\alpha}_{-\frac{1}{b}}U,\quad
0<\alpha<1,$ (15)

for some $A=A(\alpha,m)$ and $B=B(\alpha,m)$. The exact values of these
constants for different boundary conditions are summarized in Tab. 1. Note
that always $A,B\geq 0$.

Condition | $a$ | $A$ | $b=B$ | $\eta^{*}$
---|---|---|---|---
Dirichlet (12) | 0 | $1-\alpha$ | $\dfrac{\alpha}{2}$ | $\dfrac{1}{\sqrt{y(1)^{m}}}$
Neumann (13) | $\dfrac{\alpha}{m+2}$ | $1-\dfrac{m+1}{m+2}$ | $\dfrac{m+1}{m+2}\alpha$ | $\left(\dfrac{m+1}{(y^{m+1})^{\prime}(1)}\right)^{\frac{m}{m+2}}$
Robin (14) | $\dfrac{\alpha}{m}$ | $1-\dfrac{m-1}{m}\alpha$ | $\alpha$ | $\dfrac{m}{(y^{m})^{\prime}(1)}$
Table 1: Values of $a$, $b$ in (5), $A=A(\alpha,m)$, $B=B(\alpha,m)$ as in
(15), and the wetting front $\eta^{*}$ defined in (17) for different types of
boundary conditions.

To motivate our further reasoning let us focus on the Robin problem (14) for
which

$a=\frac{\alpha}{m},\quad b=\alpha.$ (16)

Other conditions have been considered in [17, 24] and lead to a similar final
integral equation. In any of the considered boundary conditions the support of
the solution will be compact (for a proof see [67] and for the classical case
[68]). Physically, this fact is equivalent to a finite speed of wetting front
propagation and arises in the degeneracy of the equation (1). Having that in
mind, let $\eta^{*}>0$ be such that

$U(\eta)=0\quad\text{for}\quad\eta>\eta^{*}.$ (17)

Moreover, by integrating (15) one can show that there is no flux through the
support’s boundary (see [24])

$U(\eta)^{m}\frac{dU}{d\eta}(\eta)=0\quad\text{for}\quad\eta>\eta^{*}.$ (18)

Note that $\eta^{*}$ is a-priori unknown and has to be determined as a part of
the solution since we are dealing with a free-boundary problem. A standard way
of proceeding is to introduce yet another transformation

$U(\eta)=Cy(z),\quad z=1-\frac{\eta}{\eta^{*}},\quad z\geq 0,$ (19)

where $C$ has to be determined. The conditions (17)-(18) now become

$y(0)=0,\quad y^{m}(0)y^{\prime}(0)=0,$ (20)

where prime denotes the differentiation with respect to $z$. The Robin
condition (14) now yields

$U^{m}(0)\frac{dU}{d\eta}(0)=-U(0)\rightarrow(y^{m})^{\prime}(1)=\frac{m\eta^{*}}{C^{m}},$
(21)

which gives us one equation to determine $\eta^{*}$ and $C$.

When we plug (19) into (15) we obtain

$\frac{C^{m}}{(\eta^{*})^{2}}\left(y^{m}y^{\prime}\right)^{\prime}=AG_{\alpha}y+B(1-z)(G_{\alpha}y)^{\prime},$
(22)

where for $a$ and $b$ as in (16) we have defined

$\begin{split}G_{\alpha}y(z)&=\frac{1}{\Gamma(1-\alpha)}\int_{(1-z)^{1/b}}^{1}(1-s)^{-\alpha}s^{a}y(1-s^{-b}(1-z))ds\\\
&=\frac{1}{b}\frac{1}{\Gamma(1-\alpha)}\left(1-z\right)^{\frac{a+1}{b}}\int_{0}^{z}\left(1-\left(\frac{1-z}{1-w}\right)^{\frac{1}{b}}\right)^{-\alpha}\frac{y(w)}{(1-w)^{1+\frac{a+1}{b}}}dw.\end{split}$
(23)

The second equality above follows from the change of variables
$w=1-s^{-a}(1-z)$. We can thus see that

$C=\left(\eta^{*}\right)^{\frac{2}{m}},$ (24)

which together with (21) closes the system and gives us an equation for
$\eta^{*}$ provided we know the solution $y$, that is $\eta^{*}$ can be found
from

$\eta^{*}=\frac{m}{(y^{m})^{\prime}(1)}.$ (25)

Now, (22) can be transformed into a nonlinear integral Volterra equation.
First, with an integration and using the condition (20) we obtain

$\begin{split}y(z)^{m}y^{\prime}(z)&=A\int_{0}^{z}G_{\alpha}y(s)ds+B\int_{0}^{z}\left(1-z)(G_{\alpha}y(s)\right)^{\prime}ds\\\
&=B(1-z)G_{\alpha}y(z)+(A+B)\int_{0}^{z}G_{\alpha}y(s)ds,\end{split}$ (26)

where we have integrated by parts. A second integration finally yields

$\begin{split}y(z)^{m+1}&=(m+1)\int_{0}^{z}\left(B(1-s)+(A+B)\left(z-s\right)\right)G_{\alpha}y(s)ds\\\
&=:\int_{0}^{z}F(z,s)G_{\alpha}y(s)ds.\end{split}$ (27)

Next, we use the definition of $G_{\alpha}$ operator (23) and change the order
of integration (Fubini’s theorem)

$y(z)^{m+1}=\frac{1}{b}\frac{1}{\Gamma(1-\alpha)}\int_{0}^{z}\left(\int_{u}^{z}F(z,s)(1-s)^{\frac{a+1}{b}}\left(1-\left(\frac{1-s}{1-u}\right)^{\frac{1}{b}}\right)^{-\alpha}ds\right)\frac{y(u)du}{(1-u)^{1+\frac{a+1}{b}}}.$
(28)

We can now substitute back $v=((1-s)/(1-u))^{1/b}$ inside the inner integral
to obtain

$\begin{split}y(z)^{m+1}&=\frac{1}{\Gamma(1-\alpha)}\int_{0}^{z}\left[\int_{\left(\frac{1-z}{1-u}\right)^{1/b}}^{1}F(z,1-s^{b}(1-u))(1-s)^{-\alpha}s^{a+b}ds\right]y(u)du\\\
&=:\int_{0}^{z}K(z,u)y(u)du,\end{split}$ (29)

where we have defined the kernel $K$. Observe that it can be written as a
combination of incomplete beta functions (and we use this fact later).
Moreover, the kernel is a continuous function since $(1-s)^{-\alpha}$ is
integrable near $s=1$ and $a>0$. We can also obtain some simple bounds on $K$
due to the fact that $F(z,s)$ is a linear function in its both variables.
Assume that $0\leq z\leq X$ for some $0<X<1$. Then,

$B(1-X)\leq F(z,s)\leq A+2B,$ (30)

and hence

$K(z,u)\leq\frac{m+1}{\Gamma(1-\alpha)}(A+2B)\int_{\left(\frac{1-z}{1-u}\right)^{1/b}}^{1}(1-s)^{-\alpha}ds=\frac{m+1}{\Gamma(2-\alpha)}(A+2B)\left(1-\left(\frac{1-z}{1-u}\right)^{1/b}\right)^{1-\alpha}.$
(31)

Now, we have

$\frac{1-z}{1-s}=1-\frac{z-s}{1-s},$ (32)

what together with elementary estimate $(1-x)^{\gamma}\geq 1-\gamma x$ for
$\gamma\geq 1$ leads to

$K(z,u)\leq\frac{m+1}{\Gamma(2-\alpha)}(A+2B)\left(\frac{1}{b}\frac{z-u}{1-X}\right)^{1-\alpha}\leq
K_{+}(z-u)^{1-\alpha},$ (33)

because $0\leq u\leq z\leq X$. Similarly, we can obtain bounds from below. To
this end, notice that

$\begin{split}K(z,u)&\geq\frac{m+1}{\Gamma(1-\alpha)}B(1-X)\left(\frac{1-z}{1-u}\right)^{a/b+1}\int_{\left(\frac{1-z}{1-u}\right)^{1/b}}^{1}(1-s)^{-\alpha}ds\\\
&=\frac{m+1}{\Gamma(2-\alpha)}B(1-X)\left(\frac{1-z}{1-u}\right)^{a/b+1}\left(1-\left(\frac{1-z}{1-u}\right)^{1/b}\right)^{1-\alpha}.\end{split}$
(34)

Using the convexity bound $(1-x)^{\gamma}\leq 1-x$ for $\gamma\geq 1$ we now
have

$K(z,u)\geq\frac{m+1}{\Gamma(2-\alpha)}B(1-X)\left(\frac{1-z}{1-u}\right)^{a/b+1}\left(\frac{z-u}{1-u}\right)^{1-\alpha}\geq
K_{-}(z-u)^{1-\alpha},$ (35)

since $0\leq u\leq z\leq X$. Whence, we have found the behaviour of the kernel
yielding

$K_{-}(z-u)^{1-\alpha}\leq K(z,u)\leq K_{+}(z-u)^{1-\alpha},$ (36)

where $K_{\pm}$ are known constants. These kernel bounds have important
consequence for the form of the solution and are crucial for our subsequent
construction of the numerical method. For example, as was shown in [69] (for a
further account see also [70, 71, 72]) the equation (29) has only one non-
trivial solution which satisfies

$C_{-}^{\frac{1}{m}}z^{\frac{2-\alpha}{m}}\leq y(z)\leq
C_{+}^{\frac{1}{m}}z^{\frac{2-\alpha}{m}},$ (37)

where $C_{\pm}>0$ can be found explicitly, however, their exact values are not
relevant for our reasoning. The calculations for other boundary conditions
presented in Tab. 1 are essentially identical. Therefore, we have shown the
following result.

######  Theorem 1.

Let

$u(x,t)=(\eta^{*})^{\frac{2}{m}}t^{a}y\left(1-\frac{1}{\eta^{*}}\frac{x}{t^{b}}\right),$
(38)

where $y=y(z)$ is a unique nontrivial solution of the Volterra equation (29)
while $a$, $b$, and $\eta^{*}$ are given in Tab. 1. Then, $u$ is a unique
self-similar solution of (1) with vanishing initial condition along with one
of the self-similar boundary conditions (12)-(14). Moreover, we have the
following estimates

$U_{-}\left(1-\frac{1}{\eta^{*}}\frac{x}{t^{b}}\right)^{\frac{2-\alpha}{m}}\leq
u(x,t)\leq
U_{+}\left(1-\frac{1}{\eta^{*}}\frac{x}{t^{b}}\right)^{\frac{2-\alpha}{m}},$
(39)

with suitable constants $U_{\pm}>0$ that can be found explicitly.

For example, when we consider the Dirichlet boundary condition and consult
Tab. 1 we obtain

$K_{-}=\frac{m+1}{\Gamma(2-\alpha)}\frac{\alpha}{2}(1-X)^{\alpha-1},\quad
K_{+}=\frac{m+1}{\Gamma(2-\alpha)}\left(\frac{\alpha}{2}(1-X)\right)^{\alpha-1},$
(40)

which by results from [69] yield

$C_{\pm}=K_{\pm}\beta\left(2-\alpha,1+\frac{2-\alpha}{m}\right),$ (41)

where $\beta$ is the Euler beta function. Eventually, the solution of (1) can
be bounded by

$U_{\pm}=\frac{1}{y(1)}C_{\pm}^{\frac{1}{m}},$ (42)

where $y(1)$ can be calculated from the solution of (29).

##  3 Numerical method

In the previous section we have shown that looking for solutions of (1) with
self-similar boundary conditions (12)-(14) is equivalent to solving (29).
Therefore, in what follows we will devise a numerical method for solving a
general class of nonlinear Volterra equations of the form

$y(z)^{m+1}=\int_{0}^{z}K(z,s)y(s)ds,\quad 0\leq z\leq 1,$ (43)

where we assume that the kernel is continuous and there exist constants
$K_{\pm}$ such that

$K_{-}(z-s)^{\gamma}\leq K(z,s)\leq K_{+}(z-s)^{\gamma},\quad\gamma\geq 0.$
(44)

Note that in this section we use the same letters for $y$ and $K$ as before,
however now we are considering a _general_ Volterra equation which may not
have anything in common with anomalous diffusion. By results from [69] the
above equation has a unique non-trivial solution satisfying

$C_{-}^{\frac{1}{m}}z^{\frac{\gamma+1}{m}}\leq y(z)\leq
C_{+}^{\frac{1}{m}}z^{\frac{\gamma+1}{m}},$ (45)

for suitable constants $C_{\pm}>0$. Therefore, the behaviour of the solution
is of power type and for a numerical treatment it is reasonable to peel it off
from the actual form of $y$. That is to say, we substitute

$y(z)=z^{\frac{\gamma+1}{m}}v(z),$ (46)

which leads to an equivalent integral equation

$v(z)=z^{-\frac{(m+1)(\gamma+1)}{m}}\int_{0}^{z}K(z,s)s^{\frac{\gamma+1}{m}}v(s)ds.$
(47)

It is crucial to note that $v$ is now bounded away from zero, that is

$0<C_{-}^{\frac{1}{m}}\leq v(z)\leq C_{+}^{\frac{1}{m}},$ (48)

what facilitates both the analysis and numerical computations. A particular
choice of the kernel, i.e. when $K(z,s)=K_{+}(z-s)^{\gamma}$ clearly leads to
a constant solution

$v(z)=\left(K_{+}\beta\left(\gamma+1,\frac{\gamma+1}{m}+1\right)\right)^{\frac{1}{m}},$
(49)

where $\beta$ is Euler beta function. This simple solution can serve as a
benchmark of various numerical methods.

###  3.1 Construction

Now we can proceed to the discretization. Introduce the grid

$z_{n}:=\frac{n}{N},\quad h:=\frac{1}{N},\quad n=0,1,...,N,$ (50)

where the total number of grid points $N$ is fixed. Further, we can discretize
the integral in (47)

$\int_{0}^{z_{n}}K(z_{n},s)s^{\frac{\gamma+1}{m}}v(s)ds=\sum_{i=0}^{n-1}\left(w_{n,i}(h)v(z_{i})+\delta_{i}(h)\right),$
(51)

where $\delta_{i}(h)$ is the local consistency error, and $w_{n,i}(h)$ are
weights that depend on the form of the kernel. When we use (51) in (47) and
truncate the remainder we obtain the following numerical scheme

$v_{n}^{m+1}=z_{n}^{-\frac{(m+1)(\gamma+1)}{m}}\sum_{i=0}^{n-1}w_{n,i}(h)v_{i},$
(52)

where by $v_{i}$ we have denoted the numerical approximation to $v(z_{i})$
while $K_{n,i}=K(z_{n},z_{i})$. Note also that in the quadrature (51) we have
not included the last point of the interval $[0,z_{n}]$. This leads to a
(half-)open quadrature and has been chosen in order to reduce the
computational cost of the method. For if the sum on the right-hand side of the
above included $v_{n}$ term we would obtain an implicit method that would
require solving a nonlinear equation of the form $x^{m+1}-a_{1}x+a_{0}=0$ in
each step. In order to keep the cost as low as possible we consider only
explicit methods.

Different quadratures would yield different values of the weights
$w_{n,i}(h)$. For example, the simplest rectangular rule in which the whole
integrand is approximated by its value at the lower terminal would yield

$w_{n,i}(h)=hK(z_{n},z_{i})z_{i}^{\frac{\gamma+1}{m}}.$ (53)

This method could be proved to be convergent (see [26]), however, it does not
solve the constant case exactly, that is to say when
$K(z,s)=K_{+}(z-s)^{\gamma}$ the numerical solution $v_{n}$ is not equal to
(49). To see this suppose that $v_{0}=v_{1}=v_{2}=C$, then from (52) at $n=2$
and (53)

$C^{m+1}=CK_{+}h^{-\frac{(m+1)(\gamma+1)}{m}+1}(2h-h)^{\gamma}h^{\frac{\gamma+1}{m}}=CK_{+},$
(54)

hence $C=(K_{+})^{1/m}$. Then, the third step yields

$v_{3}^{m+1}=K_{+}h^{-\frac{(m+1)(\gamma+1)}{m}+1}\left((3h-h)^{\gamma}h^{\frac{\gamma+1}{m}}+(3h-2h)^{\gamma}(2h)^{\frac{\gamma+1}{m}}\right)=K_{+}\left(2^{\gamma}+1\right)\neq
C^{m+1},$ (55)

for $\gamma>0$. Therefore, $v_{n}$ cannot be constant for all $n\geq 3$. An
ability of a numerical scheme to be able to solve for a constant solution can
be thought as a necessary requirement that we have to make since it would
accurately resolve the zero order Taylor series term. A family of such methods
can be devised by using an interpolating polynomial for the unknown $v_{n}$ in
the integral (51). The kernel, since it is known, is not approximated.
Although in this work we will focus only on the first degree interpolation,
i.e. a linear reconstruction, it is instructive to see how does the zeroth
order approximation look like. To this end we write

$\int_{0}^{z_{n}}K(z_{n},s)s^{\frac{\gamma+1}{m}}v(s)ds=\sum_{i=0}^{n-1}\int_{z_{i}}^{z_{i+1}}K(z_{n},s)s^{\frac{\gamma+1}{m}}v(s)ds,$
(56)

in which we expand $v$ in its Taylor series
$v(s)=v(z_{i})+v^{\prime}(\widehat{\xi}_{i})(s-z_{i})$, which leads to

$\int_{0}^{z_{n}}K(z_{n},s)s^{\frac{\gamma+1}{m}}v(s)ds=\sum_{i=0}^{n-1}\left[\left(\int_{z_{i}}^{z_{i+1}}K(z_{n},s)s^{\frac{\gamma+1}{m}}ds\right)v_{i}+\delta_{i}(h)\right],$
(57)

where by the mean value theorem the remainder is

$\delta_{i}(h)=v^{\prime}(\xi_{i})\int_{z_{i}}^{z_{i+1}}K(z_{n},s)s^{\frac{\gamma+1}{m}}(s-z_{i})ds.$
(58)

Therefore, going back to (47) gives

$v(z_{n})^{m+1}=z_{n}^{-\frac{(m+1)(\gamma+1)}{m}}\sum_{i=0}^{n-1}\left(\int_{z_{i}}^{z_{i+1}}K(z_{n},s)s^{\frac{\gamma+1}{m}}ds\right)v_{i}+\delta_{n}(h),$
(59)

where the remainder satisfies

$\begin{split}|\delta_{n}(h)|\leq
z^{-\frac{(m+1)(\gamma+1)}{m}}\sum_{i=0}^{n-1}|\delta_{i}(h)|&=z_{n}^{-\frac{(m+1)(\gamma+1)}{m}}|v^{\prime}(\xi)|\sum_{i=0}^{n-1}\int_{z_{i}}^{z_{i+1}}K(z_{n},s)s^{\frac{\gamma+1}{m}}(s-z_{i})ds\\\
&\leq
z_{n}^{-\frac{(m+1)(\gamma+1)}{m}}K_{+}|v^{\prime}(\xi)|h\sum_{i=0}^{n-1}\int_{z_{i}}^{z_{i+1}}(z_{n}-s)^{\gamma}s^{\frac{\gamma+1}{m}}ds.\end{split}$
(60)

where, once again, we have used the mean value theorem and used (36). Further,
with a substitution $s=nh\sigma$ we can evaluate the integral

$|\delta_{n}(h)|\leq
K_{+}\|v^{\prime}(\xi)\|\beta\left(\gamma+1,\frac{\gamma+1}{m}\right)h=O(h),\quad
h\rightarrow 0^{+},$ (61)

uniformly for $n\in\mathbb{N}$. Whence, the truncation error of the method is
of the first order. Here, $\beta(\cdot,\cdot)$ is the Euler beta function.
Neglecting the remainder we obtain a rectangle method for solving (47)

$v_{n}^{m+1}=z_{n}^{-\frac{(m+1)(\gamma+1)}{m}}\sum_{i=0}^{n-1}\left(\int_{z_{i}}^{z_{i+1}}K(z_{n},s)s^{\frac{\gamma+1}{m}}ds\right)v_{i}.$
(62)

In contrast with (53) the above method solves for the constant solution
exactly by the very construction (in that case $v^{\prime}\equiv 0$ and the
remainder vanishes).

It is now a straightforward exercise to develop higher order methods with the
use of the higher degree of Lagrange’s interpolating polynomial. We will
analyse only the second order method since we are able to equip it with a
suitable choice of initial conditions. As we shall see below, for higher order
methods a prescription of starting values is not a straightforward task. This
situation is similar to multistep methods for ODEs. Having that in mind we can
use a linear approximation to $v$ in the integral in (51). The important point
is to construct an explicit method by not including the terminal point of the
interval into the interpolation. That is to say, we have to partition the
interval $[0,z_{n}]$ into subintervals which contain exactly three nodes.
Then, we construct a linear interpolation based on first two nodes. Due to
this construction we have to consider separately cases when $n$ is even or
odd.

First, suppose that $n$ is even. Then, we build linear approximations based on
the first and second points of a three node interval starting with $z=0$. That
is to say, we have the partition

$\int_{0}^{z_{n}}K(z_{n},s)s^{\frac{\gamma+1}{m}}v(s)ds=\sum_{i=0}^{\frac{n}{2}-1}\int_{z_{2i}}^{z_{2i+2}}K(z_{n},s)s^{\frac{\gamma+1}{m}}v(s)ds,$
(63)

in which we use Lagrange’s interpolation

$v(s)=\left(1-\frac{s-z_{2i}}{h}\right)v(z_{2i})+\frac{s-z_{2i}}{h}v(z_{2i+1})+\frac{1}{2}v^{\prime\prime}(\widehat{\xi}_{i})(s-z_{2i})(s-z_{2i+1}).$
(64)

Here, $\widehat{\xi}_{i}$ is the intermediate point needed in the remainder.
If we use that in the integral we obtain

$\begin{split}\int_{0}^{z_{n}}K(z_{n},s)s^{\frac{\gamma+1}{m}}v(s)ds&=\sum_{i=0}^{\frac{n}{2}-1}w_{n,2i}(h)v(z_{2i})+w_{n,2i+1}(h)v(z_{2i})\\\
&+\frac{1}{2}v^{\prime\prime}(\xi)\sum_{i=0}^{\frac{n}{2}-1}\int_{z_{2i}}^{z_{2i+2}}K(z_{n},s)s^{\frac{\gamma+1}{m}}(s-z_{2i})(s-z_{2i+1})ds,\end{split}$
(65)

where we have utilized the mean value theorems for integrals and sums. The
weights are given by

$w_{n,i}(h)=\int_{z_{i}}^{z_{i+2}}K(z_{n},s)s^{\frac{\gamma+1}{m}}\left(\begin{cases}1-\frac{s-z_{i}}{h},&i\text{
even}\\\ \frac{s-z_{i}}{h},&i\text{ odd}\end{cases}\right)ds\quad\text{for
}n\text{ even and }0\leq i\leq n-2.$ (66)

A similar calculation can be conducted for odd values of $n$. In that case,
however, we first linearly interpolate between $z_{0}$ and $z_{1}$ and then
partition the rest of the interval $[z_{1},z_{n}]$ into three node
subintervals. We then have

$\int_{0}^{z_{n}}K(z_{n},s)s^{\frac{\gamma+1}{m}}v(s)ds=\int_{0}^{z_{1}}K(z_{n},s)s^{\frac{\gamma+1}{m}}v(s)ds+\sum_{i=1}^{\frac{n-1}{2}}\int_{z_{2i-1}}^{z_{2i+1}}K(z_{n},s)s^{\frac{\gamma+1}{m}}v(s)ds.$
(67)

Further, conducting analogous interpolation as in the even case leads to

$\begin{split}\int_{0}^{z_{n}}&K(z_{n},s)s^{\frac{\gamma+1}{m}}v(s)ds=w_{n,0}(h)v(0)+\sum_{i=1}^{\frac{n-1}{2}}w_{n,2i-1}(h)v(z_{2i-1})+w_{n,2i}(h)v(z_{2i})\\\
&+\frac{1}{2}v^{\prime\prime}(\xi_{1})\int_{0}^{z_{1}}K(z_{n},s)s^{1+\frac{\gamma+1}{m}}(s-z_{1})ds+\sum_{i=1}^{\frac{n-1}{2}}\int_{z_{2i-1}}^{z_{2i+1}}K(z_{n},s)s^{\frac{\gamma+1}{m}}(s-z_{2i-1})(s-z_{2i})v^{\prime\prime}(\xi_{i}(s))ds,\end{split}$
(68)

where $\xi_{s}$ is an intermediate point. In that case, the weights are given
by

$\begin{split}w_{n,0}(h)=\int_{0}^{z_{1}}K(z_{n},s)s^{\frac{\gamma+1}{m}}\left(1-\frac{s}{h}\right)ds,\\\
w_{n,1}(h)=\int_{0}^{z_{1}}K(z_{n},s)s^{\frac{\gamma+1}{m}}\frac{s}{h}ds+\int_{z_{1}}^{z_{3}}K(z_{n},s)s^{\frac{\gamma+1}{m}}\left(1-\frac{s-z_{1}}{h}\right)ds,\\\
w_{n,i}(h)=\int_{z_{i-1}}^{z_{i+1}}K(z_{n},s)s^{\frac{\gamma+1}{m}}\left(\begin{cases}\frac{s-z_{i-1}}{h},&i\text{
even}\\\ 1-\frac{s-z_{i-1}}{h},&i\text{ odd}\\\
\end{cases}\right)ds\quad\text{for }n\text{ odd and }1<i\leq n-1.\end{split}$
(69)

Therefore, the whole equation (47) using our linear reconstruction can be
written as

$v(z_{n})^{m+1}=z_{n}^{-\frac{(m+1)(\gamma+1)}{m}}\sum_{i=0}^{n-1}w_{n,i}(h)v(z_{i})+\delta_{n}(h),$
(70)

where the remainder can be estimated with the help of (65) and (68) to be

$|\delta_{n}(h)|\leq
K_{+}\|v^{\prime\prime}\|B\left(\gamma+1,\frac{\gamma+1}{m}\right)h^{2}=:Ch^{2}=O(h^{2})\quad\text{as}\quad
h\rightarrow 0^{+},$ (71)

uniformly for $n\in\mathbb{N}$. Therefore, by truncating the remainder we
obtain a consistent second order method in the form (52) where weights are
defined in (66) and (69).

Having designed a method we can move to the important question about the
initial conditions for $v$. Since (43) always possesses a trivial solution
$y\equiv 0$ we have to choose an appropriate starting value for our numerical
method in order to converge to a non-trivial one. One way of doing that is to
approximate the integral in (47) over the interval $z\in[0,h]$ with some
simple quadrature. For example, we can use the rectangle rule in which we take
the value of $v$ at the _right_ endpoint. That is, we reconstruct $v$ with a
constant function

$v(s)=v(h)-v^{\prime}(\xi_{h})(h-s),\quad s\in(0,h),$ (72)

for some $\xi_{h}\in(0,h)$. Then, from (47) we can write

$v(h)^{m+1}=h^{-\frac{(m+1)(\gamma+1)}{m}}v(h)\int_{0}^{h}K(h,s)s^{\frac{\gamma+1}{m}}ds+R(h),$
(73)

where by the mean value theorem the remainder satisfies

$|R(h)|=h^{-\frac{(m+1)(\gamma+1)}{m}}|v^{\prime}(\widehat{\xi}_{h})|\int_{0}^{h}K(h,s)(h-s)s^{\frac{\gamma+1}{m}}ds\leq
K_{+}|v^{\prime}(\widehat{\xi}_{h})|B\left(\gamma+2,\frac{\gamma+1}{m}+1\right)h,$
(74)

where we have used the assumption of kernel boundedness (36). Therefore, after
substitution $s=z\sigma$ in (73) we have

$v(h)^{m+1}-v(h)h^{-\gamma}\int_{0}^{1}K(h,h\sigma)\sigma^{\frac{\gamma+1}{m}}d\sigma=R(h).$
(75)

Now, the left-hand side of the above is bounded from (36) while the right-hand
side vanishes to zero when $h\rightarrow 0^{+}$. We can then propose that the
starting step for the numerical method is

$v_{0}=v(0)=\lim\limits_{h\rightarrow
0^{+}}\left(h^{-\gamma}\int_{0}^{1}K(h,h\sigma)\sigma^{\frac{\gamma+1}{m}}d\sigma\right)^{\frac{1}{m}},$
(76)

provided the above limit exists. We thus assume that it is the case since
otherwise, the solution could be unnecessarily difficult or even meaningless
to solve numerically. Note that, taking the above as an initial step in the
scheme does not introduce any error apart from numerical rounding. In practice
we can use for example the Gaussian quadrature to compute the integral above
with finite $h$ fixed to be the numerical method step.

The initial step (76) is sufficient to start the rectangle scheme (62)
however, in order to initialize the second order trapezoidal method we need a
guess on the value of $v(h)$. A straightforward idea is to once again use the
linear interpolation between $v(0)$ and $v(h)$ similarly as in the odd case of
method’s construction. But in present case it yields $v(h)$ implicitly

$\begin{split}v(h)^{m+1}&=h^{-\frac{(m+1)(\gamma+1)}{m}}\int_{0}^{h}K(h,s)s^{\frac{\gamma+1}{m}}v(s)ds\\\
&=h^{-\frac{(m+1)(\gamma+1)}{m}}\left[v(0)\int_{0}^{h}K(h,s)s^{\frac{\gamma+1}{m}}\left(1-\frac{s}{h}\right)ds+v(h)\int_{0}^{h}K(h,s)s^{\frac{\gamma+1}{m}}\frac{s}{h}ds\right]+R(h),\end{split}$
(77)

where, similarly as before, the remainder is $R(h)=O(h^{2})$ as $h\rightarrow
0^{+}$. After truncation, the second initial condition $v_{1}$ can be found by
solving the nonlinear equation

$h^{\frac{(m+1)(\gamma+1)}{m}}v_{1}^{m+1}-\left(\int_{0}^{h}K(h,s)s^{\frac{\gamma+1}{m}}\frac{s}{h}ds\right)v_{1}-\left(\int_{0}^{h}K(h,s)s^{\frac{\gamma+1}{m}}\left(1-\frac{s}{h}\right)ds\right)v_{0}=0.$
(78)

Observe that solving the above is required only once at the initialization
phase of the second order method. Similarly, as with the computation of
$v_{0}$ we can use Gaussian quadrature to approximate the above integrals.
Moreover, since the the above nonlinear equation is of the form
$ax^{m+1}-bx-c=0$ with a known derivative, we can readily use Newton’s method
to solve it. Therefore, (76) and (78) supply us with at worst a second order
approximation of the starting values for the scheme (52) along with weights
(66) and (69).

###  3.2 Convergence

Before we proceed to the convergence proof of (52) with weights (66) and (69)
we present two auxiliary results that will be used later. The first one is a
variation on a discrete version of the classical Grönwall-Bellman’s lemma.
These lemmas are essential in studying convergence of numerical methods for
partial differential equations in both the classical and fractional setting.
In the literature one can find many different versions of them tailored for
specific needs. For example, they are an essential tool in the important case
of subdiffusion [73] when they are usually used to solve the error inequality.
In the cited work authors prove a general case of the discrete fractional
Grönwall-Bellman’s inequality where they allow for general weights present in
the sum. This generality helps to tackle convergence of a variety of numerical
schemes based on, for example, nonuniform L1 method or Alikhanov’s
discretisations. Other interesting results can be found in [38, 40].

######  Lemma 1.

Let $\left\\{e_{n}\right\\}$, $n=1,2,...$ be a sequence of positive numbers
satisfying

$e_{n}\leq\frac{\mu}{n}\sum_{i=1}^{n-1}e_{i}+\delta,\quad n\geq 1,$ (79)

where $\mu>0$ and $\delta>0$. Then, we have

$e_{n}\leq\delta f_{n},$ (80)

where

$f_{n}=\frac{\Gamma(n+\mu)}{n!}\sum_{k=0}^{n-1}\frac{k!}{\Gamma(k+1+\mu)}.$
(81)

Moreover,

$f_{n}\begin{cases}\leq\dfrac{1}{1-\mu},&0<\mu<1,\\\ \leq\ln n+1,&\mu=1,\\\
\stackrel{{\scriptstyle
n\rightarrow\infty}}{{\sim}}\dfrac{1}{\Gamma(\mu+1)}\dfrac{n^{\mu-1}}{\mu-1},&\mu>1.\\\
\end{cases}$ (82)

###### Proof.

We proceed by mathematical induction. To begin, let us observe that
$e_{1}\leq\delta$ by the convention that $\sum_{i=1}^{0}=0$. Hence, $f_{1}=1$.
Suppose that $e_{i}\leq\delta f_{i}$ for $1\leq i\leq n-1$. Due to this
inductive assumption we have

$e_{n}\leq\delta\left(\frac{\mu}{n}\sum_{i=1}^{n-1}f_{i}+1\right).$ (83)

We claim that the sequence $f_{i}$ defined in (81) satisfies the following
difference equation

$\frac{\mu}{n}\sum_{i=1}^{n-1}f_{i}+1=f_{n},$ (84)

which combined with (83) proves the inductive assertion. When we subtract (84)
written for $n+1$ the equation for $f_{n}$ we can obtain a local equation

$(n+1)f_{n+1}-(n+\mu)f_{n}-1=0.$ (85)

Switching back to $f_{n}$ and rearranging we arrive at the recurrence

$f_{n}=\left(1+\frac{\mu-1}{n}\right)f_{n-1}+\frac{1}{n}.$ (86)

The above can be iterated to yield

$f_{n}=\prod_{i=2}^{n}\left(1+\frac{\mu-1}{i}\right)+\sum_{i=0}^{n-2}\frac{1}{n-2}\prod_{j=1}^{i}\left(1+\frac{\mu-1}{n-j+1}\right).$
(87)

Further, the products can be simplified considerably when we notice that

$\prod_{i=2}^{n}\left(1+\frac{\mu-1}{i}\right)=\frac{(1+\mu)(2+\mu)...(n+\mu-1)}{n!}=\frac{1}{n!}\frac{\Gamma(n+\mu)}{\Gamma(\mu+1)},$
(88)

and similarly for the second one. Consequently, we arrive at

$\begin{split}f_{n}&=\frac{\Gamma(n+\mu)}{n!}\left(\frac{1}{\Gamma(\mu+1)}+\sum_{k=2}^{n}\frac{1}{k}\frac{k!}{\Gamma(k+\mu)}\right)=\frac{\Gamma(n+\mu)}{n!}\left(\frac{1}{\Gamma(\mu+1)}+\sum_{k=1}^{n-1}\frac{k!}{\Gamma(k+1+\mu)}\right)\\\
&=\frac{\Gamma(n+\mu)}{n!}\sum_{k=0}^{n-1}\frac{k!}{\Gamma(k+1+\mu)},\end{split}$
(89)

which is (81) and solves the recurrence.

In order to show (82) it suffices to consider the recurrence (84). Suppose
that $0<\mu<1$, then $f_{1}=1<1/(1-\mu)$. Further, if $f_{i}\leq 1/(\mu-1)$
then

$f_{n}\leq\frac{\mu}{n}\sum_{i=1}^{n-1}\frac{1}{1-\mu}+1\leq\frac{\mu}{1-\mu}+1=\frac{1}{1-\mu}.$
(90)

Next, for $\mu=1$ we trivially have $f_{1}=1=\ln 1+1$, and further by
inductive assumption

$f_{n}\leq\frac{1}{n}\sum_{i=1}^{n-1}\left(\ln
i+1\right)+1\leq\frac{1}{n}\int_{1}^{n}\left(\ln x+1\right)dx+1,$ (91)

where we have used the fact that the sum is bounded by an appropriate
integral. Evaluating, we obtain

$f_{n}\leq\frac{1}{n}\left.x\ln x\right|_{1}^{n}+1=\ln n+1.$ (92)

Next, for the case $\mu>1$ it is convenient to use the exact formula for
$f_{n}$, that is (81). Then

$f_{n}\sim
n^{\mu-1}\sum_{k=0}^{\infty}\frac{k!}{\Gamma(k+1+\mu)}\quad\text{as}\quad
n\rightarrow\infty,$ (93)

where we have used Stirling’s formula for the prefactor and to ascertain the
series convergence: $k!/\Gamma(k+1+\mu)\sim k^{\mu}$ as $k\rightarrow\infty$.
We can find the exact form of the above sum by using the relation between beta
and gamma functions

$\sum_{k=0}^{\infty}\frac{k!}{\Gamma(k+1+\mu)}=\frac{1}{\Gamma(\mu)}\sum_{k=0}^{\infty}\frac{\Gamma(k+1)\Gamma(\mu)}{\Gamma(k+1+\mu)}=\frac{1}{\Gamma(\mu)}\sum_{k=0}^{\infty}\int_{0}^{1}(1-x)^{\mu-1}x^{k}dx,$
(94)

and by Tonelli’s theorem we can exchange the order of summation and
integration

$\sum_{k=0}^{\infty}\frac{k!}{\Gamma(k+1+\mu)}=\frac{1}{\Gamma(\mu)}\int_{0}^{1}(1-x)^{\mu-1}\sum_{k=0}^{\infty}x^{k}dx=\frac{1}{\Gamma(\mu)}\int_{0}^{1}(1-x)^{\mu-2}dx=\frac{1}{(\mu-1)\Gamma(\mu)},$
(95)

which proves the last assertion. ∎

For a thorough presentation of similar results see [74]. Next, we prove
boundedness of the numerical approximation.

######  Lemma 2.

Let $v_{n}$ be the numerical approximation of $v(z_{n})$ calculated from (52).
Fix sufficiently small number $\epsilon>0$ and choose $h>0$ small enough for
$(z_{n})^{-(m+1)(\gamma+1)/m}\sum_{i=0}^{n-1}|\delta_{i}(h)|<\epsilon$. Then,
for $n\geq 1$ we have

$0<\left(K_{-}\left(\beta\left(\gamma+1,\frac{\gamma+1}{m}+1\right)-\epsilon\right)\right)^{\frac{1}{m}}\leq
v_{n}\leq\left(K_{-}\left(\beta\left(\gamma+1,\frac{\gamma+1}{m}+1\right)+\epsilon\right)\right)^{\frac{1}{m}},$
(96)

provided that the initial step $v_{0}$ satisfies the same inequality. Here,
$\beta(\cdot,\cdot)$ denotes the beta function.

###### Proof.

We will proceed by induction. The first step is satisfied by the assumption.
Suppose further that $v_{i}$ satisfies (96) for all $0\leq i<n-1$. We will
show that the same is true for $v_{n}$. To this end, write $V_{-}^{1/m}\leq
v_{i}\leq V_{+}^{1/m}$ where $V_{\pm}$ are defined in (96). By (52) and (36)
we have

$v_{n}^{m+1}\geq
V_{-}^{\frac{1}{m}}K_{-}z_{n}^{-\frac{(m+1)(\gamma+1)}{m}}h\sum_{i=0}^{n-1}w_{n,i}(z_{n}-z_{i})^{\gamma}z_{i}^{\frac{\gamma+1}{m}}.$
(97)

Now, according to (51) we have

$h\sum_{i=0}^{n-1}w_{n,i}(z_{n}-z_{i})^{\gamma}z_{i}^{\frac{\gamma+1}{m}}=\int_{0}^{z_{n}}(z_{n}-s)^{\gamma}s^{\frac{\gamma+1}{m}}ds=z_{n}^{\frac{(\gamma+1)(m+1)}{m}}\beta\left(\gamma+1,\frac{\gamma+1}{m}+1\right)-\sum_{i=0}^{n_{1}}\delta_{i}(h).$
(98)

Combining the two above expressions we obtain

$v_{n}\geq
V_{-}^{\frac{1}{m(m+1)}}\left(K_{-}\left(\beta\left(\gamma+1,\frac{\gamma+1}{m}+1\right)-z_{n}^{-\frac{(m+1)(\gamma+1)}{m}}\delta_{n}(h)\right)\right)^{\frac{1}{m+1}}.$
(99)

Now, since by the assumption the remainder is smaller than $\epsilon$ we can
write

$v_{n}\geq
V_{-}^{\frac{1}{m(m+1)}}\left(K_{-}\left(B\left(\gamma+1,\frac{\gamma+1}{m}+1\right)-\epsilon\right)\right)^{\frac{1}{m+1}}=V_{-}^{\frac{1}{m(m+1)}}V_{-}^{\frac{1}{m+1}}=V_{-}^{\frac{1}{m}},$
(100)

what completes the induction. The proof of the upper bound proceeds
identically. ∎

As it will be seen in the following main results, the boundedness from below
plays a crucial role in the proof.

A remark concerning the nature of the nonlinearity in (43) can be revealing.
Note that if we substitute $f=y^{m+1}$ then the equation transforms into

$f(z)=\int_{0}^{z}K(z,s)f(s)^{\frac{1}{m+1}}ds,$ (101)

in which the nonlinearity is manifestly nonlipschitzian. Therefore, we cannot
use the general theory to conclude the convergence. To wit, for the Lipschitz
nonlinearity we are always able to reduce the analysis of the method’s error
$e_{n}$ to investigations of the following recurrence inequality (see [57])

$\left|e_{n}\right|\leq\mu h\sum_{i=1}^{n-1}\left|e_{i}\right|+\delta(h),$
(102)

where $\mu$ and $\nu=\nu(h)$ are independent of $n$. Using the well-known
discrete version of the Grönwall-Bellman’s lemma we can solve the above to
yield (compare with our Lemma 1)

$\left|e_{n}\right|\leq\delta(h)e^{\mu nh}.$ (103)

The term $\delta(h)$ depends on the local consistency error and vanishes for
$h\rightarrow 0$. Moreover, since we are comparing the numerical and exact
solutions at a fixed point, we have $nh\rightarrow\text{const.}$. Therefore,
$e_{n}\rightarrow 0$ and the method is convergent. The lipschitzian character
of the nonlinearity makes the proof similar as in the linear equations. For
the non-Lipschitz case we cannot in general write (102) and thus the arguments
have to be different and more subtle. Notice that this difficulty is the very
consequence of the degeneracy of the main equation (1). We summarize the
result in the following theorem.

######  Theorem 2.

Fix $z\in(0,1]$ and chose the weights of the quadrature (51) according to (66)
and (69) while the starting values from (76) and (78). Then, if

$\mu_{m}:=\frac{4K_{+}}{(m+1)V_{-}}<3,$ (104)

the scheme (52) is convergent to the nontrivial solution of (47) with an order
at least equal to

$\min\\{2,3-\mu_{m}\\}.$ (105)

Here, $V_{-}^{1/m}$ is the lower bound for both $v=v(z)$ and $v_{n}$.

###### Proof.

Let us define the error of the numerical approximation by
$e_{n}=v(z_{n})-v_{n}$, where $v(z)$ is the solution of (47) while $v_{n}$
comes from (52). Now, the difference between these two equations is

$v(z_{n})^{m+1}-v_{n}^{m+1}=z_{n}^{-\frac{(m+1)(\gamma+1)}{m}}\left(\int_{0}^{z}K(z,s)s^{\frac{\gamma+1}{m}}v(s)ds-\sum_{i=0}^{n-1}w_{n,i}(h)v_{i}\right).$
(106)

With the use of (51) and the construction leading to (70) we further obtain

$v(z_{n})^{m+1}-v_{n}^{m+1}=z_{n}^{-\frac{(m+1)(\gamma+1)}{m}}\sum_{i=1}^{n-1}w_{n,i}(h)e_{i}+\delta_{n}(h),$
(107)

where the remainder $\delta_{n}(h)$ satisfies (71) and the zero term vanishes
due to exact starting value (76). On the other hand, by the mean value theorem
we can write

$v(z_{n})^{m+1}-v_{n}^{m+1}=(m+1)\xi_{n}^{m}e_{n},$ (108)

where $\xi_{n}$ lies between $v(z_{n})$ and $v_{n}$. Next, by combining the
two above equations we have

$(m+1)\xi_{n}^{m}|e_{n}|\leq
z_{n}^{-\frac{(m+1)(\gamma+1)}{m}}\sum_{i=1}^{n-1}|w_{n,i}(h)||e_{i}|+Ch^{2},$
(109)

where $C$ is a constant explicitly defined in (71). The boundedness of both
$v(z_{n})$ and $v_{n}$ (see (48) and Lemma 2)) implies that $\xi_{n}\geq
V_{-}$ for some constant and thus

$|e_{n}|\leq\frac{z_{n}^{-\frac{(m+1)(\gamma+1)}{m}}}{(m+1)V_{-}}\sum_{i=1}^{n-1}|w_{n,i}(h)||e_{i}|+\frac{Ch^{2}}{(m+1)V_{-}},$
(110)

which is a recurrence inequality which will eventually be solved with invoking
Grönwall-Bellman’s lemma. To see this observe that by (66) and (69) we have
for example

$|w_{n,i}(h)|\leq\frac{4K_{+}}{n}z_{n}^{\frac{(m+1)(\gamma+1)}{m}},$ (111)

where the factor $4/n$ comes from the fact that $|s-z_{i}|\leq 2h$ on each
subinterval of length $2/n$, $K_{+}$ is the bound on the kernel (36), and the
integrand is not larger than $1$. This leads to

$|e_{n}|\leq\frac{4K_{+}}{(m+1)V_{-}}\frac{1}{n}\sum_{i=1}^{n-1}|e_{i}|+\frac{Ch^{2}}{(m+1)V_{-}}=:\frac{\mu_{m}}{n}\sum_{i=1}^{n-1}|e_{i}|+\delta.$
(112)

Invoking Lemma 1 thus yields

$|e_{n}|\leq\frac{Ch^{2}}{(m+1)V_{-}}f_{n},$ (113)

where $f_{n}$ is defined in (81). Now, if $0<\mu_{m}<1$, the error $|e_{n}|$
is bounded for all $n$, i.e.

$|e_{n}|\leq\frac{Ch^{2}}{(m+1)V_{-}-4K_{+}}$ (114)

Further, for $\mu_{m}>1$ we have

$|e_{n}|\leq\frac{Ch^{2}}{(m+1)V_{-}}f_{n}\sim\frac{Ch^{2}}{(\mu_{m}-1)\Gamma(\mu_{m})(m+1)V_{-}}n^{\mu_{m}-1}\sim\frac{Cz^{\mu_{m}-1}}{(\mu_{m}-1)\Gamma(\mu_{m})(m+1)V_{-}}h^{3-\mu_{m}},$
(115)

as $n\rightarrow\infty$. Here, we have used the fact that $nh\rightarrow z$.
Similarly, the marginal case $\mu=1$ yields a convergence of order $-h^{2}\ln
h$. The proof is complete. ∎

According to (105) the method retains its order for sufficiently small
$\mu_{m}$. More specifically, when $\mu_{m}<1$ the method is second order
accurate. This estimate gradually becomes worse until $\mu_{m}=3$ when the
theorem does not guarantee convergence. Note, however, that choosing
sufficiently large $m$ we always can obtain a second order method. The fact
that, according to the theorem, the order of the method can be smaller than
the order of the quadrature is probably the consequence of severe nonlinearity
of the equation. This, in turn, corresponds to the degeneracy of the original
PDE (1) what throughout the years has proved to lay significant difficulties
for both theoretical and numerical studies (see for ex. [48, 52]). The main
difference between the classical case of Lipschitzian nonlinearity is the
occurrence of $1/n$ instead of $h$ in the nonlocal recurrence (112). This
forces us to use Lemma 1) rather than the classical discrete Grönwall-
Bellman’s lemma. Since $1/n$ changes in each recurrence step the assertion is
somewhat different yielding a loss in convergence order. However, for
$0<\mu_{m}<1$ the error is $O(h^{2})$ as $h\rightarrow 0^{+}$ for any $n>0$. A
very thorough account of numerical methods for integral equations with
Lipschitzian nonlinearities can be found in [57].

##  4 Numerical examples

In this section we present several numerical examples. We start with a simple
integral equation that can be solved exactly in a closed form. Next, we
consider a more complex equation and use it to verify the order of
convergence. Further, we proceed to solving the time-fractional diffusion
problems summarized in Tab. 1 and we end this section with temporal
computational complexity estimates and comparison with finite difference
scheme. Numerical calculations have been conducted in Julia programming
language. Integrals appearing in weights (66) and (69) have been calculated
with the QuadGK package that utilizes Gauss-Kronrod adaptive quadrature.

###  4.1 Exact constant solution

We start with a simple illustration of the exactness of our numerical methods.
Let $K(z,s)=(z-s)^{\gamma}$ in (47). According to (49) the integral equation
has then a constant solution (that is, a power function is a solution of
(43)). In Tab. 2 we have collected the maximal absolute error (first norm) of
the numerical approximation (52) with second order weights (66) and (69). Note
that in calculations we have used $N=10$ which is relatively small. Our
experiments showed that increasing the number of interval divisions does not
improve the error. Notice also that the smallest error $O(10^{-16})$ appears
for the case $\gamma=0$ and $m=1$ for which the exact solution is equal to
$1/2$. The majority of simulations concluded that the error is of order of
$O(10^{-11})$ we can conclude that the method performs very well and
reproduces the exact solution with good accuracy.

| $\gamma=0$ | $\gamma=0.5$ | $\gamma=\sqrt{2}$ | $\gamma=\pi$
---|---|---|---|---
$m=1$ | $1.11\times 10^{-16}$ | $2.90\times 10^{-11}$ | $8.94\times 10^{-13}$ | $2.20\times 10^{-14}$
$m=2$ | $2.40\times 10^{-10}$ | $1.18\times 10^{-10}$ | $1.19\times 10^{-11}$ | $1.60\times 10^{-13}$
$m=10$ | $1.90\times 10^{-10}$ | $1.01\times 10^{-10}$ | $4.26\times 10^{-11}$ | $4.71\times 10^{-11}$
$m=100$ | $1.85\times 10^{-11}$ | $2.03\times 10^{-11}$ | $1.97\times 10^{-11}$ | $1.01\times 10^{-11}$
Table 2: Maximal absolute error for a numerical approximation of the solution
of (47) with $K(z,s)=(z-s)^{\gamma}$ for various $\gamma$ and $m$. The number
of interval divisions is $N=10$.

###  4.2 Order of convergence

We will calculate the empirical order of convergence of our second order
method and compare it with the rectangle scheme (62). This illustration will
be completed with the use of the integral equation (47) with

$K(z,s)=\frac{\sqrt{z-s}}{1+\sin^{2}s}.$ (116)

Obviously, the kernel satisfies our assumptions on the boundedness (36) with
$K_{-}=1/2$ and $K_{+}=1$. The order of convergence is estimated with the
extrapolation technique known as Aitken’s method (see [57])

$\text{order}\approx\log_{2}\frac{\left|v^{(2N)}_{2N}-v_{N}^{(N)}\right|}{\left|v^{(4N)}_{4N}-v_{2N}^{(2N)}\right|},$
(117)

which compares the numerical solution evaluated at $z=1$ (the worst case) with
calculations for $N$, $2N$, and $4N$ steps. The summary of obtained results is
presented in Tab. 3. We can see that the estimated order is near 2 however,
for larger $m$ it is somewhat lower. This fact probably originates in the
necessity of taking a very high order root in the equation. Moreover, for this
example, we have $\mu_{m}=0.8<3$ what satisfies the assumption of Theorem 2.
We can thus see that, as anticipated, the method is of second order. Moreover,
similar calculations have been done for the rectangle method (62) and the
results were more uniform: the estimated order was equal to $0.99$ for all
considered $m$.

$m$ | 1 | 5 | 10 | 15 | 20 | 50 | 100
---|---|---|---|---|---|---|---
order | 2.02 | 2.00 | 1.96 | 1.92 | 1.90 | 1.86 | 1.84
Table 3: Numerically calculated order of convergence of the numerical method
based on trapezoidal approximation (66), (69). The base for Aitken’s
estimation is $N=100$.

###  4.3 Anomalous diffusion

Now, we proceed to the subject of the main interest - the time-fractional
porous medium equation. We use the kernel (29) which conveniently can be
written with the help of incomplete beta functions

$\begin{split}K(z,u)&=\frac{m+1}{\Gamma(1-\alpha)}\left[B(1-u)\beta(1-\alpha,a+2b+1,1-w)\right.\\\
&\left.+(A+B)\left((1-u)\beta(1-\alpha,a+2b+1,1-w)-(1-z)\beta(1-\alpha,a+b+1,1-w)\right)\right],\end{split}$
(118)

where

$\beta(a,b,z)=\int_{0}^{z}t^{a-1}(1-t)^{b-1}dt,$ (119)

and $A$, $B$, $a$, $b$ are chosen according to the boundary conditions
summarized in Tab. 1. Some exemplary plots of the self-similar solutions of
(1) are presented on Fig. 1.

Figure 1: Exemplary self-similar solutions of the time-fractional porous
medium equation (1). Here, $\alpha=0.5$ and $m=2$.

First, we estimate the minimal value of parameter $m$, say $m_{0}$, for which
the assumption of the Theorem 2 is satisfied for the most important -
Dirichlet - boundary condition. According to (104) it is equal to a zero of a
function $m\mapsto\mu_{m}-3$. In order to estimate that we have to find
$K_{+}$ which can be a-priori calculated from (33). However, a more accurate
bound comes from the estimate

$K_{+}\geq\sup_{0\leq u\leq z\leq 1}\frac{K(z,u)}{(z-u)^{1-\alpha}},$ (120)

which can be found numerically. Note also that when we compare the formulas
(33) and (96) we can expect that the ratio $K_{+}/V_{-}$ is constant with
respect to $m$. Whence $\mu_{m}$ should be $O(m^{-1})$ as
$m\rightarrow\infty$. Results of numerical estimation are presented in Tab. 4.
We can see that usually the value of $m_{0}$ is slightly larger than $1$ and
closes to $3$ for very small $\alpha$. In applications, one usually finds
$m\approx 7-8$ with $\alpha\in[0.6,1]$ [12] and, for these important cases, we
know that Theorem 2 guarantees convergence.

$\alpha$ | 0.99 | 0.9 | 0.8 | 0.7 | 0.6 | 0.5 | 0.4 | 0.3 | 0.2 | 0.1 | 0.01
---|---|---|---|---|---|---|---|---|---|---|---
$m_{0}$ | 1.10 | 1.14 | 1.19 | 1.25 | 1.28 | 1.35 | 1.66 | 1.98 | 2.31 | 2.62 | 2.84
Table 4: Estimated critical values of $m$ for which the assumptions of Theorem
2 are satisfied for the case of anomalous diffusion.

Similarly as above we can estimate the order of convergence of the trapezoid
method. We again use Aitken’s algorithm and present the results in Tab. 5. As
can be seen the estimated order stays near $2$ especially for moderate values
of $m$ what has also been observed in the previous example.

$m$ $\alpha$ | 0.1 | 0.3 | 0.5 | 0.7 | 0.9 | 0.99
---|---|---|---|---|---|---
1 | 1.99 | 2.00 | 2.00 | 1.98 | 2.08 | 2.07
3 | 1.96 | 1.99 | 1.99 | 1.98 | 1.97 | 1.96
5 | 1.93 | 1.96 | 1.96 | 1.95 | 2.01 | 1.92
7 | 1.90 | 1.93 | 1.93 | 1.93 | 1.95 | 1.88
10 | 1.87 | 1.91 | 1.91 | 1.91 | 1.92 | 1.85
15 | 1.86 | 1.89 | 1.90 | 1.90 | 1.91 | 1.81
Table 5: Order of convergence for the trapezoidal method applied to
subdiffusion for different $\alpha$ and $m$.

As a further verification of our method we can compute the wetting front
position $\eta^{*}$ for the classical case $\alpha=1$ and compare it with
results from [75]. In that work, the values of $\eta^{*}$ were computed with a
use of power series and we can treat them as exact (up to 9 decimal places).
Our calculations for $m=2$ are summarized in Tab. 6. Results for other values
of $m$ are very similar. Notice that even for a very small number of
subdivisions we obtain an error of order $O(10^{-5})$ which is a consequence
of the second order accuracy. Further, numerical simulations show that on a
log-log scale, the error $|\eta^{*}-\eta^{*}_{exact}|$ behaves as a line with
tangent $-2$. Therefore, we can expect that the
$|\eta^{*}-\eta^{*}_{exact}|=O(N^{-2})$ as $N\rightarrow\infty$.

$N$ | 10 | 20 | 50 | 100 | 200 | 500 | 1500
---|---|---|---|---|---|---|---
$|\eta^{*}-\eta^{*}_{exact}|$ | $1.1\times 10^{-4}$ | $2.9\times 10^{-5}$ | $4.6\times 10^{-6}$ | $1.1\times 10^{-6}$ | $2.8\times 10^{-7}$ | $4.5\times 10^{-8}$ | $5.4\times 10^{-9}$
Table 6: The error in calculating wetting front position for $\alpha=1$ and
$m=2$. The exact value has been taken from [76]. The error decays at a rate
$O(N^{-2})$ as $N\rightarrow\infty$.

Wetting fronts for subdiffusive case cannot easily be compared with exact
values. However, we have some useful asymptotics. First, let us consider the
Dirichlet boundary condition for which we know that (see [26])

$\eta^{*}=O(m^{-1/2})\quad m\rightarrow\infty.$ (121)

This relation can also be obtained by combining the value for $\eta^{*}$ taken
from Tab. 1 with estimates (37) and noting that $K_{+}\propto m+1$. Our
numerical simulations are depicted on Fig. 2. Notice that the above asymptotic
behaviour is evident even for small values of $m$, i.e. in a log-log scale all
lines become parallel to $m^{-1/2}$.

A similar reasoning can be applied to Neumann and Robin conditions. To this
end we need an asymptotic behaviour of the derivative $y^{\prime}(1)$ for
$m\rightarrow\infty$. First, notice that due to (33) and (37) we have
$\|y\|=O((1+m)^{1/m})=O(1)$ for large $m$. Then, from (26) we obtain

$\|y^{\prime}\|=O\left(\frac{1}{m+1}\right)\quad\text{as}\quad
m\rightarrow\infty,$ (122)

uniformly for $\alpha\in(0,1)$. Next, using the value from Tab. 1 we can write

$\eta^{*}=\left(\frac{1}{y^{m}(1)y^{\prime}(1)}\right)^{\frac{m}{m+2}}=O\left(\frac{m+1}{m+1}\right)^{\frac{m}{m+2}}=O(1)\quad
m\rightarrow\infty.$ (123)

for Neumann condition. Similarly, the wetting front for the Robin case also is
bounded for large $m$. We can see on Fig. 2 that this observation is confirmed
with numerical simulations, that is for large $m$ wetting fronts converge to
fixed values.

Figure 2: Wetting front position $\eta^{*}$ with respect to $m$ for a
Dirichlet (top), Neumann (middle), and Robin (bottom) boundary conditions
shown in a double logarithmic scale.

###  4.4 Comparison with finite difference method

In order to compare our method with some other popular schemes we conduct a
computation cost estimates in terms of the temporal complexity. We will
contrast the scheme based on the Volterra integral equation (52) with the
usual $\theta$-weighted finite difference scheme with L1 discretization of the
fractional derivative [10, 38] (however, the comparison with finite element or
spectral methods would yield similar results). The diffusivity is linearized
according to [77] and the method for solving (1) is following

$\begin{split}&u_{j}^{i+1}-\left(1-\theta\right)\frac{h^{\alpha}}{k^{2}\Gamma(2-\alpha)}\left(D^{i+1}_{j-1/2}u^{i+1}_{j-1}+\left(D^{i+1}_{j-1/2}+D^{i+1}_{j+1/2}\right)u^{i+1}_{j}+D^{i+1}_{j+1/2}u^{i+1}_{j+1}\right)\\\
&=-\sum_{k=1}^{i}a_{k,i}u^{k}_{j}+\theta\frac{h^{\alpha}}{k^{2}\Gamma(2-\alpha)}\left(D^{i}_{j-1/2}u^{i}_{j-1}+\left(D^{i}_{j-1/2}+D^{i}_{j+1/2}\right)u^{i}_{j}+D^{i}_{j+1/2}u^{i}_{j+1}\right),\end{split}$
(124)

where the weights are defined by

$a_{k,i}=(i-k+2)^{1-\alpha}-2(i-k+1)^{1-\alpha}+(i-k)^{1-\alpha}.$ (125)

and $D^{i}_{j\pm 1/2}$ is the linearised value of the diffusion coefficient

$\begin{split}D^{i}_{j\pm
1/2}=\frac{1}{2}&\left((u^{m})^{i}_{j}+m(u^{m-1})^{i}_{j}\left(u^{i}_{j}-u^{i-1}_{j}\right)+\right.\\\
&\left.(u^{m})^{i}_{j\pm 1/2}+m(u^{m-1})^{i}_{j\pm 1/2}\left(u^{i}_{j\pm
1/2}-u^{i-1}_{j\pm 1/2}\right)\right).\end{split}$ (126)

Here, $u^{i}_{j}$ is the numerical approximation of $u(x_{j},t_{i})$, where
$x_{j}=j\Delta x$ and $t_{i}=i\Delta t$.

Suppose we would like to compute the wetting front position in the Dirichlet
problem with a tolerance $\epsilon>0$ (but any other value of $u$ would serve
the same purpose for this benchmark). According to Tab. 1 this requires
finding $y(1)=v(1)$. Note that other boundary conditions would need
$y^{\prime}(1)$ which does not change our reasoning. Since the method is of
second order, we would like to choose $N>0$ in order to satisfy

$|V_{N}-v(1)|\leq\frac{C}{N^{2}}\leq\epsilon,$ (127)

which immediately gives $N=O(\epsilon^{-1/2})$ as $\epsilon\rightarrow 0^{+}$.
Due to nonlocality of the problem (52), each value $y_{n}$ arises from all
previous ones with $N(N+1)/2$ additions and multiplications with weights
$w_{n,i}(h)$. These, in turn have to be calculated by integration for which we
assume a constant cost $c_{I}$. Finally, in each step of the iteration we take
a $m+1$-th root. The method is started with initial values (76) and (78) which
require $3$ integrations, one root taking, and one solution of the nonlinear
algebraic equation which costs $c_{E}$. Therefore, the total cost of the
scheme can be estimated as

$\text{temporal complexity of
(\ref{eqn:NumMet})}=\left(2+c_{I}\right)\frac{N(N+1)}{2}+N+1+3c_{I}+c_{E}=O(N^{2})=O(\epsilon^{-1}),$
(128)

as $\epsilon\rightarrow 0^{+}$. Hence, we obtain a linear relationship between
calculating the wetting front and the prescribed tolerance.

As for the finite difference we assume that we can choose $\theta$ in order
for the scheme to be unconditionally stable (this is not trivial since the
equation we are solving is nonlinear and degenerate). We are thus allowed to
choose $\Delta x=O(\Delta t)$ with $\Delta t=T/N$, where the final time is
denoted by $T$. The wetting front $x^{*}(T)$ is approximated by iterating the
scheme up until the final time and then finding $j^{*}$ such that
$u_{j^{*}}^{N}>0$ and $u_{j+1}^{N}=0$. Since the fractional derivative in our
finite differences is discretized using the L1 scheme, the method will have
temporal order at most equal to $2-\alpha$ (the order depends on the
regularity of the solution, see [78]), so that in the best case we have

$|u^{N}_{j^{*}}-u(x^{*}(T),T)|=|u^{N}_{j^{*}}|\leq C\left(\Delta
t^{2-\alpha}+\Delta x^{2}\right)\leq C\Delta
t^{2-\alpha}=C\left(\frac{T}{N}\right)^{2-\alpha}<\epsilon,$ (129)

so that we should have at least $N=O(\epsilon^{-1/(2-\alpha)})$ as
$\epsilon\rightarrow 0^{+}$. Now, similarly as above we can estimate the
temporal complexity of the algorithm. Let $c_{F}$ denote the fixed cost of
each time step of the finite difference scheme (124), i.e. it contains all
additions and multiplications needed to advance in time. Because of the
fractional derivative, in order to arrive at step $N$ we have to evaluate dot
products of previous solutions $u^{k}_{j}$ with the weights $a_{k,i}$. This
gives rise to $2\times N(N+1)/2$ floating point operations. Moreover, since
the method is implicit in the stable case, in each time step we have to solve
a tridiagonal system which costs $O(N)$ operations. Therefore, we can estimate
the total cost of the method

$\text{temporal complexity of
(\ref{eqn:FiniteDifference})}=N(N+1)c_{F}+N\times
O(N)=O(N^{2})=O(\epsilon^{-\frac{2}{2-\alpha}}),$ (130)

as $\epsilon\rightarrow 0^{+}$. As we can see, the computational cost of
calculating the wetting front is always higher for the finite difference case
than it is for our method. Note also that we have assumed the best case for
the former method, that is sufficient regularity of the solution guaranteeing
$2-\alpha$ temporal order, and unconditional stability that allowed us to
choose $\Delta x$ and $\Delta t$ of the same order. Thorough numerical
calculations supporting this claim has been conducted in [26] which indicate
that the method based on Volterra equation is superior.

##  5 Conclusion

We have constructed a convergent second order method for solving (43) which
can encompass self-similar solutions of a time-fractional porous medium
equation on the half-line. This fact is a consequence of a series of
transformations that changed a nonlocal nonlinear PDE into an ordinary
Volterra integral equation. The interesting feature of the latter is a non-
Lipschitz nonlinearity that possess some difficulties in numerical analysis.

Our second order method is based on a linear (trapezoidal) reconstruction has
been applied to several examples. Numerical calculations confirmed convergence
with desired accuracy. We have observed that it suffices to use a relatively
small number of interval subdivisions in order to obtain a decent
approximation of the exact solution. Moreover, the method reproduced the
asymptotic behaviour of the wetting front for large values of $m$ in three
considered boundary conditions. All calculations have been conducted on a
personal computer with a four core processor. Each simulation took at most few
tens of seconds. This can be compared with our previous observation made in
[26] stating that our numerical method is much faster that the implicit finite
difference scheme applied to (1). The Reader is refereed to that work for
concrete computation times for these two approaches. The nonlocality,
nonlinearity, stiffness, and degeneracy of the governing equation require much
computing power to resolve. Transforming the original PDE into an ordinary
Volterra equation reduces one degree of freedom and simplifies the free-
boundary problem which facilitates the method’s construction and performance.
Therefore, our trapezoidal method is an accurate and fast way of computing
self-similar solutions of the time-fractional porous medium equation on the
half-line.

## 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] R. Metzler, J. Klafter, The random walk’s guide to anomalous diffusion: a fractional dynamics approach, Physics reports 339 (1) (2000) 1–77.
  * [2] R. Klages, G. Radons, I. M. Sokolov, Anomalous transport: foundations and applications, John Wiley & Sons, 2008.
  * [3] J. Klafter, S. Lim, R. Metzler, Fractional dynamics: recent advances, World Scientific, 2012.
  * [4] S. Schaufler, W. Schleich, V. Yakovlev, Scaling and asymptotic laws in subrecoil laser cooling, EPL (Europhysics Letters) 39 (4) (1997) 383.
  * [5] D. del Castillo-Negrete, B. Carreras, V. Lynch, Nondiffusive transport in plasma turbulence: a fractional diffusion approach, Physical Review Letters 94 (6) (2005) 065003.
  * [6] M. Levandowsky, B. White, F. Schuster, Random movements of soil amebas, Acta Protozoologica 36 (1997) 237–248.
  * [7] T. Sungkaworn, M.-L. Jobin, K. Burnecki, A. Weron, M. J. Lohse, D. Calebiro, Single-molecule imaging reveals receptor–g protein interactions at cell surface hot spots, Nature 550 (7677) (2017) 543.
  * [8] A. de Pablo, F. Quiros, A. Rodriguez, J. L. Vazquez, A fractional porous medium equation, Advances in Mathematics 226 (2) (2011) 1378–1409.
  * [9] Y. Pachepsky, D. Timlin, W. Rawls, Generalized richards’ equation to simulate water transport in unsaturated soils, Journal of Hydrology 272 (1) (2003) 3–13.
  * [10] Ł. 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) (2014) 1219–1237.
  * [11] A. E.-G. El Abd, J. J. Milczarek, Neutron radiography study of water absorption in porous building materials: anomalous diffusion analysis, Journal of Physics D: Applied Physics 37 (16) (2004) 2305.
  * [12] H. Sun, M. M. Meerschaert, Y. Zhang, J. Zhu, W. Chen, A fractal Richards equation to capture the non-boltzmann scaling of water transport in unsaturated media, Advances in Water Resources 52 (2013) 292–295.
  * [13] A. Head, Dislocation group dynamics iii. similarity solutions of the continuum approximation, Philosophical Magazine 26 (1) (1972) 65–72.
  * [14] I. Podlubny, Fractional differential equations: an introduction to fractional derivatives, fractional differential equations, to methods of their solution and some of their applications, Vol. 198, Academic press, 1998.
  * [15] L. Caffarelli, J. L. Vazquez, Nonlinear porous medium flow with fractional potential pressure, Archive for Rational Mechanics and Analysis 2 (202) (2011) 537–565.
  * [16] J. L. Vázquez, The mathematical theories of diffusion: Nonlinear and fractional diffusion, in: Nonlocal and Nonlinear Diffusions and Interactions: New Methods and Directions, Springer, 2017, pp. 205–278.
  * [17] Ł. Płociniczak, Analytical studies of a time-fractional porous medium equation. derivation, approximation and applications, Communications in Nonlinear Science and Numerical Simulation 24 (1) (2015) 169–183.
  * [18] Ł. Płociniczak, Derivation of the nonlocal pressure form of the fractional porous medium equation in the hydrological setting, Communications in Nonlinear Science and Numerical Simulation 76 (2019) 66–70.
  * [19] A. El Abd, S. Kichanov, M. Taman, K. Nazarov, D. Kozlenko, W. M. Badawy, Determination of moisture distributions in porous building bricks by neutron radiography, Applied Radiation and Isotopes 156 (2020) 108970.
  * [20] P. Biler, C. Imbert, G. Karch, The nonlocal porous medium equation: Barenblatt profiles and other weak solutions, Archive for Rational Mechanics and Analysis 215 (2) (2015) 497–529.
  * [21] J.-D. Djida, J. J. Nieto, I. Area, Nonlocal time porous medium equation with fractional time derivative, Revista Matemática Complutense 32 (2) (2019) 273–304.
  * [22] Q. Fan, G.-C. Wu, H. Fu, A note on function space and boundedness of the general fractional integral in continuous time random walk, Journal of Nonlinear Mathematical Physics (2021) 1–8.
  * [23] H. Fu, G.-C. Wu, G. Yang, L.-L. Huang, Continuous time random walk to a general fractional fokker–planck equation on fractal media, The European Physical Journal Special Topics (2021) 1–7.
  * [24] Ł. Płociniczak, M. Świtała, Existence and uniqueness results for a time-fractional nonlinear diffusion equation, Journal of Mathematical Analysis and Applications 462(2) (2018) 1425–1434.
  * [25] Ł. Płociniczak, Diffusivity identification in a nonlinear time-fractional diffusion equation, Fractional Calculus and Applied Analysis 19 (2016) 883–866.
  * [26] Ł. Płociniczak, Numerical method for the time-fractional porous medium equation, SIAM Journal on Numerical Analysis 57 (2) (2019) 638–656.
  * [27] K. Diethelm, N. J. Ford, Analysis of fractional differential equations, Journal of Mathematical Analysis and Applications 265 (2) (2002) 229–248.
  * [28] D. Baleanu, K. Diethelm, E. Scalas, J. J. Trujillo, Fractional calculus: models and numerical methods, Vol. 3, World Scientific, 2012.
  * [29] C. Tadjeran, M. M. Meerschaert, H.-P. Scheffler, A second-order accurate numerical approximation for the fractional diffusion equation, Journal of Computational Physics 213 (1) (2006) 205–213.
  * [30] X. Li, C. Xu, A space-time spectral method for the time fractional diffusion equation, SIAM Journal on Numerical Analysis 47 (3) (2009) 2108–2131.
  * [31] S. B. Yuste, L. Acedo, An explicit finite difference method and a new von neumann-type stability analysis for fractional diffusion equations, SIAM Journal on Numerical Analysis 42 (5) (2005) 1862–1874.
  * [32] F. Liu, P. Zhuang, V. Anh, I. Turner, K. Burrage, Stability and convergence of the difference methods for the space–time fractional advection–diffusion equation, Applied Mathematics and Computation 191 (1) (2007) 12–20.
  * [33] T. Langlands, B. I. Henry, The accuracy and stability of an implicit solution method for the fractional diffusion equation, Journal of Computational Physics 205 (2) (2005) 719–736.
  * [34] E. Cuesta, C. Lubich, C. Palencia, Convolution quadrature time discretization of fractional diffusion-wave equations, Mathematics of Computation 75 (254) (2006) 673–696.
  * [35] Y. Lin, C. Xu, Finite difference/spectral approximations for the time-fractional diffusion equation, Journal of computational physics 225 (2) (2007) 1533–1552.
  * [36] P. Zhuang, F. Liu, I. Turner, V. Anh, Galerkin finite element method and error analysis for the fractional cable equation, Numerical Algorithms 72 (2) (2016) 447–466.
  * [37] D. Li, C. Wu, Z. Zhang, Linearized Galerkin FEMs for nonlinear time fractional parabolic problems with non-smooth solutions in time direction, Journal of Scientific Computing 80 (1) (2019) 403–419.
  * [38] D. Li, H.-L. Liao, W. Sun, J. Wang, J. Zhang, Analysis of L1-Galerkin FEMs for time-fractional nonlinear parabolic problems, Communications in Computational Physics 24 (1) (2018) 86–103.
  * [39] H.-l. Liao, D. Li, J. Zhang, Sharp error estimate of the nonuniform L1 formula for linear reaction-subdiffusion equations, SIAM Journal on Numerical Analysis 56 (2) (2018) 1112–1133.
  * [40] S. McKee, Generalised discrete gronwall lemmas, ZAMM-Journal of Applied Mathematics and Mechanics/Zeitschrift für Angewandte Mathematik und Mechanik 62 (9) (1982) 429–434.
  * [41] J. R. Webb, Weakly singular gronwall inequalities and applications to fractional differential equations, Journal of Mathematical Analysis and Applications 471 (1-2) (2019) 692–711.
  * [42] M. A. Zaky, A. S. Hendy, J. E. Macías-Díaz, Semi-implicit Galerkin–Legendre spectral schemes for nonlinear time-space fractional diffusion–reaction equations with smooth and nonsmooth solutions, Journal of Scientific Computing 82 (1) (2020) 1–27.
  * [43] F. del Teso, Finite difference method for a fractional porous medium equation, Calcolo 51 (4) (2014) 615–638.
  * [44] L. Caffarelli, L. Silvestre, An extension problem related to the fractional laplacian, Communications in Partial Differential Equations 32 (8) (2007) 1245–1260.
  * [45] F. del Teso, J. Endal, E. R. Jakobsen, Robust numerical methods for local and nonlocal equations of porous medium type. part i: Theory, arXiv preprint arXiv:1801.07148 (2018).
  * [46] F. J. Gaspar, C. Rodrigo, Multigrid waveform relaxation for the time-fractional heat equation, SIAM Journal on Scientific Computing 39 (4) (2017) A1201–A1224.
  * [47] A. Bhrawy, A Jacobi spectral collocation method for solving multi-dimensional nonlinear fractional sub-diffusion equations, Numerical Algorithms 73 (1) (2016) 91–113.
  * [48] E. DiBenedetto, D. Hoff, An interface tracking algorithm for the porous medium equation, Transactions of the American Mathematical Society 284 (2) (1984) 463–500.
  * [49] D. Stan, F. del Teso, J. L. Vázquez, Finite and infinite speed of propagation for porous medium equations with fractional pressure, Comptes Rendus Mathematique 352 (2) (2014) 123–128.
  * [50] T. Arbogast, M. F. Wheeler, A nonlinear mixed finite element method for a degenerate parabolic equation arising in flow in porous media, SIAM Journal on Numerical Analysis 33 (4) (1996) 1669–1687.
  * [51] C. Ebmeyer, Error estimates for a class of degenerate parabolic equations, SIAM Journal on Numerical Analysis 35 (3) (1998) 1095–1112.
  * [52] E. Etienne, D. Šiška, Full discretization of the porous medium/fast diffusion equation based on its very weak formulation, Communications in Mathematical Sciences 10 (4) (2012) 1055–1080.
  * [53] I. S. Pop, W.-A. Yong, A numerical approach to degenerate parabolic equations, Numerische Mathematik 92 (2) (2002) 357–381.
  * [54] M. Eisenmann, E. Hansen, Convergence analysis of domain decomposition based time integrators for degenerate parabolic equations, Numerische mathematik 140 (4) (2018) 913–938.
  * [55] H. Brunner, Volterra Integral Equations: An Introduction to Theory and Applications, Vol. 30, Cambridge University Press, 2017.
  * [56] C. T. Baker, A perspective on the numerical treatment of Volterra equations, Journal of computational and applied mathematics 125 (1-2) (2000) 217–249.
  * [57] P. Linz, Analytical and numerical methods for Volterra equations, Vol. 7, SIAM, 1985.
  * [58] K. Frischmuth, N. J. Ford, J. T. Edwards, Volterra integral equations with non-lipschitz nonlinearity, in: Rostocker Mathematisches Kolloquium, Vol. 51, 1997, pp. 65–82.
  * [59] E. Buckwar, Iterative Approximation of the Positive Solutions of a Class of Nonlinear Volterra-type Integral Equations, Logos Verlag, 1997.
  * [60] Ł. Płociniczak, H. Okrasińska-Płociniczak, Numerical method for Volterra equation with a power-type nonlinearity, Applied Mathematics and Computation 337 (2018) 452–460.
  * [61] E. Buckwar, On a nonlinear Volterra integral equation, in: Volterra equations and applications, CRC Press, 2000, pp. 157–162.
  * [62] M. Arias, R. Benítez, V. Bolós, Non-lipschitz homogeneous Volterra integral equations, in: Modern Mathematics and Mechanics, Springer, 2019, pp. 237–259.
  * [63] R. Gazizov, A. Kasatkin, S. Y. Lukashchuk, Symmetry properties of fractional diffusion equations, Physica Scripta 2009 (T136) (2009) 014016.
  * [64] R. K. Gazizov, A. A. Kasatkin, Construction of exact solutions for fractional order differential equations by the invariant subspace method, Computers & Mathematics with Applications 66 (5) (2013) 576–584.
  * [65] R. Gorenflo, Y. Luchko, F. Mainardi, Wright functions as scale-invariant solutions of the diffusion-wave equation, Journal of Computational and Applied Mathematics 118 (1-2) (2000) 175–191.
  * [66] V. S. Kiryakova, Generalized fractional calculus and applications, CRC Press, 1993\.
  * [67] J.-D. Djida, J. J. Nieto, I. Area, Nonlocal time-porous medium equation: weak solutions and finite speed of propagation, Discrete Continuous Dyn. Syst. Ser. B (2018).
  * [68] F. Atkinson, L. Peletier, Similarity profiles of flows through porous media, Archive for Rational Mechanics and Analysis 42 (5) (1971) 369–379.
  * [69] E. Buckwar, Existence and uniqueness of solutions of Abel integral equations with power-law non-linearities, Nonlinear Analysis: Theory, Methods & Applications 63 (1) (2005) 88–96.
  * [70] G. Gripenberg, Unique solutions of some Volterra integral equations, Mathematica Scandinavica 48 (1) (1981) 59–67.
  * [71] P. Bushell, W. Okrasinski, Nonlinear Volterra integral equations with convolution kernel, Journal of the London Mathematical Society 2 (3) (1990) 503–510.
  * [72] W. Okrasinski, On nontrivial solutions to some nonlinear ordinary differential equations, Journal of Mathematical Analysis and Applications 190 (2) (1995) 578–583.
  * [73] H.-l. Liao, W. McLean, J. Zhang, A discrete gronwall inequality with applications to numerical schemes for subdiffusion problems, SIAM Journal on Numerical Analysis 57 (1) (2019) 218–237.
  * [74] W. F. Ames, B. Pachpatte, Inequalities for differential and integral equations, Vol. 197, Academic press, 1997.
  * [75] W. Okrasiński, On approximate solutions to some nonlinear diffusion problems, Zeitschrift für angewandte Mathematik und Physik ZAMP 44 (4) (1993) 722–731.
  * [76] W. Okrasiński, S. Vila, Power series solutions to some nonlinear diffusion problems, Zeitschrift für angewandte Mathematik und Physik ZAMP 44 (6) (1993) 988–997.
  * [77] J. Crank, Free and moving boundary problems, Oxford University Press, 1987.
  * [78] N. 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) (2019) 2135–2155.