Adaptive discretization of an integro-differential equation with a weakly singular convolution kernel

Abstract

An integro-differential equation involving a convolution integral with a weakly singular kernel is considered. The kernel can be that of a fractional integral. The integro-differential equation is discretized using the discontinuous Galerkin method with piecewise constant basis functions. Sparse quadrature is introduced for the convolution term to overcome the problem with the growing amount of data that has to be stored and used in each time-step. A priori and a posteriori error estimates are proved. An adaptive strategy based on the a posteriori error estimate is developed. Finally, the precision and effectiveness of the algorithm are demonstrated in the case that the convolution is a fractional integral. This is done by comparing the numerical solutions with analytical solutions.

Introduction

Fractional order operators (integrals and derivatives) have proved to be very suitable for modeling memory effects of various materials and systems of technical interest. In particular, they are very useful when modeling viscoelastic materials, see, e.g., [3], [4], [10]. The drawback of these models is that, when the response is integrated numerically, the whole previous stress or strain history must be included in each time-step. Rather few algorithms for integrating viscoelastic responses (integral equations with singular kernels) are available. Most of them are based on the Lubich convolution quadrature for fractional order operators, see [11] and, e.g., [8]. The Lubich convolution quadrature requires uniformly distributed time-steps. This is a cumbersome restriction, in particular, when analyzing non-linear viscoelastic responses. Furthermore, it is not possible to use adaptivity and goal oriented error estimates. It also restricts the possibility to use sparse time history. In the present work we discuss these difficulties in the context of the following integro-differential equation

$$\text{u}{}_{\mathrm{t}}\text{(t)+∫}{}_0{}^{\mathrm{t}}\text{β(t−s)Au(s)}\text{d}\text{s=f(t),}\text{t∈(0,T),}$$$$\text{u(0)=u}{}_0\text{,}$$

where u_t=du/dt and A is a self-adjoint, positive definite, linear operator on a separable Hilbert space H with inner product (·,·) and norm ∥·∥. The operator A may be unbounded with domain of definition D(A)⊂H, but we assume that it has compact inverse, so that the spectral theorem applies. In applications H could typically be $\text{L}{}_2\text{(Ω)}$ for some spatial domain $\text{Ω}$, and A an elliptic partial differential operator with respect to the spatial variables, or the finite element approximation of such an operator. The abstract Hilbert space framework makes it possible to discuss time discretization without going into details about the spatial approximation. We assume that the data u₀∈H and f(t)∈H are such that the equation has a unique, appropriately regular solution.

The kernel function β may be weakly singular but integrable. More precisely, we assume that β is real-valued, belongs to L₁(0,T), and positive definite in the sense that, for any T⩾0,

$$\text{∫}{}_0{}^{\mathrm{T}}\text{∫}{}_0{}^{\mathrm{t}}\text{β(t−s)ϕ(s)ϕ(t)}\text{d}\text{s}\text{d}\text{t⩾0,}\text{∀ϕ∈L}{}_2\text{(0,T).}$$

Our chief example is

$$\text{β(t)=}\text{1}\text{Γ(α)}\text{1}\text{t}{}^{\mathrm{1-\alpha }}\text{,}\text{0<α⩽1,}$$

for which the convolution integral can be interpreted as an integral of fractional order α, see, e.g., the textbook [15], and (1.1) can be written as

$$\text{u}{}_{\mathrm{t}}\text{(t)+D}{}^{\mathrm{-\alpha }}\text{[Au](t)=f(t),}\text{t∈(0,T),}$$$$\text{u(0)=u}{}_0\text{,}$$

where the operator D^−α is the fractional order integral. In the limit α=0 the kernel β approaches the delta function and we obtain the parabolic equation

$$\text{u}{}_{\mathrm{t}}\text{+Au=f;}\text{u(0)=u}{}_0\text{,}$$

while specializing to α=1 results in

$$\text{u}{}_{\mathrm{t}}\text{+∫}{}_0{}^{\mathrm{t}}\text{Au(s)}\text{d}\text{s=f,}$$

which is equivalent to the hyperbolic equation

$$\text{u}{}_{\mathrm{tt}}\text{+Au=f}{}_{\mathrm{t}}\text{;}\text{u(0)=u}{}_0\text{,}\text{u}{}_{\mathrm{t}}\text{(0)=f(0).}$$

This means that by varying the fractional integral exponent we obtain a link between parabolic behavior, (e.g., heat conduction) and hyperbolic behavior, (e.g., wave propagation).

Here we develop an adaptive algorithm with a priori and a posteriori error estimates for solving (1.1). The a posteriori error estimate forms the basis for the adaptive strategy. For the numerical integration we adopt the discontinuous Galerkin method with piecewise constant basis functions. To overcome the problem with the growing amount of data, that has to be stored and used in each time-step, we introduce sparse quadrature for the convolution integral. Sparsely distributed time-steps are used in the distant part of the history while small steps are used in the most recent part. The idea is to break up the convolution structure by using piecewise linear interpolants between the large steps in the distant part of the history. This was first studied in [16], [18].

The present work is a further development of the study by McLean et al. [13]. Indeed, they allowed variable time-steps and introduced sparse quadrature, but their study was limited to a priori error estimates.

Similar considerations apply to fractional order differential equations, D^αu+Au=f, which can be rewritten as Volterra integral equations of the second kind by applying the integral operator D^−α, see [6]. The kernel function in the resulting integral equation is then of the same kind as in the present study. Fractional differential equations are studied in this way in [5], [6], [7]. The numerical schemes are based on the Lubich convolution quadrature and therefore need equispaced grids, or alternatively logarithmically distributed time-steps, as in [9]. This prevents the use of adaptivity and sparse quadrature. These issues will be addressed in the forthcoming paper [2] using the methods of the present work.