DISCRETE MAXIMAL REGULARITY FOR VOLTERRA EQUATIONS AND NONLOCAL TIME-STEPPING SCHEMES

. In this paper we investigate conditions for maximal regularity of Volterra equations deﬁned on the Lebesgue space of sequences (cid:96) p ( Z ) by using Bl¨unck’s theorem on the equivalence between operator-valued (cid:96) p -multipliers and the notion of R -boundedness. We show suﬃcient conditions for maximal (cid:96) p − (cid:96) q regularity of solutions of such problems solely in terms of the data. We also explain the signiﬁcance of kernel sequences in the theory of viscoelasticity, establishing a new and surprising connection with schemes of approximation of fractional models.

1. Introduction. In this paper we study a discrete time formulation for a class of integro partial differential equations with delay whose prototype equation is, where A is a closed linear operator defined on a Banach space X. This model appears in linear viscoelasticity theory, more precisely in the study of viscoelastic fluids, heat conduction with memory and electrodynamics with memory, among others [23]. In such applications the operator A is typically the negative Laplacian in X = L 2 (Ω), the elasticity operator, the Stokes operator, or the biharmonic ∆ 2 , equipped with suitable boundary conditions. Equation (1) includes many important models in PDEs. For instance, taking a(t) ≡ 1 we arrive at ∂ β t u(t) = Au(t) + g(t, u(t)), t > 0, (2) where β depends on the value of α. The existence and uniqueness of solutions to (1), or to the equivalent initial value problem a(t − s)Au(s)ds + (Qu)(t) + f (t, u(t)) where (Qu)(t) := 0 −∞ a(t − s)Aϕ(s)ds has been studied extensively in recent years. Numerical methods for time discretization of Volterra type equations (3) have been proposed by various authors. Recent developments on optimal p − q timespace estimates for the corresponding linearized equations ask for time discretizations that preserve maximal regularity [4,11,12,15,18]. Indeed, it has been proved that this kind of regularity results in a discrete setting provide more precise error estimates for the numerical analysis of nonlinear PDEs [11,13,17]. This motivates the following first question: (Q1) Given a closed linear operator A and a kernel a, for which time discretization scheme of (3) there is maximal regularity in the p -setting?.
On the other hand, it is well known in linear viscoelasticity theory that the kernel a(t) describes the behavior of the material under compression, which is a combination of two symbols called spring and dashpot in mechanical engineering. Some standards models are: Hookean solid, which is the prototype of spring; Newtonian fluid, being the prototype of dashpot; Kelvin-Voigt solid, that can be viewed as a spring and a dashpot in parallel; Maxwell fluid, which corresponds to a spring and a dashpot in series; Poynting Thompson solid, corresponding to Hookean solid and Maxell fluid in parallel; and power type materials, that represent an infinite series of springs and dashpots [8]. For instance, the kernel a(t) := g α (t) := t α−1 Γ(α) is the power type material function [23, pp.130.131] which coincides with a Newtonian fluid (dashpot) when α = 1 and a Hookean solid (spring) when α = 2.
Our first key observation in this paper is that standard numerical time-stepping schemes of approximation for fractional evolution equations have the form [11]: where τ is the step size, b is a scalar-valued sequence and the star denotes finite convolution.
Examples of schemes that fit into this framework are the Backward Euler scheme, the second order backward difference formula and the L1-scheme, among others Lubich's quadrature methods. See [11] for more details.
In such a case, the discrete model where f is given (but in concrete examples it corresponds to a time discretization of the nonlinear term g(·, u) in (2)) is equivalent to the abstract Volterra equation for the explicit sequence kernel b τ (n) = τ α k α (n) and where, in relation to (4), A novel characteristic of our model (5) is that we will assume zero-padding in the past time, i.e. u(n) = 0 for all n = −1, −2, ... This is a concept used in digital signal processing and refers to adding zeros to a time domain signal to increase its length. In physical terms, it refers to the notion of causality. It is remarkable that this hypothesis for (5) coincides with the fact that Lubich's convolution quadrature method implicitly considers zero-padding in the negative real axis, see [22, p.131, lines 1-2].
A second important observation is that by means of the recently defined Poisson transformation [2,19], which represents a way of sampling time continuous systems into discrete time models, we have the remarkable relation [19,Example 3.3] It reveals a new and surprising association between (5) and the Volterra equation This relationship shows that time-stepping schemes and material functions are correlated, which in terms of the above analysis can be stated as follows: Power type material function ⇐⇒ Fractional backward Euler scheme. To be more precise, the above expression means that whereas in continuous time the kernel g α (t) represents the power material function of Volterra equations in viscoelasticity theory [23, p.131 (vi)], we have that in discrete time the sequence kernel k −α (n) is in correspondence to the characteristic function of the fractional backward Euler scheme by means of its generating function, i.e. we have, see [11, identities (3.2)]. In particular, when α = 1 we have that δ(ξ) = (1 − ξ) is the characteristic function of the Backward Euler scheme method [11,Section 3.1] and, on the other hand, the function g 1 (t) ≡ 1 is the kernel that represents a Newtonian fluid for Volterra equations in viscoelasticity theory [23, p.130 (ii)]. This correspondence is expressed as: Newtonian fluid ⇐⇒ Backward Euler scheme. In this way, we reveal new insights about a surprising relation between time-stepping schemes and linear viscoelasticity theory that, as far as we know, has not been previously observed in the existing literature. Our findings are also supported and coincident with very recent research in mechanical engineering, where methodologies to identify the desired viscoelastic behavior before choosing a real material has been developed [8]. In terms of time schemes, the connection stated in this paper can be interpreted as a methodology to identify the desired (and probably best) time stepping scheme in terms not only of the mathematical model but also of the characteristics of the real material that it models.
In view of the above considerations, we address the following second question: (Q2) Which class of material type functions are time-stepping schemes correlated with?.
Whereas question (Q2) is completely new, question (Q1) has been answered in many papers from different perspectives. For example, for the model (4) when α = 1, Ashyralyev, Piskarev and Weis [6] showed the discrete maximal regularity using the implicit Euler method; Kovács, Li and Lubich [15] proved the discrete maximal regularity for the Crank-Nicholson, BDF and A-stable Runge-Kutta methods; Kemmochi and Saito [12,13] proved the maximal p -regularity for the θ-method, and Leykekhman and Vexler [16] proved the maximal regularity for Galerkin methods. In the fractional case, the first work is due to Lizama [20]. Recently, a further important step was given by Jin, Li and Zhou [11]. They have provided an analysis for several time-stepping schemes, including the convolution quadratures generated by the implicit Euler method and second-order backward difference formula, the L1 scheme, the explicit Euler method and a fractional variant of the Crank-Nicolson method.
The paper is organized as follows: In Section 2, we introduce some basic concepts and results in the existing literature related to discrete time Fourier transforms and operator-valued multiplier theorems that will be later needed. In Section 3, we answer question (Q1). We introduce the new concept of 1-regular sequences that together with the R-boundedness of the symbol of the equation ensures maximal regularity of the (linearized) discretization scheme (3) with zero-padding in the past time, in the p -setting. Using this new concept of 1-regularity, our first main contribution in this paper is to unify earlier approaches to (Q1) in a simpler way. More concretely, we prove the following theorem: Theorem. Let X be a U M D space, 1 < p < ∞, and b : Z → C such that b(n) = 0 for all n ∈ Z − . Suppose that b is 1-regular and satisfies that the z-transform of b converges on the complex unit disk and their radial limit exists for all t ∈ T 0 := (−π, π) \ {0}. Moreover, b(t) = 0 for all t ∈ T 0 . Assume that Then the solution of the equation defined as 0 by negative values of Z, satisfies the following maximal regularity estimate u p (Z;X) + b Au p (Z;X) ≤ C f p (Z;X) .
We also obtain an easier computable condition for the Hilbert space case, where R-boundedness can be replaced by norm boundedness. Our strategy of proof follows closely the recent works [11,15,19,21] and employs the (discrete) operator-valued Fourier multiplier technique developed by Blunck [3,7].
In Section 3, we show that our model (6) includes an important class of schemes of approximation given by the convolution quadrature method introduced by C. Lubich in 1988 [22], by identifying kernels b(n) with characteristic functions of each scheme by means of the discrete time Fourier transform. For example, we obtain the Explicit and Implicit Euler scheme, the Backward Euler and the Fractional Crank-Nicholson scheme among others and, as an application of our results, we provide maximal regularity estimates for the aforementioned time-stepping schemes. Moreover, we address question (Q2) and we obtain a direct correlation between time-stepping schemes and material type functions in viscoelasticity theory.
2. Analytical framework and notation. In this section, we present some notations and preliminary results that will be needed throughout this work. Except for Lemma 2.3 below, all the quoted results are contained in the existing literature on the subject.
Let X be a Banach space. We denote by S(Z; X) the space of vector-valued sequences f : Z → X such that for each k ∈ N 0 there exists a constant C k > 0 satisfying p k (f ) := sup n∈Z |n| k f (n) < C k and if X = R we denote the space as S(Z). We write as C n per (R; X), n ∈ N 0 , the space of all 2π-periodic X-valued and n-times continuously differentiable functions defined in R. In what follows, we will denote T := (−π, π) and T 0 : and it follows that T f ∈ S (Z, X) = {T : S(Z) → X : T is linear and continuous}.
Remark 1. This mapping identifies p (Z; X) with a subspace of S (Z; X). Thus, a function f ∈ p (Z; X) will be identified with T f ∈ S (Z, X).
There exists another natural mapping that identifies C ∞ per (T; X) with a subspace of D (T; X) = {T : C ∞ per (T) → X : T is linear and continuous} which assigns to each S ∈ C ∞ per (T; X) the linear map and we have L S ∈ D (T; X).
and the inverse transform is given by where ϕ ∈ C ∞ per (T; X). The discrete time Fourier transform (DTFT) between the spaces of distributions S (Z; X) and D (T; X) then follows as: whose inverse F −1 : D (T; X) → S (Z; X) is given by We finally present a technical lemma which will be necessary througout the paper. We first need the following definition.
Definition 2.2. Given u : Z → X and v : Z → C the convolution product between u and v is defined as whenever the series converges. Moreover, the convolution of a distribution T ∈ S (Z, X) with a sequence a : N 0 → C is defined by where Observe that a • ϕ ∈ S(Z). We denote the Heaviside sequence. For a : the finite convolution. We will need the following Lemma that corresponds to a weaker version of [21, Lemma 2.2].
Lemma 2.3. Let u, v : Z → X be given and a : N 0 → C which is defined by 0 for negative values of n. We assume that the series converges on the complex unit disk, and that the radial limitâ(t) = lim r 1ã Proof. Given ϕ ∈ C ∞ per (T) and by hypothesis, we obtain for all n ∈ Z. Hence Choosing ϕ k (t) := e −ikt , k ∈ Z we find that proving the Lemma.
We recall the notion of R-bounded sets and p -multipliers in the space B(X, Y ) of bounded linear operators from X into Y endowed with the uniform operator topology.
for all See [3, Section 2.2] and [9] for more information about R-boundedness. We next recall the concept of p -multiplier. (15) for all f ∈ p (Z; X) and all ϕ ∈ C ∞ per (T). Here We now recall the following Fourier multiplier theorem for operator-valued symbols given by S. Blunck [7,3]. This theorem provides sufficient conditions to ensure when an operator-valued symbol is a multiplier, and allows to establish an equivalence between p -multipliers and the notion of R-boundedness for the U M D class of Banach spaces. For more information about these spaces see [ The converse of Blunck's theorem also holds for any Banach spaces X, Y as follows: 3. Abstract setting: A characterization of maximal p -regularity. Let b : N 0 → C be given and X be a Banach space. We denote: For a given vector-valued sequence f : N 0 → X we consider the abstract discrete equation where A is a closed linear operator with domain D(A) defined in a Banach space X. Observe that using the notation of the finite convolution we have u = b Au + f = Au b + f . We also recall that by [D(A)] we denote the domain of A endowed with the graph norm.
For a sequence b : N 0 → C extended to negative subscripts n by 0, we assume the following hypothesis: (H) b : The z-transformb of b converges on the complex unit disk: We introduce the following definition. Observe that, in some sense, it corresponds to the discrete counterpart of the notion of k-regularity introduced in the paper [14]. See also [23].
Remark 2. It follows immediately from the definition that ifb(t) is 1-regular and there exists a sequence a such thatâ(t)b(t) = 1 then a is 1-regular, too. In the following result, we show the equivalence between the R-boundedness of the operator valued symbol of the equation (16) defined by (1 −b(t)A) −1 and the fact that it is an p -multiplier.
and denote M (t) := (1 −b(t)A) −1 , then the following assertions are equivalent: Therefore, (i) =⇒ (ii) By hypothesis we have that there exists a bounded operator T such that (15) holds. Now, (ii) holds as a consequence of Theorem 2.7.
Proof. We observe that N (t) =b(t)M (t) and AN (t) = M (t) − I. Therefore and using (18) we get and since b is 1-regular, we conclude from Theorem 2.6 the assertion.
The following theorem shows sufficient conditions for maximal p 0 -regularity. It corresponds to the main abstract result of this paper.
Then equation (16) has maximal p 0 -regularity. In addition, u, b Au ∈ p (Z; X) and there exists a constant C > 0 ( independent of f ∈ p (Z; X)) such that the following maximal regularity estimate holds u p (Z;X) + b Au p (Z;X) ≤ C f p (Z;X) .
Proof. Let f ∈ p 0 (N 0 ; X) be given. By hypothesis, there exists u ∈ p (Z; [D(A)]) such that for all ϕ ∈ C ∞ per (T). In particular, the closedness of A yields Define N (t) =b(t)M (t) and note that AN (t) = M (t) − I. By Theorem 3.5 there for all ψ ∈ C ∞ per (T) where Observe that by hypothesis ϕ(t) := ψ(t) b(−t) ∈ C ∞ per (T). Setting ϕ in (21), from (22) and (23) we get From Lemma 2.3 we conclude from the above identity and the fact that b(n) = 0 for negative values of n that that is v ∈ p 0 (N 0 ; X). Since AN (t) = M (t) − I, after multiplication by e int ϕ(t) and integration over the interval (−π, π), we have By replacing (22) and (20) in the above identity and taking into account (24) we obtain for all ϕ ∈ C ∞ per (T). Choosing ϕ(t) = e −ikt for all k ∈ Z, we conclude that u ∈ p 0 (N 0 ; [D(A)]) and satisfies the equation (16). We have proven the existence of a solution. It remains to show the solution is unique.
Let u ∈ p 0 (N 0 ; [D(A)]) be a solution of (16) with f ≡ 0. Since b is 1-regular and satisfies the hypothesis (H) b , we obtain the following identities: valid for any S ∈ C ∞ per (T 0 , B(X, Y )). For all ϕ ∈ C ∞ per (T), using the hypothesis, the fact that M ∈ C ∞ per (T 0 , B(X, where in the last equality we have used (25) with S = ϕ · M − . Therefore using (10) we get

[D(A)]), and the identity
Taking ϕ k (t) := e −ikt , k ∈ Z we obtain u ≡ 0 and then the claim is proven. Finally, the last assertion follows from the closed graph theorem. Indeed, for each f ∈ p 0 (N 0 ; X) we have already proved the existence of a unique sequence u : N 0 → [D(A)] such that u ∈ p 0 (N 0 ; [D(A)]) and it satisfies: We now define the operator T : p 0 (N 0 ; X) → p 0 (N 0 ; [D(A)]) as T (f ) = u. Using the closed graph theorem it follows that T is a bounded operator and moreover, there exists a constant C > 0 such that As a immediate consequence of the fact that R-boundedness is equivalent to boundedness in Hilbert spaces we obtain the following corollary.

Applications to schemes of approximation.
The following examples of sequences are related to several classes of discrete operators that correspond to time discretizations of fractional evolution equations. This kind of nonlocal operators has experimented a rapid expansion over the past few years and provides a new paradigm in modelling partial differential equations. Applications range from subdiffusions processes observed in molecular physics, anomalous transport in cell biology, to long range dynamics from nano-scale to geophysical systems. In linear viscoelasticity theory, fractional models are frequently used to describe the behavior of polymeric materials [23, p.132].
We start with the observation that standard numerical time-steeping schemes of approximation for fractional evolution equations have the form where τ is the step size, b is a scalar-valued sequence and the star denotes finite convolution, see [11] and the examples below.
We denote by C the following set where δ i,j denotes the Kronecker delta, i.e.
Consider the abstract fractional evolution equation where A is a closed linear operator defined in a Banach space X, and the following scheme for its discretization where f (n) is a time discretization of F (t), e.g. f (n) := F (τ n). See also [11, scheme (6.1)]. Examples of schemes that fit into this framework are the Backward Euler scheme, the second order backward difference formula and the L1-scheme. See [11] for more details. Suppose that b ∈ C. Then, convolving (29) by a we obtain the equivalent system where g(n) := τ α (a f )(n). Consequently, the abstract setting of Section 3 applies.
It is important to observe that in contrast with previous papers (see e.g. [11] and [15]), where the discretization operator -i.e. the coefficients of the methodare implicitly given by a characteristic function of the so-called generating formula in this work we are able to find explicitly a sequence b ∈ C for some approximation schemes. In other words, we prove that for some approximation schemes of abstract fractional evolution equations, the condition holds. This is possible thanks to the knowledge of a crucial sequence which is defined as follows: For any α ∈ R \ {0, −1, −2, ...} we define the sequence where k 0 (n) := δ 0,n . This sequence has been considered in several papers in connection with fractional differences of order α [1,2,19,20,21]. From [1, formula (2.1)] we have the generation formula Among others, the following important property holds: k α * k β = k α+β [19,2]. Our approach reveals a new connection between schemes of approximation, discrete Volterra equations and material functions in linear viscoelasticity theory. We recall from [23, pp.130-131] that standards models are: Hookean solid, Newtonian fluid, Kelvin-Voigt solid, Maxwell fluid, Poynting Thompson solid and power type materials. In this paper, we point out that different schemes defined by means of kernels b are in direct relation with different materials and consequently it could give us new insights about the appropriate scheme to be used in concrete applications. We note that this interesting information is concentrated in the kernel, i.e. a function a(t), of Volterra equations. We will reveal the connection between the kernels a(t) and b(n) in the following examples.
is the characteristic function of the backward Euler scheme (BE) in time t, with a constant time step size τ > 0. It has been used for the discretization of the abstract fractional evolution equation recently studied in [11, Section 3.1]. Comparing (30) with (32) and using (31) we Then the radial limit a(t) = τ α (1 − e −it ) −α exists for all t ∈ T. We can easily check that Therefore the sequence a(n) is 1-regular. We observe that (a * b)(n) = (k α * k −α )(n) = k 0 (n) = δ 0,n which means that b ∈ C, i.e. the condition (H) a holds. For example, in case α = 1 we have a(n) ≡ τ. This means that for u(0) given, the difference equation is equivalent to the abstract Volterra equation In the general case of α > 0, the difference-differential equation , is equivalent to the abstract Volterra equation which is the discrete version of the Volterra equation We observe that the kernel a(t) = g α (t) is a power type material function in linear viscoelasticity theory [23, p.131 (vi)]. The correspondence with the kernel k α (n) is justified by the identity valid for α > 0. See [19]. In other words, the sequences k α (n) are in direct correspondence, by means of the Poisson transformation [2,19], with the functions g α (t) that have the meaning of the shear modulus of a homogeneous isotropic linear material in viscoelasticity theory [ defines the fractional second-order backward difference scheme (BDF2) for discretizing the model 33. See [11,Section 3.2]. Comparing (30) with (35) and using (31) and the identity In the case α = 1, we consider the sequence a(n) defined by We can easily check thatâ(t) = , t ∈ T.
On the other hand, from [23, Chapter I, Section 5.2 (iv)] it is known that the kernel a(t) = 1 − e −at is the standard material function for a Maxwell fluid in linear viscoelasticity theory. We note the following identity: It shows that the sequence (36) is linked by means of the Poisson transformation with the material function a(t) = 1 − e −2t . Observe that in viscoelasticity theory we have a = µ/ν where µ, ν have physical meanings. In particular, it should be interesting to have a more closed view between the dependence of the second order difference scheme and these two numbers. We summarize this example in the following correspondence Second-order backward difference scheme ⇐⇒ Maxwell fluid.
In particular, with a = 2 we get ∞ 0 p n (t)g α (t)e −2t dt = k α (n) 3 n+α and with a = −1/2 we have ∞ 0 p n (t)g α (t)e − 1 2 t dt = 2 α k α (n). Therefore, the convolution property of the Poisson transform [19] produces This identity shows a correspondence between a α (n) and the kernel which, to our knowledge, cannot be clearly identified with any material function in viscoelasticity theory. However, this is in some sense similar to a material function of the form Power type × Maxell fluid.
defines the explicit Euler difference scheme as it has been studied in [11,Section 5] for the discretization of the fractional equation (33). Define the sequence Using the generation formula (31) we obtain Note that which shows 1-regularity of the sequence a(n) for all α > 0. Also, we have δ(t) = 1 which means that condition (H) a is verified. Compared with Example 4.1, the Volterra model has in this example the slightly different form but concerning the comparison with viscoelasticity theory, the same correlation than in Example 4.1 applies. This example shows that same material functions can define somewhat different approximation schemes.
Example 4.4. (Fractional Crank-Nicholson scheme) Let 0 < α < 2 be given. The symbol [11,Section 6]. We define the sequence For α = 1 it corresponds to the Crank-Nicholson scheme with step τ , that is δ(t) = 1−ξ = − 1 2 + 1 1−ξ shows that a(n) = − τ 2 δ 0,n + τ, n ∈ N 0 , and 0 otherwise. Therefore, this scheme behaves essentially as the backward Euler scheme except by the gap in n = 0. Hence, it should correspond to a Newtonian fluid. In case α > 0 the description of the sequence kernel in (39) shows the same behavior as a power type material function. Compare it with Example 4.1.
Using the generation formula (31) we obtain after a simple computation Note that [ a(t)] a(t) which shows 1-regularity of the sequence a(n) for all 0 < α < 2. Moreover, it is clear that the condition (H) is satisfied. In fact, using the Cauchy product, we find for each n ∈ N 0 and defined by 0 otherwise. We summarize this example with the following correspondences Crank-Nicholson scheme ⇐⇒ Newtonian fluid and Fractional Crank-Nicholson scheme ⇐⇒ Power type material function.
We remark that the above correspondence of the Crank-Nicholson scheme, that shows coincidence with the Euler scheme, is not at all surprising as it was also observed by Jin, Li and Zhen [11, Section 1] but in the context of p -maximal regularity.
In this section we give some applications of Theorem 3.6 that are connected with the different schemes of approximation considered so far. They also improve analogous results of earlier papers on this subject; cf. [11].
Example 5.1. We verify the conditions provided in Theorem 3.6 in order to show the existence and uniqueness of p 0 (N 0 , L 2 (R)) solutions for the following equation where we set in equation (16) b(n) = 1 2 n , A = ∂ ∂x with domain D(A) = W 1,2 (R), and f ∈ p 0 (N 0 , L 2 (R)). Observe by Remark 3.3 that b is k-regular for all k ∈ N and b(t) = 2(2 − e −it ) −1 = 0 for all t ∈ T 0 . See also [10] where this kernel b has been previously considered. Moreover, since the operator Au = u with domain D(A) = W 1,2 (R) generates a contraction C 0 -group on L 2 (R) there exists M > 0 such that the following estimate holds: , for all z such that (z) > 0.
Remark 3. If A is an R-sectorial operator on X of angle απ/2 (where 0 < α < 2) then the hypothesis on A of the above corollary are automatically satisfied.
Remark 4. From Example (4.3) we have that the explicit Euler method represented by the sequence kernel a(n) = τ α k α (n−1) satisfies all the hypothesis of Theorem 3.6. Consequently, an analogous corollary about maximal p 0 -regularity for the explicit Euler scheme still holds.
The following theorem is correllated with [11,Theorem 6]. It makes use of Example 4.2. 3 n−l k α (n − l)A(k α u)(l) + f (n), n ∈ N 0 , has a unique solution u ∈ p 0 (N 0 ; X) that satisfies the following maximal p 0 -regularity estimate where the constant C is independent of f, τ, α and A.