High order fast algorithm for the Caputo fractional derivative

In the paper, we present a high order fast algorithm with almost optimum memory for the Caputo fractional derivative, which can be expressed as a convolution of $u'(t)$ with the kernel $(t_n-t)^{-\alpha}$. In the fast algorithm, the interval $[0,t_{n-1}]$ is split into nonuniform subintervals. The number of the subintervals is in the order of $\log n$ at the $n$-th time step. The fractional kernel function is approximated by a polynomial function of $K$-th degree with a uniform absolute error on each subinterval. We save $K+1$ integrals on each subinterval, which can be written as a convolution of $u'(t)$ with a polynomial base function. As compared with the direct method, the proposed fast algorithm reduces the storage requirement and computational cost from $O(n)$ to $O((K+1)\log n)$ at the $n$-th time step. We prove that the convergence rate of the fast algorithm is the same as the direct method even a high order direct method is considered. The convergence rate and efficiency of the fast algorithm are illustrated via several numerical examples.

Similar Publications

Quadratically-constrained basis pursuit has become a popular device in sparse regularization; in particular, in the context of compressed sensing. However, the majority of theoretical error estimates for this regularizer assume an a priori bound on the noise level, which is usually lacking in practice. In this paper, we develop stability and robustness estimates which remove this assumption. Read More

We provide two new methods for computing lower bounds of eigenvalues of symmetric elliptic second-order differential operators with mixed boundary conditions of Dirichlet, Neumann, and Robin type. The methods generalize ideas of Weinstein's and Kato's bounds and they are designed for a simple and straightforward implementation in the context of the standard finite element method. These lower bounds are obtained by a posteriori error estimators based on local flux reconstructions, which can be naturally utilized for adaptive mesh refinement. Read More

In this article, we consider exactly divergence-free H(div)-conforming finite element methods for time-dependent incompressible viscous flow problems. This is an extension of previous research concerning divergence-free H1-conforming methods. For the linearised Oseen case, the first semi-discrete numerical analysis for time-dependent flows is presented here whereby special emphasis is put on pressure- and Reynolds-semi-robustness. Read More

The efficient numerical integration of large-scale matrix differential equations is a topical problem in numerical analysis and of great importance in many applications. Standard numerical methods applied to such problems require an unduly amount of computing time and memory, in general. Based on a dynamical low-rank approximation of the solution, a new splitting integrator is proposed for a quite general class of stiff matrix differential equations. Read More

As already done for the matrix case for example in [Joe Harris, Algebraic Geometry - A first course, p.256] we give a parametrization of the Bouligand tangent cone of the variety of tensors of bounded TT rank. We discuss how the proof generalizes to any binary hierarchical format. Read More

Computing rational minimax approximations can be very challenging when there are singularities on or near the interval of approximation - precisely the case where rational functions outperform polynomials by a landslide. We show that far more robust algorithms than previously available can be developed by making use of rational barycentric representations whose support points are chosen in an adaptive fashion as the approximant is computed. Three variants of this barycentric strategy are all shown to be powerful: (1) a classical Remez algorithm, (2) a "AAA-Lawson" method of iteratively reweighted least-squares, and (3) a differential correction algorithm. Read More

We are concerned with the computation of the ${\mathcal L}_\infty$-norm for an ${\mathcal L}_\infty$-function of the form $H(s) = C(s) D(s)^{-1} B(s)$, where the middle factor is the inverse of a meromorphic matrix-valued function, and $C(s),\, B(s)$ are meromorphic functions mapping to short-and-fat and tall-and-skinny matrices, respectively. For instance, transfer functions of descriptor systems and delay systems fall into this family. We focus on the case where the middle factor is very large. Read More

We introduce the Harmonic Virtual Element Method (HVEM), a modification of the Virtual Element Method (VEM) for the approximation of the 2D Laplace equation using polygonal meshes. The main difference between HVEM and VEM is that in the former method only boundary degrees of freedom are employed. Such degrees of freedom suffice for the construction of a proper energy projector on (piecewise harmonic) polynomial spaces. Read More

In this article, we derive a Bayesian model to learning the sparse and low rank PARAFAC decomposition for the observed tensor with missing values via the elastic net, with property to find the true rank and sparse factor matrix which is robust to the noise. We formulate efficient block coordinate descent algorithm and admax stochastic block coordinate descent algorithm to solve it, which can be used to solve the large scale problem. To choose the appropriate rank and sparsity in PARAFAC decomposition, we will give a solution path by gradually increasing the regularization to increase the sparsity and decrease the rank. Read More

Many inverse problems involve two or more sets of variables that represent different physical quantities but are tightly coupled with each other. For example, image super-resolution requires joint estimation of image and motion parameters from noisy measurements. Exploiting this structure is key for efficiently solving large-scale problems to avoid, e. Read More