Computation of quadrature rules for integration with respect to refinable functions on assigned nodes

Applied Numerical Mathematics 90 (2015) 168–189

Contents lists available at ScienceDirect

Applied Numerical Mathematics www.elsevier.com/locate/apnum

Computation of quadrature rules for integration with respect to refinable functions on assigned nodes Francesco Calabrò a,∗ , Carla Manni b,1 , Francesca Pitolli c,2 a b c

DIEI, Università degli Studi di Cassino e del Lazio Meridionale, Italy Dipartimento di Matematica, Università di Roma “Tor Vergata”, Italy Dipartimento SBAI, Università di Roma “La Sapienza”, Italy

a r t i c l e

i n f o

a b s t r a c t

Article history: Received 16 April 2014 Received in revised form 15 September 2014 Accepted 19 November 2014 Available online 18 December 2014

Integrals involving refinable functions are of interest in several applications ranging from discretization of PDEs to wavelet analysis. We present a procedure to construct quadrature rules with assigned nodes for these integrals. The process requires in input the refinement mask coefficients and the sequence of nodes only. The corresponding weights are computed by an iterative procedure that does not involve the solution of linear systems. The proposed approach is deeply based on the strong connection between balanced measures and integrals of refinable functions. © 2014 IMACS. Published by Elsevier B.V. All rights reserved.

Keywords: Refinable functions Quadrature rules Isogeometric analysis Wavelets

1. Introduction A refinable function φ is a solution of the refinement equation

φ(x) =



a j φ(2x − j ),

(1.1)

j ∈Z

where the set of real numbers a := {a j , j ∈ Z} is the so-called (refinement) mask, see [6,9] for details. Very often in applications the mask a is compactly supported, i.e. it has only finitely many elements different from zero. If a is compactly supported, without any loss of generality, for some m ∈ N we can assume

a j = 0,

j∈ / [0, m],

a0 am = 0,

(1.2)

so that (1.1) reads

φ(x) =

m 

a j φ(2x − j ).

j =0

* 1 2

Corresponding author. Tel.: +39 0776 2993621. E-mail addresses: [email protected] (F. Calabrò), [email protected] (C. Manni), [email protected] (F. Pitolli). Tel.: +39 06 72594660; fax: +39 06 72594699. Tel.: +39 06 49766631; fax: +39 06 4957647.

http://dx.doi.org/10.1016/j.apnum.2014.11.010 0168-9274/© 2014 IMACS. Published by Elsevier B.V. All rights reserved.

(1.3)

F. Calabrò et al. / Applied Numerical Mathematics 90 (2015) 168–189

169

The most popular (compactly supported) refinable functions are surely the cardinal B-splines, i.e., B-splines on integer knots, and (father) functions for Daubechies orthonormal wavelets, see Section 4 and references therein. Here we are interested in the computation of integrals of the form

+∞ f (x)φ(x) dx,

(1.4)

−∞

where φ is a given refinable function. From (1.2)–(1.3) it follows that if a is compactly supported as in (1.2) also φ is compactly supported on [0, m] and so the integral (1.4) becomes

m f (x)φ(x) dx.

(1.5)

0

Integrals of type (1.4) or (1.5) are required in several contexts as pointed out below. 1.1. Applications and related areas Among the numerous applications that require the evaluation of integrals involving refinable functions, we shortly outline the three following relevant cases.

• Isogeometric analysis (IgA) is a recently proposed computational technique for the solution of boundary value problems involving PDEs. The key ingredient of IgA is the use of the same functions adopted in Computer Aided Design (CAD) both to describe the geometry of the domain and to obtain the numerical approximation of the solution, [7]. Therefore B-splines, and their rational counterpart (NURBS), are the main building blocks in IgA. When a Galerkin approach is considered in NURBS-based IgA, it is imperative to compute the discrete problem of the approximation of the PDEs via quadrature rules able to properly integrate the underlining B-spline space functions. In order to fix the ideas, consider the calculation of the mass matrix B = (blk ), the stiffness matrix S = (slk ) and the linear advection matrix T = (tlk ). As discussed in [1,7] the calculations are made in the parametric domain considering the following univariate integrals:

blk :=

+∞ β(x)ϕl (x)ϕk (x)dx,

(1.6)

−∞

+∞

σ (x)ϕl (x)ϕk (x)dx,

(1.7)

+∞ ϑ(x)ϕl (x)ϕk (x)dx, tlk :=

(1.8)

slk := −∞

−∞

where {ϕl }l∈ I is a basis for the considered discretization space and β, σ , ϑ are known functions taking into account the coefficients of the investigated partial differential equation and the problem geometry. In the case of3 a periodic knot vector4 with uniform spacing all the ϕl can be chosen as the translated copies of the cardinal B-spline of degree n, which is supported in [0, n + 1], see Section 4. Quadrature rules for the evaluation of (1.4) are well suited to approximate the entries of the mass and of the advection matrix. Indeed, a quadrature formula exact for B-splines can be used in order to integrate the product of the first two functions instead of applying a standard quadrature on the product of the three functions. Recalling that the derivative of a B-spline can be expressed in terms of linear combinations of B-splines of lower degree, see [3], also the evaluation of the elements of the stiffness matrices can be written as integrals of the form (1.4). As for all Galerkin methods, the numerical evaluation of the expressions in (1.6)–(1.8) is a bottle neck of the process and corresponding efficient quadrature rules are crucial for the efficiency of the whole IgA Galerkin approximation process. This motivated the large interest recently received by this topic [33]. Optimal quadrature rules for the mass and stiffness matrices in uniform B-spline IgA discretizations are constructed in [23]. The computation of the nodes and weights results in a global, non-linear system over the whole domain of

3 4

For all definitions adopted in this section we refer to [1]. More attention has to be paid when an open knot vector is considered, because of boundary functions.

170

F. Calabrò et al. / Applied Numerical Mathematics 90 (2015) 168–189

the B-spline space. The spline space of the product of two uniform B-spline basis functions is further addressed in [1], leading to a rule which can be obtained as the solution of a local non-linear system. A variety of quadrature rules for isogeometric discretizations of low degree based on local Bézier representations have been investigated in [36]. However, up to our knowledge, so far the refinability property of the cardinal B-spline has not been explored in order to optimize the calculations in (1.6)–(1.8), see also [5]. • Many applications, such as adaptive numerical schemes for solving PDEs or wavelet analysis, require the evaluation of integrals of the type

+∞ f (x)ψ(x − k) dx,

k∈Z

(1.9)

−∞

where ψ is a wavelet. Well known examples are compactly supported wavelets associated with orthogonal Daubechies refinable functions and biorthogonal B-spline refinable bases, see, for instance, [9]. Usually, wavelets are constructed starting from a given refinable function φ . Indeed a wavelet can be obtained as a linear combination of the integer translates of φ(2·), i.e.:

ψ(x) =



b j φ(2x − j ),

(1.10)

j ∈Z

where the sequence b = {b j , j ∈ Z} is the wavelet mask. We notice that if φ is compactly supported, ψ can be compactly supported, too. In the salient cases only finitely many elements of the sequence b are non-zero. Inserting the wavelet equation (1.10) in the integral (1.9), the evaluation of (1.9) can be decomposed in subsequent evaluations of integrals of the form (1.4) or (1.5). • Due to the refinability property of φ , integrals of the type (1.4) can be also seen as refinable functionals. Indeed, (1.4) is a linear functional on a suitable space of functions. From (1.3) we have

+∞ f (x)φ(x)dx =

m 

+∞ f (x)φ(2x − j )dx =

aj

j =0

−∞

 +∞  m  aj y+ j j =0

−∞

2

f

−∞

2

φ( y )d y .

(1.11)

Therefore, (1.4) satisfies

+∞ f (x)φ(x)dx =

 +∞  m  aj x+ j j =0

−∞

2

f

−∞

2

φ(x)dx,

(1.12)

which is usually referred to as balancing equation, see [24], and can be also taken as a definition for refinable functionals, see [31]. 1.2. Problem setting Due to the relevance of the problem in the applications in this paper we focus on the evaluation of (1.5) where φ is a refinable function defined by (1.1) and a is a compactly supported mask as in (1.2). Furthermore, we assume, as it is usually the case, that φ is normalized so that

m M0 :=

φ(x)dx = 1.

(1.13)

0

Since M0 = 0, from (1.1) it follows m 

a j = 2.

(1.14)

j =0

Existence and uniqueness of solutions of (1.1) with the normalization (1.13) are guaranteed under suitable hypotheses on the mask a, see [9]. In practice, integrals of type (1.5) are approximated by using a suitable quadrature formula, defined by pairs of nodes and weights {(ξi , ωi ), i = 0, . . . , N }, such that N  i =0

m

ωi f (ξi ) ≈

f (x) φ(x) dx. 0

(1.15)

F. Calabrò et al. / Applied Numerical Mathematics 90 (2015) 168–189

171

In this paper we are interested in constructing quadrature rules for approximating (1.5) when the nodes

ξi ,

i = 0, . . . , N

are given in input. This is often the case in practical applications whenever f is known or can be (efficiently) evaluated at some specific sites only. As usual, the weights in (1.15) will be determined by means of the so-called exactness equations. Therefore, we aim to solve the following problem. Problem 1.1. Assume that φ solves the refinement equation (1.3) for a suitable mask coefficients as in (1.2). Let the N + 1 points

{ξi , i = 0, . . . , N },

0 ≤ ξ i < ξ i +1 ≤ m ,

(1.16)

be given, find the weights {ωi , i = 0, . . . , N } such that N 

m

ωi (ξi ) p =

i =0

x p φ(x)dx,

p = 0, . . . , N .

(1.17)

0

It is well know that the above problem has a unique solution for any selection of distinct nodes. 1.3. Related literature Quadrature formulae for integrals involving refinable functions have been addressed by several authors. Most of them focused either on (almost) equally spaced nodes or on Gaussian nodes. Quadrature rules over equally spaced nodes are obtained in [29] by means of the inversion of a Vandermonde like matrix. In [8] the integral (1.5) is approximated by a quadrature rule on equally spaced nodes with shift, while the weights are given by the mask coefficients. Equally spaced nodes with shift are considered in [37] too, where the shift is selected to obtain degree of exactness N and the construction is carried out by explicitly computing the so-called modified moments (see Section 2). Such an approach is investigated in the context of composite rules in [26]. Composite rectangular rules ensuring degree of exactness n for cardinal B-splines of degree n are addressed in [39]. Positive refinable functions have been considered as weights of Gaussian quadrature rules by several authors. The case of B-splines has been addressed in [34], where an approximate procedure is used to compute associated orthogonal polynomials via the classical three-term recurrence relation. In [17] the approach is extended to a larger class of positive refinable functions, the so-called GP refinable ripplets [20]. The coefficients of the three-term recurrence relation for the corresponding orthogonal polynomials are computed by an approximate procedure based on the mask coefficients only. Iterative computation for coefficients of the three-term recurrence relation is performed in [2] for B-splines and in [31] for general masks by using a procedure previously presented in [32], see also [30]. It is well known that, if the nodes (1.16) are given, the highest degree of exactness, which is at least N, is provided by the (unique) interpolatory quadrature rule. From the theoretical point of view, the corresponding weights could be computed by solving a linear system. The choice of the basis in the space where exactness is required determinates the linear system to be solved both for the matrix and the right hand side. The straightforward construction uses the monomial bases. In this case the right hand side reduces to ordinary moments of φ , while the involved matrix is a Vandermonde like one, which is ill-conditioned, [15,25]. Other choices of the basis require the computation of the so-called modified moments of φ . We will refer to these methods as exactness-based methods. Another method for the calculation of the quadrature weights in the case of preassigned nodes is the one proposed by Kautsky in [28] and fully described in [13]. This method calculates the interpolating quadrature rule once the three term relation defining the polynomials that are orthogonal with respect to the weight function involved is given. The method is purely algebraic and uses the relation between the weights of a quadrature rule and the principal vectors of some matrices, as for the case of the calculation of Gaussian quadrature via eigenvalues and eigenvectors of the Jacobi matrix itself, see [16]. To the best of our knowledge, the problem of constructing quadrature rules for refinable functions with general given nodes has been previously addressed in [27] only. Nevertheless, the computation of the weights in [27] is reduced to the evaluation of derivatives of Lagrange polynomials, so it requires again the solution of a Vandermonde like linear system. 1.4. Proposed approach and comparisons The aim of this work is threefold: (i) We introduce a procedure for the calculation of the modified moments. This procedure can be profitably used in order to apply the exactness-based methods reviewed in the previous subsection to solve the Problem 1.1. Moreover, we notice that Kautsky’s method can be applied merging the procedure proposed in [13] with the calculation made in [31,32]. Up to our knowledge these methods have not been yet adapted for the calculation of the integrals in (1.5) with arbitrary but a priori assigned nodes.

172

F. Calabrò et al. / Applied Numerical Mathematics 90 (2015) 168–189

(ii) We propose a new iterative procedure to compute the required weights via modified moments. Our method closely follows the classical construction of Gaussian quadrature rules by means of orthogonal polynomials and three-term recurrence relations: the main difference consists in the measure we deal with, actually a discrete measure concentrated on the nodes (1.16). More precisely, our procedure consists in the following steps: – to determine the family of orthogonal polynomials with respect to a discrete measure defined by the given nodes; – to compute the modified moments of φ with respect to the above polynomials; – to compute the required weights in term of the above modified moments and orthogonal polynomials. Since in our case the nodes are given as input, our approach can be seen as a generalization of the method presented in [30] and [31]. On this concern we remark that our procedure provides an efficient tool to construct sequences of interpolatory quadrature rules on prescribed nodes such as Chebychev nodes, see [38]. (iii) Finally we test the methods in several cases. First we present the range of applicability of these numerical procedures. Since we are considering a quadrature on assigned nodes, the existence of the rule is guaranteed no matter if the φ is sign changing, while for the methods that deal with the construction of weighted type Gauss rules this feature has to be carefully considered. Then we perform convergence tests. The new procedure can be compared with the method by Kautsky: both of them deeply exploit the three term recurrence relation for suitable families of orthogonal polynomials. As testified by several numerical experiments the two approaches behave in a completely similar way. Nevertheless, the new procedure is especially tailored for integrals involving refinable functions and it is extremely simple to be described and coded. As a consequence it offers a different perspective and constitutes a valid alternative with respect to [13,28]. Since the goal of this paper is to present a procedure to construct the weights of a quadrature rule for integrals involving refinable functions assuming the nodes are given in input, we do not address here the problem of the selection of the nodes and we test the procedures for standard choices. The remaining of the paper is divided into five sections. Section 2 describes the details of the recursive computation of modified moments. The new proposed method for the computation of the quadrature weights is addressed in Section 3. In order to focus on the wide range of applicability, we collect in Section 4 the masks of some of the most common refinable functions summarizing their distinguish properties. Section 5 presents some calculated quadratures and numerical tests that put in evidence the behavior of the proposed construction also compared with the other techniques. Finally, in Section 6 we present some concluding remarks and we outline some possible generalizations to more general masks. 2. Modified moments There exist different approaches to solve Problem 1.1. As stated in Section 1.3, the so-called exactness-based methods directly address the linear system defined by (1.17) to compute the weights. Besides monomials, suitable bases can be considered to formulate the exactness conditions in order to improve the conditioning of the resulting matrix. Among the others orthogonal polynomials are of great interest. Let { P j , j = 0, . . .} be a system of algebraic polynomials, with P j of degree equal to j, mutually orthogonal on the interval [0, m] with respect to a measure that we will denote by μ. With such a choice (1.17) are equivalent to N 

ωi P j (ξi ) = L j ,

j = 0, . . . , N ,

(2.1)

i =0

where

m Li :=

P i (x)φ(x) dx,

i = 0, . . .

(2.2)

0

denote the modified moments of φ see, for example, [16, Section 2.1]. The aim of this section is to compute the modified moments Li by a recursive procedure. A procedure for these calculations is available in [30]. Because these modified moments are a key ingredient in our procedure we present in this section complete proofs in slightly different hypotheses. We recall that the moments

m xi φ(x)dx,

Mi :=

i = 1, . . .

0

(see the exactness equations (1.17)) can be explicitly calculated via a recursive relation. This result is well known (see ad ex. [8, Section 4]) and we report it for completeness.

F. Calabrò et al. / Applied Numerical Mathematics 90 (2015) 168–189

173

Proposition 2.1. Let M0 be known. Then the moments Mi can be calculated by the following recursive relation: m 

1

Mi =

2(2i − 1)

aj

i 

  jq

q =1

j =1

i

q

Mi −q ,

i = 1, . . . .

Here we want to follow a similar strategy to compute the general modified moments (2.2) by an iterative procedure. The procedure we are describing can be derived from the one described in [30]. Similar ideas were also used in [31,32,37]. We recall that the family { P i , i = 0, . . .} of orthogonal polynomials is connected by the following three-term relation5 [16]:

P −1 (x) = 0,

P 0 (x) = 1,

P i (x) = (x − αi −1 ) P i −1 (x) − βi −1 P i −2 (x),

(2.3)

αi and βi are known for a large class of measures. In general, they can be calculated by the following

The coefficients expressions6 :

αi =

i = 1, . . . .

xP i , P i μ , P i , P i μ

i = 0, . . . ;

βi =

P i , P i μ , P i −1 , P i −1 μ

i = 1, . . . ;

(2.4)

where the symbol ·, · μ denotes the scalar product with respect to the measure μ. For the following computations we need the next lemma, involving shifted orthogonal polynomials. Lemma 2.2. Let { P i , i = 0, . . . , N } be a family of orthogonal polynomials as in (2.3). Then,

 Pi

x+ j

 =

2

i 1  ( j) cl,i P l (x), 2i

(2.5)

l =0

( j)

where the coefficients cl,i can be computed in a recursive manner. More precisely, if polynomials (2.3)–(2.4), then:

αi , βi are the recurrence coefficients for the

( j)

c i , i = 1, ( j)

( j)

( j)

( j)

( j)

cl,i = cl−1,i −1 + ( j − 2αi −1 + αl )cl,i −1 + βl+1 cl+1,i −1 − 4βi −1 cl,i −2 ,

l = 0, . . . , i − 1 ,

(2.6)

where coefficients with at least one negative index are set to be zero. ( j)

Proof. First, we notice that c i ,i = 1 because the polynomials are monic. ( j) Then, using the definition of the coefficients cl,i we can write:

 P i −1

x+ j

 =

2

1 2 i −1

i −1  l =0

( j)

cl,i −1 P l (x),

 P i −2

x+ j



2

=

1 2 i −2

i −2  l =0

( j)

cl,i −2 P l (x).

(2.7)

Now, we recall that the polynomials satisfy the following three term relation, see (2.3):

P −1 (x) = 0,

P 0 (x) = 1,

P i (x) = (x − αi −1 ) P i −1 (x) − βi −1 P i −2 (x),

i = 1, . . . .

(2.8)

From (2.8) we have:

 Pi

x+ j 2



 =

x+ j 2



     x+ j x+ j − α i −1 P i −1 − β i −1 P i −2 . 2

2

Substituting (2.7) into (2.9) we obtain:

5 6

As well known, the relation is defined up to a multiplicative constant. Here, we will consider the monic case. For the sake of simplicity x denotes either the independent variable or the identity function.

(2.9)

174

F. Calabrò et al. / Applied Numerical Mathematics 90 (2015) 168–189

 Pi



x+ j

 =

2

=

x+ j

 − α i −1

2 i −1 1 

2i

l =0



i −1 

1 2 i −1



( j)

cl,i −1 xP l (x) +

j 2

l =0

1

( j)

cl,i −1 P l (x) − βi −1



− α i −1

i −1 1 

2 i −1

l =0

2 i −2

i −2  l =0

( j)

cl,i −2 P l (x) 1

( j)

cl,i −1 P l (x) − βi −1

2 i −2

i −2  l =0

( j)

cl,i −2 P l (x).

Now, we use (2.8) to express xP l (x) = P l+1 (x) + αl P l (x) + βl P l−1 (x). Substituting and changing the summation index we finally obtain:

 Pi



x+ j 2



i i −1  1  ( j) ( j) = i cl−1,i −1 P l (x) + cl,i −1 αl P l (x) 2 l =1

i −2 

+

l =0

l =0

( j) cl+1,i −1 βl+1 P l (x) + ( j

− 2α i −1 )

i −1  l =0

( j) cl,i −1 P l (x) − 4βi −1

i −2  l =0

 ( j) cl,i −2 P l (x)

.

( j)

2

By (2.5) the left side term can be written in terms of cl,i , thus the previous equation gives the recursive relation (2.6).

With the previous lemma, we can use the balancing property (1.12), replacing f by the polynomial P i in order to obtain a recursive relation for the modified moments. The following proposition extends to a general mask the results in [37, Section 2.5.3] using the procedures of [32,31], compare also with what described in [30, Section 3]. ( j)

Proposition 2.3. Let cl,i be as in Lemma 2.2 and set L0 := M0 . Then, the modified moments Li can be calculated recursively by:

Li =

m i −1  

1 2(2i − 1)

 ( j) a j c l ,i

i = 1, . . . .

Ll ,

(2.10)

j =0

l =0

Proof. We start from definition (2.2) and the balancing property (1.12) to write:

+∞ Li =

P i (x)φ(x)dx =

 +∞  m  aj x+ j j =0

−∞

2

Pi

−∞

2

φ(x)dx.

(2.11)

Therefore, from (2.5) and (1.14), we get:

Li =

+∞  i m  aj 1 j =0

=

2

−∞

2i

i m 1  aj

2i

l =0 j =0

2

Hence we obtain (2.10).

l =0

( j)

1 ( j) cl,i P l (x)φ(x)dx = i 2

cl,i Ll =

1 2i

 Li +

m i  aj  j =0

m i −1   aj l =0 j =0

2

2

( j) c l ,i

l =0



( j)

cl,i Ll =

m P l (x)φ(x)dx 0

1

1

2

2 i +1

Li + i

i −1  m 

( j)

a j cl,i Ll .

l =0 j =0

2

3. The new procedure In this section we propose an alternative procedure to solve Problem 1.1. To this aim we consider a well established approach that uses exact integration of an appropriate approximation of the integrand function f . More precisely, we fix a discrete measure concentrated on the nodes (1.16), i.e.

γ :=

N 

c i δx−ξi ,

(3.1)

i =0

where c i are so far positive real numbers and δx−ξ denotes the Kronecker delta symbol: δx−ξ = 1 if x = ξ , δx−ξ = 0 if x = ξ . Then we have the following: Proposition 3.1. Let { Q j , j = 0, . . .}, be a system of algebraic polynomials, mutually orthogonal on the interval [0, m] with respect to the discrete measure γ as in (3.1). Then the weights {ωi , i = 0, . . . , N } that solve Problem 1.1 can be written as:

F. Calabrò et al. / Applied Numerical Mathematics 90 (2015) 168–189

ωi =

N  L j c i Q j (ξi ) , N 2 l=0 cl Q j (ξl ) j =0

i = 0, . . . , N ,

175

(3.2)

where L j are the modified moments (2.2) constructed with respect to the discrete measure γ , i.e. taking Q j instead of P j . Proof. Write the following approximation for the function f , [35, Section 10.1]:

f (x) ≈

N  f , Q j γ Q j (x), Q j , Q j γ j =0

where

m f , g γ :=

f (x) g (x)dγ . 0

Therefore, we consider the following approximation of the integral (1.5):

m f (x)φ(x)dx ≈

m N  f , Q j γ Q j (x)φ(x) dx. Q j , Q j γ j =0

0

0

By construction, such an approximation provides the exact value of less than or equal to N. N m From (3.1) we have 0 g (x)dγ = i =0 c i g (ξi ), therefore

m 0

f (x)φ(x) dx whenever f is a polynomial of degree

 N  m N N N     L j c i Q j (ξi ) f , Q j γ f , Q j γ Q j (x)φ(x)dx = Lj = f (ξi ), N 2 Q j , Q j γ Q j , Q j γ l=0 cl Q (ξl ) j =0

0

so that, setting

j =0

i =0

j =0

j

ωi as in (3.2), we obtain the coefficients solving Problem 1.1. 2

We note that expression (3.2) above can be simplified by choosing in (3.1) all the coefficients c i equal to a given value. With this choice we get

ωi =

N  L j Q j (ξi ) , N 2 l=0 Q j (ξl ) j =0

i = 0, . . . , N .

(3.3)

Finally, in order to explicitly construct the required quadrature formula we need an efficient procedure to construct the set of polynomials

{ Q j , j = 0, . . . , N }

(3.4)

of degree j that are orthogonal with respect to the discrete measure7 γ given in (3.1). The family (3.4), follows a three-term recurrence relation as in (2.3)–(2.4). If the sampler measure μ is discrete, the coefficients can be calculated directly by (2.4). This is just our case, since the measure μ we are dealing with is the measure γ defined in (3.1). We notice that the problem can be solved in an exact manner if exact computation is available. By the other side, it has been noticed that the three-term relation can be unstable, thus there is a wide literature on how to construct the coefficients αi , βi , see [16, Section 2.2]. For details on our choices when implementing the procedure, see Section 5. Summarizing, the new procedure for the construction of the weights of the quadrature rule (1.15) is as follows: 1. Get in input the points {ξi , i = 0, . . . , N }; 2. Fix a discrete measure on the nodes as in (3.1); 3. Calculate the recurrence coefficients for the (discrete) orthogonal polynomials { Q j , j = 0, . . . , N } via8 (2.4); ( j)

4. Calculate the recursive coefficients cl,i according to (2.6); 5. Calculate the modified moments {L j , j = 0, . . . , N } according to (2.10); 6. Calculate the weights according to (3.2).

7 8

Actually, we can also construct a polynomial of degree N + 1 with the same properties, but it vanishes at all the ξi , i = 0, . . . , N. Recall that (2.4) can be explicitly calculated because γ is discrete.

176

F. Calabrò et al. / Applied Numerical Mathematics 90 (2015) 168–189

Extensive numerical experiments to show the properties of such an approach will be presented in Section 5. Remark 3.2. As stated in the Introduction, the weights in Problem 1.1 can be computed by solving a linear system where each equation corresponds to exactness of the quadrature rule on an element of a basis of the space of polynomials of degree less than or equal to N. Besides monomials, suitable bases can be considered to improve the conditioning of the resulting matrix. The procedure we have presented can be seen as a useful procedure to solve the above mentioned linear system in the case when the family of orthogonal polynomials is determined by the discrete measure (3.1). Moreover we notice that the choice of the coefficients c i in the discrete measure (3.1) will not affect the definition of the quadrature rule. Indeed, the quadrature rule in (1.15) is uniquely determined by the nodes (1.16) and the exactness equations (1.17). 4. Some compactly supported masks For the sake of completeness in this subsection we collect some well known (compactly supported) masks corresponding to refinable functions with different features. The properties of the mask a and the associated refinable function φ , determine the behavior of the space of integer translates of φ , that is the properties of functions of the form:



qk φ(x − k),

qk ∈ R.

(4.1)

k

Among the various properties we recall that a mask is called variation diminishing if the number of sign changes of (4.1) is bounded above by the number of sign changes in the sequence of the coefficients {qk }. The mask is interpolatory if for some integer r we have

d2 j = δ0− j ,

where d j := a j +r , j ∈ Z.

This implies that the function in (4.1) interpolates the sequence (k + r , qk ), k ∈ Z. Both the above properties are of salient interest in approximation theory and computer aided geometric design. 4.1. Cardinal B-splines The most popular refinable functions are the cardinal B-splines, i.e. minimally supported piecewise polynomials on integer knots, see [3] for details. The cardinal B-spline of degree n is compactly supported on [0, n + 1], strictly positive on (0, n + 1), centrally symmetric and belongs to C n−1 . The corresponding refinement mask, which is not interpolatory but variation diminishing, is [6]:

 a = a(n) := a j (n), j = 0, . . . , n + 1 ,

a j (n) =



n+1

1 2n



.

j

(4.2)

We will denote the corresponding refinable functions by Nn . In particular, in the following section we will construct quadrature rules for the cardinal B-splines of degree n = 2 and n = 3, i.e. N2 and N3 . Their behavior is displayed in Fig. 1 and the corresponding masks are



a(2) =

1 3 3 1

, , ,

4 4 4 4



,

a(3) =

1 1 3 1 1

, , , ,

8 2 4 2 8

(4.3)

.

4.2. Orthogonal Daubechies refinable functions Daubechies refinable functions are compactly supported orthogonal functions corresponding to masks satisfying the orthogonality condition

     p a (ξ )2 +  p a (ξ + π )2 = 2,

(4.4)

where the trigonometric polynomial

p a (ξ ) =

m 

a j e −i j ξ ,

j =0

is the so-called mask symbol, see [9, Chapter 6] for details. The coefficients of the Daubechies masks are obtained by solving the polynomial equation (4.4) with some suitable additional assumptions. The mask coefficients9 determining refinable functions associated to compactly supported wavelets

9

Note that the mask coefficients in [9] are normalized so that they sum up to

√

2.

F. Calabrò et al. / Applied Numerical Mathematics 90 (2015) 168–189

177

Fig. 1. The cardinal B-splines of degree n = 2 (left) and n = 3 (right). The corresponding masks are reported in Eq. (4.3).

Fig. 2. The Daubechies refinable functions for p = 2 (right) and p = 6 (left). For the first the corresponding mask is reported in Eq. (4.5).

with maximum number of vanishing moments, p, for their support width and extremal phase choice can be found in [9, Table (6.1)]. In this case m = 2p − 1. The coefficients do not have an explicit expression but for p = 1 and p = 2, and parameterizing the mask by p, that is a = a( p ), we have

a(1) = {1, 1},

a(2) =

1+

√

4

3 3+

,

√

4

3 3−

,

√

4

3 1−

,

√

4

3

.

(4.5)

When p = 1, the corresponding refinable function is the characteristic function of the interval [0, 1), which is the unique positive and symmetric orthogonal refinable function. For p > 1 the refinable functions are non-symmetric oscillating functions and do not have an explicit expression. The Daubechies refinable functions are displayed in Fig. 2 for p = 2 and p = 6. As the figure shows, their smoothness increases as the length of the refinement mask increases. Approximatively, they are of class C 0.19( p −1) for p sufficiently large, see [9, Chapter 7]. 4.3. Refinable ripplets A ripplet is a function f such that its integer translates { f (· − k), k ∈ Z} provide a totally positive system, i.e., all the minors of any collocation matrix are non-negative. Therefore, ripplets are non negative and variation diminishing, see [18]. Refinable ripplets, i.e., ripplets satisfying a refinement equation, were considered in [19], where conditions on the refinement mask, in order the corresponding refinable function is a ripplet are given. A well known example of refinable ripplets are the cardinal B-splines, see Section 4.1. A further example is given by the class of GP refinable functions introduced in [20] and associated with the masks

178

F. Calabrò et al. / Applied Numerical Mathematics 90 (2015) 168–189

Fig. 3. The GP refinable ripplet for n = 3, h = 5. The corresponding mask is reported in (4.6).





a = a(n, h) := a j (n, h), j = 0, . . . , n + 1 , a j (n, h) :=

n

1



n+1

2h

h ≥ n , n = 2, . . . ,    n − 1 + 4 2h−n − 1 , j−1





j

where we set l = 0 if 0 < l or l > n. The corresponding refinable ripplets are compactly supported on [0, n + 1], centrally symmetric, strictly positive on (0, n + 1) and belong to C n−2 . The real number h is a shape parameter. For h = n the GP refinable function reduces to the cardinal B-spline of degree n and class C n−1 . When h increases the GP refinable function resembles more and more to the cardinal B-spline of degree n − 2. The GP refinable ripplet corresponding to the mask



a(3, 5) =

1

1

1

1

2

2

2

25

,1 − 5

,1 − 5

, 5

(4.6)

is plotted in Fig. 3. 4.4. Interpolatory 4-point scheme In [11] the one parameter family of symmetric masks

a = a( w ) := − w , 0,

1 2

+ w , 1,

1 2

+ w , 0, − w ,

w ∈ R,

(4.7)

was introduced, see also [12] where the properties of the mask (4.7) were analyzed in the context of geometric modeling. The mask is interpolatory and the corresponding refinable functions are symmetric functions compactly supported on [0, 6]. For w = 0√we recover the mask of the cardinal B-spline of degree 1, i.e., the hat function, properly translated in [2, 4]. For all 0 < w < 58−1 ≈ 0.1545 the corresponding refinable functions belong to C 1 . In Fig. 4 the refinable functions corresponding to the masks

 

a

1 5



1 7 7 1 = − , 0, , 1 , , 0, − , 5

10

10

5



a

1 16





1 9 9 1 = − , 0, , 1 , , 0, − 16

16

16

16

(4.8)

are displayed. 5. Numerical tests In this section we present some numerical results on quadrature formulae of the type (1.15). Our aim is to present the resulting quadratures on assigned nodes and to test this quadratures when the number of nodes increases. The performed numerical tests give a numerical evidence that the quadrature rules based on Chebychev type nodes work as well as the Gaussian quadrature rules. For the calculation of the required quadratures we will consider three different methods: CMP: Our procedure as described in Section 3. K: The method proposed by Kautsky, see [13,28].

F. Calabrò et al. / Applied Numerical Mathematics 90 (2015) 168–189

179

Fig. 4. The interpolatory refinable functions corresponding to the 4-point mask with w = 1/5 (left) and w = 1/16 (right). The corresponding masks are reported in (4.8).

V: The method obtained by solving the exactness equations (1.17) that requires the inversion of a Vandermonde matrix. Method K can be applied combining the techniques in [31,32] for the calculation of the recursive coefficients of the orthogonal polynomials with the procedure in [13] that allows to compute the interpolating quadrature rule. Method CMP needs the computation of the discrete-orthogonal polynomials that has not been yet discussed. For the calculations made in the present work, we have used the routine available by W. Gautschi at the URL http://www.cs.purdue.edu/archives/2002/wxg/ codes/lanczos.m. This routine, described in [16, Section 2.2.3], is based on an algorithm by W.B. Gragg and W.J. Harrod [22]. Other procedures are available for the same purpose, and we have performed some tests in reference cases in order to fix our choice comparing, e.g., the technique described in [14] or the simple Stjties algorithm. In all of our tests the best results are always achieved by Gautschi’s routine. All presented tests are made in Matlab environment using a laptop PC. Remark 5.1. The numerical tests we performed show some instabilities due to the fixed (double) precision used by Matlab. We point out that in order to test the construction of quadrature formulae in terms of stability with respect to round-off errors this choice should be discarded, and some symbolic calculation performed (as done ad ex. in [31]). Stability tests of this type are beyond the scope of the present paper. On the other hand, the aim of our tests is to show how to use these formulae practically in a software environment widely employed in mathematical simulations. 5.1. Computed quadrature rules In this subsection, we present some results on calculated quadrature formulae with respect to the refinable functions seen in Section 4. We have considered two families of quadrature nodes: OEq: Open-equispaced nodes



m ( j + 1) N +2

, j = 0, . . . , N .

(5.1)

Interpolating quadrature on equispaced nodes suffers of the (usual) instability10 due to the oscillations of the interpolating polynomial [25]; for this reason we will present results only in the case of low degrees. Cheb: Chebychev type11 nodes

m

1 + cos(π ( j + 1)/( N + 2)) 2 

, j = 0, . . . , N .

(5.2)

The resulting quadrature rule is such that i |ωi | is rapidly increasing. Notice that it is not the usually used Chebychev rule, indeed it is called Fejer second rule, see [38] and references therein. This choice suits better to functions which vanish at the endpoints, as is the case for compactly supported continuous refinable functions. 10

11

180

F. Calabrò et al. / Applied Numerical Mathematics 90 (2015) 168–189

Table 1 Some calculated rules for B-splines of degree 2 and 3 and nodes as in (5.1)–(5.2). 7 nodes of open-equispaced type B-spline of degree 2

B-spline of degree 3

ω

Nodes ξ

Weights

0.375000000000000 0.750000000000000 1.12500000000000 1.50000000000000 1.87500000000000 2.25000000000000 2.62500000000000

0.0289909348894260 0.0943671478650902 0.250501883423694 0.252280067643579 0.250501883423695 0.0943671478650902 0.0289909348894260

ω

Nodes ξ

Weights

0.500000000000000 1 1.50000000000000 2 2.50000000000000 3 3.50000000000000

0.0100529100529101 0.0841269841269840 0.239682539682540 0.332275132275133 0.239682539682540 0.0841269841269840 0.0100529100529100

7 nodes of Chebychev type B-spline of degree 2

B-spline of degree 3

ω

Nodes ξ

Weights

0.114180701233070 0.439339828220179 0.925974851452365 1.50000000000000 2.07402514854764 2.56066017177982 2.88581929876693

0.00181968822010757 0.0385524854660656 0.237124723563811 0.445006205500033 0.237124723563811 0.0385524854660656 0.00181968822010759

ω

Nodes ξ

Weights

0.152240934977427 0.585786437626905 1.23463313526982 2 2.76536686473018 3.41421356237310 3.84775906502257

0.000715767513804665 0.0160714285714284 0.225474708676672 0.515476190476191 0.225474708676672 0.0160714285714286 0.000715767513804790

15 nodes of Chebychev type B-spline of degree 2

B-spline of degree 3

ω

Nodes ξ

Weights

0.0288220793951544 0.114180701233070 0.252795581546182 0.439339828220179 0.666644650470597 0.925974851452365 1.20736451697581 1.50000000000000 1.79263548302419 2.07402514854764 2.33335534952940 2.56066017177982 2.74720441845382 2.88581929876693 2.97117792060485

1.99964819253518e−05 0.000752751577043554 0.00517944610667290 0.0202136694688080 0.0541435579935992 0.117370137006536 0.191955949027195 0.220728984676440 0.191955949027195 0.117370137006536 0.0541435579935992 0.0202136694688080 0.00517944610667289 0.000752751577043515 1.99964819253805e−05

ω

Nodes ξ

Weights

0.0384294391935391 0.152240934977427 0.337060775394909 0.585786437626905 0.888859533960796 1.23463313526982 1.60981935596774 2 2.39018064403226 2.76536686473018 3.11114046603920 3.41421356237310 3.66293922460509 3.84775906502257 3.96157056080646

3.44899618193146e−07 8.97837123781060e−05 0.00139009286027726 0.00930007492507507 0.0382858335518381 0.110426089303495 0.209926903291441 0.261161754911755 0.209926903291441 0.110426089303495 0.0382858335518382 0.00930007492507497 0.00139009286027729 8.97837123780856e−05 3.44899618208148e−07

The distribution of the Chebychev nodes is similar to the Gauss ones: this property leads to an interpolatory quadrature that well compares with the Gauss one also in the case of singular measures [4,38]. We have chosen to use this family of nodes in order to test our routine also for high degrees. First of all, we consider the B-spline case because of their popularity and also because of their representation as piecewise polynomials, a property that we will use in the next section in order to construct test cases with known exact solutions. In Table 1 we show the resulting quadrature rules in the case of the B-splines of degree 2 and 3. In Figs. 5–6, the quadrature weights are plotted as a function of the given quadrature nodes. The plots on the left refer to 7 nodes chosen as in (5.1) and (5.2); the plots on the right refer to Chebychev nodes with N = 2k − 2, k = 3, . . . , 6. In this way the resulting nodes are nested and this favorable property can be very useful if an automatic quadrature routine has to be implemented. For larger N the magnitude of the weights is so small that the resulting plots are not easily readable and thus are not showed. Some comparison results when considering higher degrees are presented in the next subsection. Figs. 7–11 show the calculated rules when the refinable functions reported in the previous subsection are used. 5.2. Convergence simulations In this section we test our procedure in order to evaluate its performances when dealing with high degrees. The tests are made on six functions taken from [4,38], properly translated and rescaled:

F. Calabrò et al. / Applied Numerical Mathematics 90 (2015) 168–189

Fig. 5. Quadrature formulae for the B-spline of degree 2.

Fig. 6. Quadrature formulae for the B-spline of degree 3.

Fig. 7. Quadrature formulae for the Daubechies mask with p = 2.

181

182

F. Calabrò et al. / Applied Numerical Mathematics 90 (2015) 168–189

Fig. 8. Quadrature formulae for the Daubechies mask with p = 6.

Fig. 9. Quadrature formulae for the ripplet with n = 3, h = 5.

Fig. 10. Quadrature formulae for the 4-point scheme with w = 1/5.

F. Calabrò et al. / Applied Numerical Mathematics 90 (2015) 168–189

183

Fig. 11. Quadrature formulae for the 4-point scheme with w = 1/16.

f 1 (x) = x˜ 20 ,

f 2 (x) = e x˜ ,

f 3 (x) = e −˜x ,

f 4 (x) =

2

f 5 (x) = e −1/˜x , 2

1

1 + 16x˜ 2  3 f 6 (x) = x˜ .

, with x˜ := 2

x m

− 1,

In Fig. 12 the test functions are plotted in the case m = 4. First, we consider the B-splines of degree 2 and 3. In Figs. 13–14 we plot the resulting convergence histories: the absolute value of the error in log scale versus the number of nodes. True solutions were calculated in an exact manner using the piecewise polynomial representation of the B-spline functions. In order to compare the available techniques, we have used the methods seen at the beginning of the present section, and compared the results obtained with the Gaussian quadrature12 calculated via the algorithm in [32]. The selected combination reported are: Gaussian, CMP + Cheb, K + Cheb, CMP + OEq, V + Cheb. The first three perform as in the case f (x) dx [38, Figure 3.1]: convergence is guaranteed and, in particular, the three families have (almost) the same convergence history when the integrand function is less regular. The other two fail to converge: in the case of equispaced nodes this is due to the oscillating quadrature weights while for matrix inversion this is for the numerical instability of the procedure used to compute the weights. From the simulations we can conclude that the procedure for the calculation of weights that we have introduced performs well in all examined cases. The calculated errors decrease up to machine precision and a saturation is observed when this precision is reached. Our simulations also put in evidence the big trouble of the quadrature on equally spaced nodes. In particular, the Chebychev-type distribution of nodes is the easiest to choose between those that result effective in our simulations. To make a better comparison and further test the performance of the method, in Tables 2–3 the results of integration with respect to some refinable functions and in two test cases are reported up to 103 quadrature points. In particular, in the tests we use f 1 and the oscillating function f 7 (x) = sin(20π x/m) too. 6. Conclusion In this paper we have addressed the problem of computing quadrature rules with given nodes for integrals involving a compactly supported refinable function. Modified moments can be used in known methods for the calculation of the required quadrature formulae. We have proposed a new method in order to solve the problem. The method does not require inversion of Vandermonde like matrices. It is deeply based on the refinement properties (1.1) and just needs the mask coefficients a to construct the required formulae. The algorithm we have proposed is an extension of those in [30,32], where the case of Gaussian quadrature has been considered. Therefore it keeps the same computational cost. From the performed tests, one can conjecture that also the stability properties are the same. The here-proposed algorithm gives very good results from the stability viewpoint also for high degrees, where matrix inversion fails. It can be compared with Kautsky’s method [13,28]: as testified by several numerical experiments the two approaches behave in a completely similar way.

12

Gaussian quadrature rules can be calculated in these cases because B-splines are positive.

184

F. Calabrò et al. / Applied Numerical Mathematics 90 (2015) 168–189

Fig. 12. Functions used in order to test quadrature rules: see Figs. 13–14.

By the other side, when considering interpolation-based high-degree quadrature rules one has to carefully choose the nodes, otherwise simple composite rules are to be preferred. This is evident in the tests when equispaced points are considered. On the other hand, the numerical results in Section 5 show that the sequence of quadrature rules on Chebychev nodes perform as well as Gauss quadrature – as previously stated in [4,38] – also when applied in the considered cases.

F. Calabrò et al. / Applied Numerical Mathematics 90 (2015) 168–189

185

Fig. 13. Convergence history in semilog scale: the case of B-spline of degree 2. Five methods are compared. Results for Gaussian quadratures calculated via the algorithm proposed in [32] are plotted with blue circles. The ones for interpolating quadratures on Chebychev nodes (5.2) calculated via our procedure are in red stars. Black circles depict the results provided by the interpolating quadratures for the same nodes (5.2) constructed by the procedure in [13]. Results for the interpolating quadratures on equispaced nodes (5.1) calculated via our procedure are reported with green bullets. Finally, the ones for the interpolating quadrature on Chebychev nodes calculated via matrix inversion are plotted by magenta stars. In each panel, the logarithm of the absolute error is plotted versus the number of taken nodes. The six plots refer, in the same order, to the six test functions as in Fig. 12. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

186

F. Calabrò et al. / Applied Numerical Mathematics 90 (2015) 168–189

Fig. 14. Convergence history in semilog scale: the case of B-spline of degree 3. For the description of used methods we refer to Fig. 13.

As we have already discussed in the Introduction, our procedure can be of interest for several relevant applications ranging from wavelet analysis to numerical approximation of solutions of PDE’s (Isogeometric Analysis). The results obtained on the quadrature rules based on Chebychev type nodes are particularly promising for the development of a stable quadrature rules for Isogeometric Analysis. We end by mentioning some possible interesting generalizations of the presented approach.

F. Calabrò et al. / Applied Numerical Mathematics 90 (2015) 168–189

187

Table 2 Absolute errors when calculating integrals with respect to the cardinal B-spline of degree 2. On the columns the three considered quadrature rules. In the upper panel the function is f 1 (x) = (2x/3 − 1)20 , in the lower f 7 (x) = sin(20π x/3). Reference solution for the first function was calculated using the piecewise-polynomial representation for B-splines, while the second considered function integrates to zero because it is centrally anti-symmetric while the B-spline is centrally symmetric. The obtained results, that are in all cases in the order of the machine precision, are thus a test of exactness (upper panel) and symmetry on calculated weights (lower panel). In both cases the proposed method performs as well as (up to machine precision) Mantica’s algorithm, see [32], as can be noticed comparing the first two columns. The third confirms that the Chebychev nodes also can be profitably used instead of the Gauss ones. N +1

Gaussian points + procedure in [32]

Gaussian points + our algorithm

Chebychev points + our algorithm

7 15 31 63 127 255 511 1000

7.72374865473971e−05 2.16840434497101e−19 5.42101086242752e−19 4.66206934168767e−18 1.72388145425195e−17 1.46367293285543e−17 3.57786716920217e−18 1.35525271560688e−17

7.72374865473962e−05 3.15502832193282e−17 1.28044276570538e−16 4.39101879856629e−17 1.77917576504871e−16 8.56519716263549e−17 8.44918753017954e−16 7.08634539936526e−16

0.000187088513845203 3.82032777754954e−07 3.79470760369927e−17 1.11672823766007e−17 3.67110855603592e−16 2.98914538954254e−16 5.68989300120393e−16 1.26092712660064e−16

7 15 31 63 127 255 511 1000

4.86763407359092e−15 6.92294528994997e−14 1.36608330238695e−15 7.57419239579057e−16 5.31847351505061e−15 1.62271219537889e−15 3.83324384424749e−15 1.15701036834310e−14

5.00641195166907e−15 6.89414616974332e−14 1.46790598000933e−15 6.90423721999807e−16 7.54279190661789e−17 5.11785207464186e−16 6.91496257153481e−16 1.13209756950807e−15

1.48427277413266e−15 1.78654495168983e−16 1.05180353944739e−15 8.30477301851773e−16 6.36886748853670e−17 4.28711311828816e−16 1.05935236610071e−15 1.24413342370742e−15

Table 3 Calculated integrals when applying our interpolating quadrature rule on Chebychev nodes with respect to three different refinable functions. In the first column the results are calculated with respect to the Daubechies refinable function with p = 2 as plotted in Fig. 2 on the left. In the second column the refinable function considered is the GP refinable ripplet plotted in Fig. 3. In the third column integration is made w.r.t. the refinable function associated to the 4-points scheme with w = 1/16, see Fig. 4 on the right. In the upper panel the integrated function is f 1 (x) = (2x/m − 1)20 , in the lower f 7 (x) = sin(20π x/m). The exact integral is available for the second case when a centrally symmetric refinable function is considered because, as before, the exact value is zero. The only case where the algorithm fails to give a result is for the 4-points scheme at 1000 nodes. N +1

Daubechies mask with p = 2

GP ripplet with n = 3, h = 5

7 15 31 63 127 255 511 1000

0.0179481003318603 0.0191114518270956 0.0191098163293648 0.0191098163293647 0.0191098163293647 0.0191098163293654 0.0191098163293647 0.0191098163293649

0.000649031570979579 5.58163604431815e−06 6.71666556451113e−06 6.71666556454636e−06 6.71666556456796e−06 6.71666556396310e−06 6.71666556342346e−06 6.71666556396273e−06

7 15 31 63 127 255 511 1000

0.477305932022489 0.459348598600497 0.0228074635659197 0.0331652996085855 0.0331652996085844 0.0331652996085855 0.0331652996085871 0.0331652996085832

−4.92010945873922e−16 1.53668717290875e−17 −1.62767739763161e−15 −1.32304642449846e−15 −1.10060516035535e−15 −8.70100847740203e−17 1.20839920471711e−15 7.34680863286460e−16

4-point scheme with w = 1/16 0.00131800055573950 7.09896217793944e−06 4.45523895983711e−06 4.45523895977289e−06 4.45523895959817e−06 4.45523895948272e−06 4.45523895904834e−06 –

−4.65382940517678e−15 −1.59710772710230e−15 −1.91694963776912e−15 −1.28198269820252e−15 6.14908327391952e−16 2.19938513762878e−15 −7.76191012503920e−16 –

6.1. Dilation greater than 2 A more general form of (1.1) is as follows, see [10]:

φ(x) =



(M )

aj

φ( Mx − j ),

j ∈Z

where the dilation factor M is an integer greater than 2. In this case (1.14) is replaced by



(M )

aj

= M.

j ∈Z

It is clear that the arguments in Sections 3 and 2 can be immediately extended to the case of (compactly supported) refinable functions with general dilation. Relevant examples of refinable functions with dilation M > 2 are:

188

F. Calabrò et al. / Applied Numerical Mathematics 90 (2015) 168–189

(3)

(3)

Fig. 15. Graphs of φ1,h for h = 1/5 (left) and of φ2,h for h = 7/12 (right).

(M )

• the cardinal B-splines of degree n which are refinable for any integer dilation M ≥ 2. The mask coefficients {ak,n , k = 0, . . . , ( M − 1)(n + 1)} can be obtained by the equality ( M − 1)(n+1)

(M )

ak,n zk =

k =0

1  Mn

1 + z + z 2 + . . . + z M −1

n+1

.

For M = 3 the previous relation provides, see [21], (3 ) ak,n

=

k 

1 3n



n+1

j =[(k+1)/2]

j





j

k− j

,

k = 0, . . . , 2(n + 1);

(6.1)

• the refinable ripplets with dilation M = 3 introduced in [21]. The explicit expression of the associated masks is given by (3 )



(3 )



an,h = ak,n,h , j = 0, . . . , 2(n + 1) , n = 1, . . . , h (3 ) 3 − 2h (3) h (3 ) (3 ) ak,n,h := ak,n−1 + ak−1,n−1 + ak−2,n−1 , 3 3 3 (3)

(3)

where ar ,n−1 are the B-spline mask coefficients given in (6.1) – we assume ar ,n−1 = 0 if r < 0 or r > 2n – and 0 < h < (3) 3/2 is a (real) shape parameter. The corresponding 3-refinable functions φn,h are nonnegative, have support [0, n + 1], are centrally symmetric and belong to C n−1 . In Fig. 15 two examples of these refinable ripplets are displayed. 6.2. Non-compactly supported masks In the literature there are some examples of refinable functions associated with masks with infinite support, see, for instance, the fundamental interpolatory spline functions [6] or the fractional splines [40]. Even if in this case the refinable function has infinite support, nevertheless suitable decay to infinity is guaranteed under suitable decaying conditions on the mask coefficients. As a consequence, integrals of type (1.4) are well-defined and our procedure can be also used in this case providing that the mask is truncated discarding the coefficients with sufficiently small values. References [1] F. Auricchio, F. Calabrò, T. Hughes, A. Reali, G. Sangalli, A simple algorithm for obtaining nearly optimal quadrature rules for NURBS-based isogeometric analysis, Comput. Methods Appl. Mech. Eng. 249 (2012) 15–27. [2] A. Barinka, T. Barsch, S. Dahlke, M. Konik, Some remarks on quadrature formulas for refinable functions and wavelets, Z. Angew. Math. Mech. 81 (2001) 839–855. [3] C. de Boor, A Practical Guide to Splines, revised edition, Appl. Math. Sci., vol. 27, Springer-Verlag, New York, 2001. [4] F. Calabrò, A. Corbo Esposito, An evaluation of Clenshaw–Curtis quadrature rule for integration w.r.t. singular measures, J. Comput. Appl. Math. 229 (2009) 120–128. [5] F. Calabrò, C. Manni, The choice of quadrature in NURBS-based isogeometric analysis, in: M. Papadrakakis, M. Kojic, I. Tuncer (Eds.), Proceedings SEECCM III, 2013.

F. Calabrò et al. / Applied Numerical Mathematics 90 (2015) 168–189

[6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40]

189

C. Chui, An Introduction to Wavelets, Wavelet Analysis and Its Applications, vol. 1, Academic Press Inc., Boston, MA, 1992. J. Cottrell, T. Hughes, Y. Bazilevs, Isogeometric Analysis. Toward Integration of CAD and FEA, Wiley, 2009. W. Dahmen, C. Micchelli, Using the refinement equation for evaluating integrals of wavelets, SIAM J. Numer. Anal. 30 (1993) 507–537. I. Daubechies, Ten Lectures on Wavelets, CBMS-NSF Regional Conference Series in Applied Mathematics, vol. 61, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 1992. I. Daubechies, J. Lagarias, Two-scale difference equations. I. Existence and global regularity of solutions, SIAM J. Math. Anal. 22 (1991) 1388–1410. S. Dubuc, Interpolation through an iterative scheme, J. Math. Anal. Appl. 114 (1986) 185–204. N. Dyn, D. Levin, Subdivision schemes in geometric modelling, Acta Numer. 11 (2002) 73–144. S. Elhay, J. Kautsky, Algorithm 655: IQPACK: FORTRAN subroutines for the weights of interpolatory quadratures, ACM Trans. Math. Softw. 13 (1987) 399–415. H.J. Fischer, On generating orthogonal polynomials for discrete measures, Z. Anal. Anwend. 17 (1998) 183–205. G. Galimberti, V. Pereyra, Solving confluent Vandermonde systems of Hermite type, Numer. Math. 18 (1971) 44–60. W. Gautschi, Orthogonal Polynomials: Computation and Approximation, Numerical Mathematics and Scientific Computation, Oxford Science Publications, Oxford University Press, New York, 2004. W. Gautschi, L. Gori, F. Pitolli, Gauss quadrature for refinable weight functions, Appl. Comput. Harmon. Anal. 8 (2000) 249–257. T. Goodman, Total positivity and the shape of curves, in: Total Positivity and Its Applications, Jaca, 1994, in: Mathematics and Its Applications, vol. 359, Kluwer Acad. Publ., Dordrecht, 1996, pp. 157–186. T. Goodman, C. Micchelli, On refinement equations determined by Pólya frequency sequences, SIAM J. Math. Anal. 23 (1992) 766–784. L. Gori, F. Pitolli, A class of totally positive refinable functions, Rend. Mat. Appl. 20 (2000) 305–322. L. Gori, F. Pitolli, E. Santi, Refinable ripplets with dilation m = 3, Jaen J. Approx. 3 (2011) 173–191. W. Gragg, W. Harrod, The numerically stable reconstruction of Jacobi matrices from spectral data, Numer. Math. 44 (1984) 317–335. T. Hughes, A. Reali, G. Sangalli, Efficient quadrature for NURBS-based isogeometric analysis, Comput. Methods Appl. Mech. Eng. 199 (2010) 301–313. J. Hutchinson, Fractals and self-similarity, Indiana Univ. Math. J. 30 (1981) 713–747. D. Huybrechs, Stable high-order quadrature rules with equidistant points, J. Comput. Appl. Math. 231 (2009) 933–947. D. Huybrechs, S. Vandewalle, Composite quadrature formulae for the approximation of wavelet coefficients of piecewise smooth and singular functions, J. Comput. Appl. Math. 180 (2005) 119–135. B. Johnson, J. Modisette, P. Nordlander, J. Kinsey, Quadrature integration for orthogonal wavelet systems, J. Chem. Phys. 110 (1999) 8309. J. Kautsky, S. Elhay, Calculation of the weights of interpolatory quadratures, Numer. Math. 40 (1982) 407–422. S.G. Kwon, Quadrature formulas for wavelet coefficients, J. Korean Math. Soc. 34 (1997) 911–925. D. Laurie, J. de Villiers, Orthogonal polynomials and Gaussian quadrature for refinable weight functions, Appl. Comput. Harmon. Anal. 17 (2004) 241–258. D. Laurie, J. de Villiers, Orthogonal polynomials for refinable linear functionals, Math. Comput. 75 (2006) 1891–1903. G. Mantica, A stable Stieltjes technique for computing orthogonal polynomials and Jacobi matrices associated with a class of singular measures, Constr. Approx. 12 (1996) 509–530. A. Mantzaflaris, B. Jüttler, Integration by interpolation and look-up for Galerkin-based isogeometric analysis, Comput. Methods Appl. Mech. Eng. 284 (2015) 373–400. J. Phillips, R. Hanson, Gauss quadrature rules with B-spline weight functions, Math. Comput. 28 (1974) 126. A. Quarteroni, R. Sacco, F. Saleri, Numerical Mathematics, second edition, Texts in Applied Mathematics, vol. 37, Springer-Verlag, Berlin, 2007. D. Schillinger, S. Hossain, T. Hughes, Reduced Bézier element quadrature rules for quadratic and cubic splines in isogeometric analysis, Comput. Methods Appl. Mech. Eng. 277 (2014) 1–45. W. Sweldens, R. Piessens, Quadrature formulae and asymptotic error expansions for wavelet approximations of smooth functions, SIAM J. Numer. Anal. 31 (1994) 1240–1264. L. Trefethen, Is Gauss quadrature better than Clenshaw–Curtis?, SIAM Rev. 50 (2008) 67–87. ´ One point quadrature rule with cardinal B-spline, J. Math. Anal. Appl. 360 (2009) 432–438. Z. Udoviˇcic, M. Unser, T. Blu, Fractional splines and wavelets, SIAM Rev. 42 (2000) 43–67.