Guaranteed error bounds for stabilised space-time IgA approximations to parabolic problems

# Guaranteed error bounds for stabilised space-time IgA approximations to parabolic problems

## Abstract

The paper is concerned with space-time IgA approximations of parabolic initial-boundary value problems. We deduce guaranteed and fully computable error bounds adapted to special features of IgA approximations and investigate their applicability. The derivation method is based on the analysis of respective integral identities and purely functional arguments. Therefore, the estimates do not contain mesh-dependent constants and are valid for any approximation from the admissible (energy) class. In particular, they provide computable error bounds for norms associated with stabilised space–time IgA approximations as well as imply efficient error indicators enhancing the performance of fully adaptive solvers. The last section of the paper contains a series of numerical examples where approximate solutions are recovered by IgA techniques. The mesh refinement algorithm is governed by a local error indicator generated by the error majorant. Numerical results discussed in the last section illustrate both reliability, as well as the quantitative efficiency of the error estimates presented.

## 1 Introduction

Time-dependent systems governed by parabolic partial differential equations (PDEs) are typical models in scientific and engineering applications. This triggers their active investigation in modelling, mathematical analysis and numerical solution. By virtue of fast development of parallel computers, treating time as yet another dimension in space in evolutionary equations became quite natural. The space-time approach is not affected by the disadvantages of time-marching schemes. Its various versions can be useful in combination with parallelisation methods, e.g., those discussed in [12, 13, 30].

Investigation of effective adaptive refinement methods is crucial for constructing fast and efficient solvers for PDEs. At the same time, scheme localisation is strongly linked with reliable and quantitatively efficient a posteriori error estimation tools. These tools are intended to identify the areas of a computational domain with relatively high discretisation errors and by that provide a fully automated refinement strategy in order to reach the desired accuracy level for the current approximation. Local refinement tools in IgA such as T-splines, THB-splines, and LR-splines were combined with various a posteriori error estimation techniques, e.g., error estimates using hierarchical bases [9, 41], residual-based [19, 42, 5, 24, 14], and goal-oriented error estimates [40, 8, 25, 26]. Below we use a different (functional) method providing fully guaranteed error estimates in various weighted norms equivalent to the global energy norm. These estimates include only global constants (independent of the mesh characteristic ) and are valid for any approximation from the admissible functional space. Functional type error estimates (so-called majorants and minorants of deviation from the exact solution) were introduced in [37, 38] and later applied to different mathematical models [34, 31]. They provide guaranteed, sharp, and fully computable upper and lower bounds of errors. This approach, in combination with the IgA approximations generated by tensor-product splines, was proposed and investigated in [23] for elliptic boundary value problems (BVP).

In this paper, we derive new functional type a posteriori error estimates for time-dependent problems (cf. also [36]) in the context of the stabilised space-time IgA scheme introduced in [30]. This scheme exploits the time-upwind test function based on the space-time streamline diffusion method (see, e.g., [17, 20, 21]) and the approximations provided by the IgA framework. The obtained functional estimates, in turn, provide bounds for the error measured in new stronger norm induced by alternative stabilised space-time variational formulation of the parabolic problem. By exploiting the universality and efficiency of considered error estimates as well as the smoothness of IgA approximations, we aim to construct fully-adaptive, fast and efficient parallel space-time methods that could tackle complicated problems inspired by industrial applications. We also study the numerical properties of newly derived error bounds and compare their performance to the behaviour of the bounds known from [36] on an extensive set of examples.

The outline of the paper is as follows: Section 2 states the problem, discusses its solvability and provides an overview of the existing error control tools for initial BVPs (I–BVPs). In Section 3, we deduce new functional type a posteriori error estimates using a stabilised formulation of parabolic I–BVPs. Our analysis is based on a series of transformations performed on a stabilised variational setting; the result of these transformations defines respective generalised solutions. Section 4 presents a stabilised space-time IgA scheme with its main properties along with an overview of main ideas and definitions of the IgA framework. Section 5 is devoted to the algorithmic realisation of an adaptive procedure based on the a posteriori error estimates discussed above. Finally, in Section 6 we present and discuss obtained numerical results that demonstrate the efficiency of several majorants and the error identity for a comprehensive range of examples.

## 2 Parabolic model problem

Let , , denote a space-time cylinder, where , , is a bounded Lipschitz domain with a boundary , and is a given time interval with the final time , . Here, the boundary of the space-time cylinder is defined as with , , and . We discuss our approach to guaranteed error control of space-time approximations within the paradigm of a classical linear parabolic I–BVP: find satisfying the parabolic PDE, the boundary condition, and the initial condition

 ∂tu−Δxu=finQ,u=0onΣ,andu=u0on¯¯¯¯Σ0, (1)

respectively, where is a time derivative, denotes the Laplace operator in space, is a given source function, and is prescribed initial data. Here, is a space of square-integrable functions over equipped with the usual norm and the scalar product denoted respectively by

 Extra open brace or missing close brace

By , , we denote standard Sobolev spaces of functions having derivatives of the order with respect to (w.r.t.) space and time. Next, we introduce the Sobolev spaces

 H1,00(Q) :={u∈L2(Q) :∇xu∈[L2(Ω)]d,u∣∣Σ=0}, (2) V10:=H10(Q) :={u∈H1(Q) :u∣∣Σ=0}, H10,¯0(Q) :={u∈V10 :u∣∣ΣT=0}, V10,0–:=H10,0–(Q) :={u∈V10 :u∣∣Σ0=0}, VΔx0:=HΔx0(Q) :={u∈V10 :Δxu∈L2(Q)}, VΔx0,0–:=HΔx0,0–(Q) :={u∈VΔx0 :u∣∣Σ0=0},with the norm∥w∥2VΔx0,0–:=∥Δxw∥2Q+∥∂tw∥2Q, V∇x∂t0,0–=H∇x∂t0,0–(Q) :={w∈VΔx0,0– :∇x∂tw∈L2(Q)}.

For the vector-valued functions, we additionally introduce the Hilbert spaces

 Hdivx,0(Q) :={y∈[L2(Q)]d:divxy∈L2(Q)}and Hdivx,1(Q) :={y∈Hdivx,0(Q):∂ty∈[L2(Q)]d}.

It follows from [28] that, if and , then problem (1) is uniquely solvable in , and the solution depends continuously on in the -norm. Moreover, according to [28, Remark 2.2], the norm is an absolutely continuous function of for any . If , then the problem is uniquely solvable in a wider class , and meets the modified variational formulation

 (∇xu,∇xw)Q−(u,∂tw)Q=:a(u,w)=l(w):=(f,w)Q+(u0,w)Σ0 (3)

for all , where . According to well-established arguments (see [27, 28, 43, 6, 7]), without loss of generality, we can ‘homogenise’ the problem, i.e., consider (3) with .

Our main goal is to establish fully computable estimates for space-time approximations of this class of problems. For this purpose, we use a functional approach to derive a posteriori error estimates. The first and the simplest forms of such estimates have been derived in [36] for (1). The paper [36] provides the upper bound of the norm

 |||e|||2(ν1,ν2):=ν1∥∇xe∥2Q+ν2∥e∥2ΣT,νi≥0, (4)

where is the difference between the exact solution and any approximation in the respective energy class . Assuming for simplicity that the initial condition is satisfied exactly, it is shown that for any approximating and any , we have the following inequality:

 |||e|||2(2−ν,1) :=(2−ν)∥e∥2Q+∥e∥2ΣT ≤¯¯¯¯¯MI(v,y;β):=1ν((1+β)∥y−∇xv∥2Q+(1+1β)C2F∥divxy+f−∂tv∥2Q), (5)

where is an auxiliary parameter, and stands for the constant in the Friedrichs inequality [11]

 ∥w∥Q≤CF∥∇xw∥Q,∀w∈H1,00(Q). (6)

The numerical properties of (5) w.r.t. the time-marching and space-time methods are discussed in [33, 32, 18] in the framework of finite-difference and finite-element schemes. The advanced upper error-bound (valid for the same error norm) introduced in [36] contains an additional auxiliary function . For the same and , as well as any , any fixed parameters and , an alternative majorant has the form

 (2−ν)∥∇xe∥2Q+(1 −1γ)∥e∥ΣT≤¯¯¯¯¯MII(v,y,η):=γ∥η∥2ΣT+∥u0−v∥2Σ0−2(η,u0−uh)Σ0+2F(v,η) +1ν{(1+β)∥y−∇xv+∇xη∥2Q+C2F(1+1β)∥divxy+f−∂tv−∂tw∥2Q}, (7)

where

 F(v,η):=(∇xv,∇xη)+(∂tv−f,η).

The optimal parameters , , and are easily defined in each particular case by minimisation of the respective majorant.

Finally, we note that for the case where , the heat equation (1) imposes the error identity (see [1]):

 ∥Δxe∥2Q+∥∂te∥2Q+∥∇xe∥2ΣT=:|||e|||2L≡EId2(v):=∥∇x(u0−v)∥2Σ0+∥Δxv+f−∂tv∥2Q. (8)

The numerical performance of estimates and , and of the identity (8) is studied in Section 6.

## 3 Error majorants

In this section, we derive error majorants of the functional type for a stabilised weak formulation of parabolic I–BVPs. They provide guaranteed and fully computable upper bounds of the distance between the exact solution and some approximation . Functional nature of the majorants allows us to obtain a posteriori error estimates for any conforming approximation .

We begin by testing (1) with the time-upwind test function

 λw+μ∂tw,w∈V∇x∂t0,0–,λ,μ≥0, (9)

and arrive at the stabilised weak formulation for , i.e.,

 (∂tu,λw+μ∂tw)Q+(∇xu,∇x(λw+μ∂tw))Q=:as(u,w)=ls(w):=(f,λw+μ∂tw)Q,∀w∈V∇x∂t0,0–. (10)

Then, the error , ( this condition is required to ensure the existence of the term ), is measured in terms of the norm generated by the bilinear form , i.e.,

 |||e|||2s,νi:=ν1∥∇xe∥2Q+ν2∥∂te∥2Q+ν3∥∇xe∥2ΣT+ν4∥e∥2ΣT, (11)

where are the positive weights introduced in the derivation process.

To obtain guaranteed error bounds of , we apply the method similar to the one developed in [36, 33] for parabolic I–BVPs. As a starting point, we consider the space of functions (see (2)) equipped with the norm

 ∥w∥2V∇x∂t0,0–:=supt∈[0,T]∥∇xw(⋅,t)∥2Ω+∥w∥2VΔx0,0–,

where , which is dense in . According to [28, Remark 2.2], the norms and are equivalent. Below, we exploit the density this property to derive the majorants (11) in Theorems 1 and 2.

###### Theorem 1

For any and , the following estimate holds:

 (2−1γ)(λ∥∇xe∥2Q +μ∥∂te∥2Q)+λ∥e∥2ΣT+μ∥∇xe∥2ΣT =:|||e|||I,2s≤¯¯¯¯¯MIs,h(v,y;γ,β,α):=γ{λ¯¯¯¯¯MI(v,y;β)+μ˜MI(v,y;α)}, (12)

where is the majorant defined in (5) with , i.e.,

 ¯¯¯¯¯MI(v,y;β):=(1+β)∥rId∥2Q+(1+1β)C2F∥rIeq∥2Q

and

 ˜MI(v,y;α):=(1+α)∥divxrId∥2Q+(1+1α)∥rIeq∥2Q.

Here, is the Friedrichs constant (6), and are the residuals defined by relations

 rIeq(v,y) :=f−∂tv+divxyandrId(v,y):=y−∇xv, (13)

are the weights introduced in (9), , and

Proof: Let be a sequence in such that in for . It satisfies the identity (cf. (10))

 as(un,w)=(fn,λw+μ∂tw)Q,wherefn=(un)t−Δxun∈L2(Q). (14)

Next, we consider a sequence approximating in the sence that in for . By subtracting from the left- and right-hand side (LHS and RHS, respectively) of (14) and by setting the test function to be the following difference , we arrive at the error identity

 λ∥∇xen∥2Q+μ∥∂ten∥2Q+12(μ∥∇xen∥2ΣT+λ∥en∥2ΣT)=λ((fn−∂tvn,en)Q−(∇xvn,∇xen)Q)+μ((fn−∂tvn,∂ten)Q−(∇xvn,∇x∂ten)Q). (15)

First, we modify the RHS of (15) by means of the relation

 (divxy,λen+μ∂ten)Q+(y,∇x(λen+μ∂ten))Q=0.

The obtained result can be presented as follows:

 λ∥∇xen∥2Q+μ∥∂ten∥2Q +12(μ∥∇xen∥2ΣT+λ∥en∥2ΣT) =λ((fn−∂tvn+divxy,en)Q+(y−∇xvn,∇xen)Q) +μ((fn−∂tvn+divxy,∂ten)Q+(y−∇xv,∇x∂ten)Q) =λ((rIeq(vn,y),en)Q+(rId(vn,y),∇xen)Q) +μ((rIeq(vn,y),∂ten)Q+(rId(vn,y),∇x∂ten)Q). (16)

Integrating by parts in the term leads to the identity

 μ(rId,∇x(∂ten))Q=−μ(divx(y−∇xvn),∂ten)Q=−μ(divxy−Δxvn,∂ten)Q.

Using density arguments, i.e., , in , and in , as , we arrive at the identity formulated for with , i.e.,

 λ∥∇xe∥2Q+μ∥∂te∥2Q +12(μ∥∇xe∥2ΣT+λ∥e∥2ΣT) Missing \left or extra \right (17)

By means of the Hölder, Friedrichs, and Young inequalities with positive scalar-valued parameters , , and , we deduce estimate (12).

The next theorem assumes higher regularity on the approximation and the auxiliary function .

###### Theorem 2

For any and , we have the inequality

 (2−1ζ)(λ∥∇xe∥2Q +μ∥∂te∥2Q)+μ(1−1ϵ)∥∇xe∥2ΣT+λ∥e∥2ΣT=:|||e|||II,2s≤¯¯¯¯¯MIIs,h(v,y;ζ,α,ϵ,β) :=ϵμ∥rId∥2ΣT+ζ(λ((1+α)¯¯¯¯¯MI(v,y;β)+(1+1α)μ2λ2∥∂trId∥2Q)+μ∥rIeq∥2Q),

where is the majorant defined in (5) with , is the Friedrichs constant in (6), and are the residuals defined in (13), are the parameters from (9), , , and .

## 4 Stabilised formulation and its IgA discretisation

For the reader convenience, we recall the general concept of the IgA approach, the definitions of B-splines (NURBS), and their use in geometrical representation of the space-time cylinder as well as in the construction of the IgA trial spaces, which are used to approximate the solution of the variational problem (3).

Throughout the paper, denotes the degree of polynomials used for the IgA approximations, whereas denotes the number of basis functions used to construct a -spline curve. A knot-vector is a non-decreasing set of coordinates in the parameter domain, written as , , where and . The knots can be repeated, and multiplicity of the -th knot is indicated by . In what follows, we consider only open knot vectors, i.e., . For the one-dimensional parametric domain , denotes a locally quasi-uniform mesh, where each element is constructed by distinct neighbouring knots. The global size of is denoted by , where

The univariate B-spline basis functions are defined by means of Cox-de Boor recursion formula

 ˆBi,p(ξ):=ξ−ξiξi+p−ξiˆBi,p−1(ξ)+ξi+p+1−ξξi+p+1−ξi+1ˆBi+1,p−1(ξ),ˆBi,0(ξ):={1%ifξi≤ξ≤ξi+10otherwise, (18)

and are -times continuously differentiable across the -th knot with multiplicity .

The multivariate B-splines on the parameter domain , , are defined as a tensor-product of the corresponding univariate ones. In the multidimensional case, we define the knot-vector dependent on the coordinate direction , , where indicates the direction (in space or in time). Furthermore, we introduce a set of multi-indices , and a multi-index indicating the order of polynomials. The tensor-product of univariate B-spline basis functions generates multivariate B-spline basis functions

 ˆBi,p(ξ):=d+1∏α=1ˆBiα,pα(ξα),ξ=(ξ1,...,ξd+1)∈ˆQ.

The univariate and multivariate NURBS basis functions are defined in the parametric domain by means of B-spine basis functions, i.e., for the given and any , the NURBS basis functions are defined as

 ˆRi,p(ξ):=wiˆBi,p(ξ)W(ξ). (19)

Here, is the weighting function where .

The physical space-time domain is defined by the geometrical mapping of the parametric domain :

 Q:=Φ(ˆQ)⊂\mathdsRd+1,Φ(ξ):=∑i∈IˆRi,p(ξ)Pi, (20)

where are the control points. For simplicity, we assume the same polynomial degree for all coordinate directions, i.e., for all . By means of geometrical mapping (20), the mesh discretising is defined as . The global mesh size is denoted by

 h:=maxK∈Kh{hK},hK:=∥∇Φ∥L∞(K)^hˆK. (21)

Moreover, we assume that is a quasi-uniform mesh, i.e., there exists a positive constant independent of , such that

The finite dimensional spaces on are constructed by a push-forward of the NURBS basis functions, i.e.,

 Vh:=span{ϕh,i:=ˆRi,p∘Φ−1}i∈I, (22)

where the geometrical mapping is invertible in , with smooth inverse on each element , see, e.g., see [39, 3]. The subspace

 V0h:=Vh∩VΔx,10,0–

is introduced for the functions satisfying homogeneous boundary and initial conditions.

A stable space-time IgA scheme for (1) was presented and analysed in [30], where the authors proved its efficiency for fixed and moving spatial computational domains. In our analysis, we use spline bases of sufficiently high order, so that . In order to provide an efficient discretisation method, we consider (10), where and in (9) with some positive parameter and the global mesh-size defined in (21). It implies the stabilised space-time IgA scheme: find satisfying

 (∂tuh,vh+δh∂tvh)Q+(∇xuh,∇x(vh+δh∂tvh))Q=:as,h(uh,vh)=ls,h(vh):=(f,vh+δh∂tvh)Q. (23)

for all . The -coercivity of w.r.t. the norm

 |||vh|||2s,h:=∥∇xvh∥2Q+δh∥∂tvh∥2Q+∥vh∥2ΣT+δh∥∇xvh∥2ΣT (24)

follows from [30, Lemma 1]. As was noted in [30], coercivity implies that the IgA solution of (23) is unique. Moreover, since the IgA scheme (23) is posed in the finite dimensional space , uniqueness yields existence of the solution in . Moreover, following [30, 29], we can show boundedness of the bilinear form in (23) with respect to appropriately chosen norms. Combining coercivity and boundedness properties of with the consistency of the scheme and approximation results for IgA spaces implies a corresponding a priori error estimate presented in Theorem 3 below.

###### Theorem 3

Let , , , be the exact solution to (3), and let be the solution to (23) with some fixed parameter . Then, the following a priori error estimate

 ∥u−uh∥s,h≤Chr−1∥u∥Hr(Q) (25)

holds, where , and is a generic constant independent of .

Proof: See, e.g., [30, Theorem 8].

Corollary 1 presents a posteriori error majorants for and , where , .

###### Corollary 1

(i) If and , Theorem 1 yields the estimate

 (2−1γ)(∥∇xe∥2Q+δh∥∂te∥2Q) +∥e∥2ΣT+δh∥∇xe∥2ΣT=:|||e|||I,2s,h ≤¯¯¯¯¯MIs,h(v,y;γ,β,α):=γ(¯¯¯¯¯MI(v,y;γ,β)+δh˜MI(v,y;α)), (26)

where and are defined in Theorem 1 and the best and are given by relations

 β=CF∥rIeq∥Q∥