Raytcho Lazarov

Bangti Jin, Raytcho Lazarov, Zhi Zhou

The subdiffusion equation with a Caputo fractional derivative of order $α\in(0,1)$ in time arises in a wide variety of practical applications, and it is often adopted to model anomalous subdiffusion processes in heterogeneous media. The L1 scheme is one of the most popular and successful numerical methods for discretizing the Caputo fractional derivative in time. The scheme was analyzed earlier independently by Lin and Xu (2007) and Sun and Wu (2006), and an $O(τ^{2-α})$ convergence rate was established, under the assumption that the solution is twice continuously differentiable in time. However, in view of the smoothing property of the subdiffusion equation, this regularity condition is restrictive, since it does not hold even for the homogeneous problem with a smooth initial data. In this work, we revisit the error analysis of the scheme, and establish an $O(τ)$ convergence rate for both smooth and nonsmooth initial data. The analysis is valid for more general sectorial operators. In particular, the L1 scheme is applied to one-dimensional space-time fractional diffusion equations, which involves also a Riemann-Liouville derivative of order $β\in(3/2,2)$ in space, and error estimates are provided for the fully discrete scheme. Numerical experiments are provided to verify the sharpness of the error estimates, and robustness of the scheme with respect to data regularity.

Bangti Jin, Raytcho Lazarov, Zhi Zhou

We consider the initial boundary value problem for the homogeneous time-fractional diffusion equation $\partial^α_t u - \De u =0$ ($0< α< 1$) with initial condition $u(x,0)=v(x)$ and a homogeneous Dirichlet boundary condition in a bounded polygonal domain $Ω$. We shall study two semidiscrete approximation schemes, i.e., Galerkin FEM and lumped mass Galerkin FEM, by using piecewise linear functions. We establish optimal with respect to the regularity of the solution error estimates, including the case of nonsmooth initial data, i.e., $v \in L_2(Ω)$.

Bangti Jin, Raytcho Lazarov, Zhi Zhou

Over the past few decades, there has been substantial interest in evolution equations that involving a fractional-order derivative of order $α\in(0,1)$ in time, due to their many successful applications in engineering, physics, biology and finance. Thus, it is of paramount importance to develop and to analyze efficient and accurate numerical methods for reliably simulating such models, and the literature on the topic is vast and fast growing. The present paper gives a concise overview on numerical schemes for the subdiffusion model with nonsmooth problem data, which are important for the numerical analysis of many problems arising in optimal control, inverse problems and stochastic analysis. We focus on the following aspects of the subdiffusion model: regularity theory, Galerkin finite element discretization in space, time-stepping schemes (including convolution quadrature and L1 type schemes), and space-time variational formulations, and compare the results with that for standard parabolic problems. Further, these aspects are showcased with illustrative numerical experiments and complemented with perspectives and pointers to relevant literature.

Beiping Duan, Raytcho Lazarov, Joeseph Pasciak

In this paper, we develop and study algorithms for approximately solving the linear algebraic systems: $\mathcal{A}_h^αu_h = f_h$, $ 0< α<1$, for $u_h, f_h \in V_h$ with $V_h$ a finite element approximation space. Such problems arise in finite element or finite difference approximations of the problem $ \mathcal{A}^αu=f$ with $\mathcal{A}$, for example, coming from a second order elliptic operator with homogeneous boundary conditions. The algorithms are motivated by the recent method of Vabishchevich, 2015, that relates the algebraic problem to a solution of a time-dependent initial value problem on the interval $[0,1]$. Here we develop and study two time stepping schemes based on diagonal Padé approximation to $(1+x)^{-α}$. The first one uses geometrically graded meshes in order to compensate for the singular behavior of the solution for $t$ close to $0$. The second algorithm uses uniform time stepping but requires smoothness of the data $f_h$ in discrete norms. For both methods, we estimate the error in terms of the number of time steps, with the regularity of $f_h$ playing a major role for the second method. Finally, we present numerical experiments for $\mathcal{A}_h$ coming from the finite element approximations of second order elliptic boundary value problems in one and two spatial dimensions.

Emilia Bazhlekova, Bangti Jin, Raytcho Lazarov et al.

We study the Rayleigh-Stokes problem for a generalized second-grade fluid which involves a Riemann-Liouville fractional derivative in time, and present an analysis of the problem in the continuous, space semidiscrete and fully discrete formulations. We establish the Sobolev regularity of the homogeneous problem for both smooth and nonsmooth initial data $v$, including $v\in L^2(Ω)$. A space semidiscrete Galerkin scheme using continuous piecewise linear finite elements is developed, and optimal with respect to initial data regularity error estimates for the finite element approximations are derived. Further, two fully discrete schemes based on the backward Euler method and second-order backward difference method and the related convolution quadrature are developed, and optimal error estimates are derived for the fully discrete approximations for both smooth and nonsmooth initial data. Numerical results for one- and two-dimensional examples with smooth and nonsmooth initial data are presented to illustrate the efficiency of the method, and to verify the convergence theory.

Stanislav Harizanov, Raytcho Lazarov, Pencho Marinov et al.

In this paper we consider efficient algorithms for solving the algebraic equation ${\mathcal A}^α{\bf u}={\bf f}$, $0< α<1$, where ${\mathcal A}$ is a symmetric and positive definite matrix obtained form finite difference or finite element approximations of second order elliptic problems in ${\mathbb R}^d$, $d=1,2,3$. The method is based on the best uniform rational approximation of the function $t^{β-α}$ for $0 < t \le 1$ and natural $β$, and the assumption that one has at hand an efficient method (e.g. multigrid, multilevel, or other fast algorithm) for solving equations like $({\mathcal A} +c {\mathcal I}){\bf u}= {\bf f}$, $c \ge 0$. The provided numerical experiments on model problems with ${\mathcal A}$ obtained by finite element approximation of elliptic equations in one and three spacial dimensions confirm the efficiency of the proposed algorithms.

Bangti Jin, Raytcho Lazarov, Zhi Zhou

In this work, we develop variational formulations of Petrov-Galerkin type for one-dimensional fractional boundary value problems involving either a Riemann-Liouville or Caputo derivative of order $α\in(3/2, 2)$ in the leading term and both convection and potential terms. They arise in the mathematical modeling of asymmetric super-diffusion processes in heterogeneous media. The well-posedness of the formulations and sharp regularity pickup of the variational solutions are established. A novel finite element method is developed, which employs continuous piecewise linear finite elements and "shifted" fractional powers for the trial and test space, respectively. The new approach has a number of distinct features: It allows deriving optimal error estimates in both $L^2(D)$ and $H^1(D)$ norms; and on a uniform mesh, the stiffness matrix of the leading term is diagonal and the resulting linear system is well conditioned. Further, in the Riemann-Liouville case, an enriched FEM is proposed to improve the convergence. Extensive numerical results are presented to verify the theoretical analysis and robustness of the numerical scheme.

Bangti Jin, Raytcho Lazarov, Xiliang Lu et al.

We consider a boundary value problem involving a Riemann-Liouville fractional derivative of order $α\in (3/2,2)$ on the unit interval $(0,1)$. The standard Galerkin finite element approximation converges slowly due to the presence of singularity term $x^{α-1}$ in the solution representation. In this work, we develop a simple technique, by transforming it into a second-order two-point boundary value problem with nonlocal low order terms, whose solution can reconstruct directly the solution to the original problem. The stability of the variational formulation, and the optimal regularity pickup of the solution are analyzed. A novel Galerkin finite element method with piecewise linear or quadratic finite elements is developed, and $L^2(D)$ error estimates are provided. The approach is then applied to the corresponding fractional Sturm-Liouville problem, and error estimates of the eigenvalue approximations are given. Extensive numerical results fully confirm our theoretical study.

Guido Kanschat, Raytcho Lazarov, Youli Mao

We apply geometric multigrid methods for the finite element approximation of flow problems governed by Darcy and Brinkman systems used in modeling highly heterogeneous porous media. The method is based on divergence-conforming discontinuous Galerkin methods and overlapping, patch based domain decomposition smoothers. We show in benchmark experiments that the method is robust with respect to mesh size and contrast of permeability for highly heterogeneous media.

Bangti Jin, Raytcho Lazarov, Joseph Pasciak et al.

We investigate semi-discrete numerical schemes based on the standard Galerkin and lumped mass Galerkin finite element methods for an initial-boundary value problem for homogeneous fractional diffusion problems with non-smooth initial data. We assume that $Ω\subset \mathbb{R}^d$, $d=1,2,3$ is a convex polygonal (polyhedral) domain. We theoretically justify optimal order error estimates in $L_2$- and $H^1$-norms for initial data in $H^{-s}(Ω),~0\le s \le 1$. We confirm our theoretical findings with a number of numerical tests that include initial data $v$ being a Dirac $δ$-function supported on a $(d-1)$-dimensional manifold.

Raytcho Lazarov, Petr Vabishchevich

We consider the homogeneous equation ${\mathcal A} u=0$, where ${\mathcal A}$ is a symmetric and coercive elliptic operator in $H^1(Ω)$ with $Ω$ bounded domain in ${\mathbb R}^d$. The boundary conditions involve fractional power $α$, $ 0 < α<1$, of the Steklov spectral operator arising in Dirichlet to Neumann map. For such problems we discuss two different numerical methods: (1) a computational algorithm based on an approximation of the integral representation of the fractional power of the operator and (2) numerical technique involving an auxiliary Cauchy problem for an ultra-parabolic equation and its subsequent approximation by a time stepping technique. For both methods we present numerical experiment for a model two-dimensional problem that demonstrate the accuracy, efficiency, and stability of the algorithms.

Stanislav Harizanov, Raytcho Lazarov, Pencho Marinov et al.

We discuss, study, and compare experimentally three methods for solving the system of algebraic equations $\mathbb{A}^α\bf{u}=\bf{f}$, $0< α<1$, where $\mathbb{A}$ is a symmetric and positive definite matrix obtained from finite difference or finite element approximations of second order elliptic problems in $\mathbb{R}^d$, $d=1,2,3$. The first method, introduced by Harizanov et.al, based on the best uniform rational approximation (BURA) $r_α(t)$ of $t^{1-α}$ for $0 \le t \le 1$, is used to get the rational approximation of $t^{-α}$ in the form $t^{-1}r_α(t)$. Here we develop another method, denoted by R-BURA, that is based on the best rational approximation $r_{1-α}(t)$ of $t^α$ on the interval $[0,1]$ and approximates $t^{-α}$ via $r^{-1}_{1-α}(t)$. The third method, introduced and studied by Bonito and Pasciak, is based on an exponentially convergent quadrature scheme for the Dundord-Taylor integral representation of the fractional powers of elliptic operators. All three methods reduce the solution of the system $\mathbb{A}^α\bf{u}=\bf{f}$ to solving a number of equations of the type $(\mathbb{A} +c\mathbb{I})\bf{u}= \bf{f}$, $c \ge 0$. Comprehensive numerical experiments on model problems with $\mathbb A$ obtained by approximation of elliptic equations in one and two spatial dimensions are used to compare the efficiency of these three algorithms depending on the fractional power $α$. The presented results prove the concept of the new R-BURA method, which performs well for $α$ close to $1$ in contrast to BURA, which performs well for $α$ close to $0$. As a result, we show theoretically and experimentally, that they have mutually complementary advantages.

Johannes Kraus, Raytcho Lazarov, Maria Lymbery et al.

In this paper we propose and analyze a preconditioner for a system arising from a finite element approximation of second order elliptic problems describing processes in highly het- erogeneous media. Our approach uses the technique of multilevel methods and the recently proposed preconditioner based on additive Schur complement approximation by J. Kraus (see [8]). The main results are the design and a theoretical and numerical justification of an iterative method for such problems that is robust with respect to the contrast of the media, defined as the ratio between the maximum and minimum values of the coefficient (related to the permeability/conductivity).

Beiping Duan, Bangti Jin, Raytcho Lazarov et al.

We present and analyze a space-time Petrov-Galerkin finite element method for a time-fractional diffusion equation involving a Riemann-Liouville fractional derivative of order $α\in(0,1)$ in time and zero initial data. We derive a proper weak formulation involving different solution and test spaces and show the inf-sup condition for the bilinear form and thus its well-posedness. Further, we develop a novel finite element formulation, show the well-posedness of the discrete problem, and establish error bounds in both energy and $L^2$ norms for the finite element solution. In the proof of the discrete inf-sup condition, a certain nonstandard $L^2$ stability property of the $L^2$ projection operator plays a key role. We provide extensive numerical examples to verify the convergence of the method.

Bangti Jin, Raytcho Lazarov, Zhi Zhou

We consider the initial/boundary value problem for the fractional diffusion and diffusion-wave equations involving a Caputo fractional derivative in time. We develop two "simple" fully discrete schemes based on the Galerkin finite element method in space and convolution quadrature in time with the generating function given by the implicit backward Euler method/second-order backward difference method, and establish error estimates optimal with respect to the regularity of the initial data. These two schemes are first and second-order accurate in time for nonsmooth initial data. Extensive numerical experiments for one and two-dimensional problems confirm the convergence analysis. A detailed comparison with several popular time stepping schemes is also performed. The numerical results indicate that the proposed fully discrete schemes are accurate and robust for nonsmooth data, and competitive with existing schemes.

Bangti Jin, Raytcho Lazarov, Vidar Thomée et al.

We consider three types of subdiffusion models, namely single-term, multi-term and distributed order fractional diffusion equations, for which the maximum-principle holds and which, in particular, preserve nonnegativity. Hence the solution is nonnegative for nonnegative initial data. Following earlier work on the heat equation, our purpose is to study whether this property is inherited by certain spatially semidiscrete and fully discrete piecewise linear finite element methods, including the standard Galerkin method, the lumped mass method and the finite volume element method. It is shown that, as for the heat equation, when the mass matrix is nondiagonal, nonnegativity is not preserved for small time or time-step, but may reappear after a positivity threshold. For the lumped mass method nonnegativity is preserved if and only if the triangulation in the finite element space is of Delaunay type. Numerical experiments illustrate and complement the theoretical results.

Bangti Jin, Raytcho Lazarov, Dongwoo Sheen et al.

In this work, we consider the numerical solution of an initial boundary value problem for the distributed order time fractional diffusion equation. The model arises in the mathematical modeling of ultra-slow diffusion processes observed in some physical problems, whose solution decays only logarithmically as the time $t$ tends to infinity. We develop a space semidiscrete scheme based on the standard Galerkin finite element method, and establish error estimates optimal with respect to data regularity in $L^2(D)$ and $H^1(D)$ norms for both smooth and nonsmooth initial data. Further, we propose two fully discrete schemes, based on the Laplace transform and convolution quadrature generated by the backward Euler method, respectively, and provide optimal convergence rates in the $L^2(D)$ norm, which exhibits exponential convergence and first-order convergence in time, respectively. Extensive numerical experiments are provided to verify the error estimates for both smooth and nonsmooth initial data, and to examine the asymptotic behavior of the solution.