• No results found

A continuous wavelet-Galerkin method for the linear wave equation

N/A
N/A
Protected

Academic year: 2021

Share "A continuous wavelet-Galerkin method for the linear wave equation"

Copied!
20
0
0

Loading.... (view fulltext now)

Full text

(1)

FOR THE LINEAR WAVE EQUATION ∗

HOANG NGUYEN† AND GIANPAOLO EVANGELISTA

Abstract. We consider the continuous space-time Galerkin method for the linear second-order wave equation proposed by French and Peterson in 1996. A bottleneck for this approach is how to solve the discrete problems effectively. In this paper, we tackle this bottleneck by essentially employ-ing wavelet bases in space. We show how to decouple the correspondemploy-ing linear system and we prove that the resulting subsystems can be uniformly preconditioned by simple diagonal preconditioners, leading to efficient iterative solutions.

Key words. wave equation, continuous Galerkin, finite element, wavelet, preconditioning AMS subject classifications. 65T60, 65M60, 65F35

1. Introduction. The continuous Galerkin method, or continuous space-time Galerkin method, is a variational technique dealing with evolution problems, using approximation spaces of continuous functions in both space and time. An advantage of this approach is that methods of any order of accuracy are relatively easy to for-mulate. For wave propagation problems, this approach is particularly appropriate, since the discrete problems inherit energy conservation properties from the continu-ous problems. This is especially useful when dealing with nonlinear wave problems, see [11] and the references therein.

A bottleneck for the practical computation of variational techniques in space-time, in particular the continuous Galerkin method for solving the wave equation, is that, except when the accuracy in time is of the lowest order, the resulting linear system is, in general, a coupled system, and it is not obvious how to solve it efficiently. This is not a problem when one invokes the method of lines with spatial Galerkin discretization, leading e.g. to the celebrated Runge-Kutta Discontinuous Galerkin schemes, see [4] and the references therein. One drawback of such approaches is that the Courant-Friedrichs-Levy (CFL) condition has to be taken into account, cf. [14]. Recently, (semi) explicit discontinuous Galerkin methods have been proposed in [9] and [14]. In those latter approaches, by generating appropriate space-time meshes in a time-stepping manner, one obtains merely small-sized uncoupled linear subsystems, while the CFL condition is effectively avoided. However, the mesh generation would pose some difficulties when the spatial dimension becomes, e.g., larger than two.

In this paper, we reconsider the continuous space-time Galerkin method pro-posed in [11], tackling the bottleneck mentioned above. In particular, we show how the linear system can be decoupled and how the resulting subsystems can be uni-formly preconditioned, meaning that the condition numbers of the preconditioned subsystems are bounded independently on both spatial and temporal stepsizes. The decoupling merely relies on decompositions of temporal system matrices, while the uniform preconditioning mainly relies on the compressibility of operators in wavelet coordinates and the Riesz properties of wavelet bases (the latter means that wavelets generate Riesz bases for a scale of Sobolev spaces). Moreover, all of the resulting preconditioned subsystems are symmetric positive definite (SPD), leading to efficient

Received . . . ; . . .

Department of Science and Technology, Link¨oping University, Campus Norrk¨oping, SE-60174

Norrk¨oping, Sweden (hoang@itn.liu.se, giaev@itn.liu.se). 1

(2)

iterative solutions in the sense that the computational complexity is of linear order with respect to the number of degrees of freedom. Hence, our method can be viewed as a semi explicit method as well.

Although non-wavelet bases, in combination with preconditioning techniques such as those discussed in [1] and the references therein, might also offer uniformly well conditioned (sub)systems, we have opted for wavelet bases because of the following. Firstly, wavelet bases offer diagonal preconditioners in many applications, thanks to their Riesz properties (see e.g. [6, 5], and [15] for a related problem). It turns out that we are able to obtain diagonal preconditioners for the wave problem also, and our numerical results indicate satisfactory performances. Secondly, sparse wavelet approximations to parabolic problems have recently been proposed in [19], dramat-ically reducing the number of degrees of freedom while retaining the order of accu-racy. Finally, wavelet bases have successfully been used in adaptive approximations to operator equations, such as partial differential equations or (boundary) integral equations, see e.g. [6, 5]. Hence, this paper can also be viewed as an intermediate step toward other wavelet applications to the wave equation, which will be objects of future work. It should be noted that although we focus on nodal bases for the temporal discretization, our approach would apply to other temporal bases as well. Moreover, our approach would also work well with other approaches such as the dis-continuous Galerkin methods for second-order hyperbolic problems proposed in [13] or other Galerkin methods for parabolic problems.

The outline of this paper is as follows: In the rest of this section, we specify our notation. In Section 2, we formulate the wave propagation problem and we recall the continuous Galerkin approximation to its solution as well as the discretization error estimates given in [11]. In addition, we give a matrix vector representation of the re-sulting coupled linear system. In Section 3, we discuss the concept of multiresolution approximation and biorthogonal wavelets. In addition, we prove some spectral prop-erties of the representation of operators in wavelet coordinates, In Section 4, we show how the coupled linear system from Section 2 can be decoupled and preconditioned, employing the results of Section 3. Our numerical results are presented in Section 5. We begin with some basic notation and definitions that will be used throughout this paper. First of all, the time derivative will be denoted by the dot notation. To make the notation not unnecessarily complicated, we will drop references to the un-derlying domain or index sets whenever there is no risk of confusion, i.e., we will write L2 for L2(Ω) etc.. Unless otherwise stated, h·, ·i and k · k will denote some canonical

inner product and norm, e.g., the Euclidean inner product, the spectral norm, the L2-norm, etc.. In order to avoid repeated use of generic but unspecified constants, by C . D we mean that C can be bounded by a multiple of D independently of parameters on which C and D may depend. In cases of interest, we will state the (in)dependence explicitly. Obviously, C & D means D . C, and C h D means C. Dand C & D.

We will adopt the following compact notation from the literature (cf. [6]). For Σ being a countable collection of functions in some separable Hilbert space H, equipped with some inner product h·, ·iH and norm k · kH, we will formally identify Σ with a

column vector (of functions in H): for c = (cσ)σ∈Σ being a column vector of scalars,

cTΣ will denote the formal seriesP

σ∈Σcσσ; and likewise for C = (cς,σ)ς,σ∈Σ being

a matrix, CΣ will denote the column vector (P

σ∈Σcς,σσ)ς∈Σ. Such a collection Σ is

called a Riesz system (in H) if kcTΣk

(3)

where l2(Σ) := c = (cσ)σ∈Σ : kck = kckl2(Σ) :=

P

σ∈Σc2σ

12

<∞ . In addition, if such a system is a basis for H then it is called a Riesz basis. Likewise, for Σ and

˜

Σ being two countable collections of functions in H, hΣ, ˜ΣiH will denote the matrix

(hσ, ˜σiH)σ∈Σ,˜σ∈ ˜Σ, and so, for A and ˜Abeing two matrices of appropriate dimensions,

hAΣ, ˜A ˜ΣiH= AhΣ, ˜ΣiHA˜T.

Finally, for s ≥ 0, Hs and Hs

0 will denote the usual Sobolev spaces on a

d-dimensional domain Ω, whereas for s < 0, Hs and Hs

0 will denote the dual of H−s

and H0−srespectively, i.e. Hs= (H−s)′ and H0s= (H0−s)′, with (H0)′= H0= H00=

(H00)′= L2being the pivot space. For H being any of the Sobolev spaces above, h·, ·iH

and k · kH will denote the usual inner product and the usual norm on H respectively.

2. Continuous Galerkin approximation and error estimates. In this sec-tion, we slightly simplify the formulation and error estimates in [11] so that they bet-ter fit our setting. The coupled linear system resulting from the continuous Galerkin approach will also be given.

2.1. The linear wave equation. Let IR ∋ T > 0 and let Ω be an open, bounded d-dimensional domain with sufficiently smooth boundary ∂Ω. We consider the following initial boundary value problem: Given f = f (t, x), U0 = U0(x) and

V0= V0(x), find U = U (t, x) such that:

¨ U− ∆U = f in (0, T ) × Ω, U = 0 on (0, T ) × ∂Ω, U(0, ·) = U0 in Ω, ˙ U(0, ·) = V0 in Ω, (2.1)

with ∆ being the Laplacian operator. Note that, with V = ˙U, Y = U V  , Y0= U0 V0  , F = 0 f  , and A =  0 −I −∆ 0  , (2.1) is equivalent to ˙ Y + AY = F in (0, T ) × Ω, Y = 0 on (0, T ) × ∂Ω, Y(0, ·) = Y0 in Ω. (2.2)

Under standard regularity assumptions on the given data, results on the existence, uniqueness and regularity of the solution of (2.1) are known. Further, with A being considered as an operator in H1

0 × L2 with domain D(A) = (H2 ∩ H01) × H01, A

generates a strongly continuous semigroup, see e.g. [17]. The latter is essential for the derivation of the results in the next section. In this paper, we tacitly assume that the given data are sufficiently regular so that the corresponding results or formulations are meaningful, cf. [11]. As mentioned in [11], although we confine ourselves to (2.1), the results therein can be easily generalized to the case where −∆ is replaced by any time-independent, selfadjoint, second-order uniformly elliptic operator. Our results can also be generalized to that case.

2.2. Continuous Galerkin approximation. Continuous Galerkin approxima-tions to the solution of (2.1) are defined in the tensor product of the temporal and spatial approximation spaces as follows. For the spatial discretization, with j and p being non-negative discretization parameters, the approximation space

(4)

is assumed to satisfy the following direct estimates, or Jackson estimates: inf

vj∈Vj

kw − vjkL2 .2−jskwkHs (2.3)

for w ∈ Hs∩ H1

0 and 0 ≤ s ≤ p + 1. In addition, with Px : H01 → Vj being the

Galerkin projection from H1

0 onto Vj defined by

h∇(w − Pxw), ∇vjiL2 = 0 ∀vj∈ Vj, (2.4)

we assume that

kw − PxwkHm .2−j(s−m)kwkHs (2.5)

for w ∈ Hs∩ H1

0, m ∈ {0, 1} and m ≤ s ≤ p + 1. For the temporal discretization,

with Pq(I) being the space of polynomials of degree at most q on some interval I, let

Vh= Vh(q) := {vh= vh(t) : t ∈ [0, T ], vh∈ C([0, T ]), (vh)|In ∈ Pq(In)}

where 0 = t0< t1<· · · < tN = T is a partition of (0, T ) into intervals In:= (tn−1, tn)

of length hn := tn− tn−1, and

h:= max

n hn. (2.6)

Throughout this paper, we will assume that p, q ≥ 1. Based on the formulation (2.2), with

Sn:= Ω × In, Vhj := Vh⊗ Vj, u0:= Px(U0) and v0:= Px(V0),

the approximations u and v to U and V respectively are the solution of the following variational problem: Find u, v ∈ Vhj such that u(0, ·) = u0, v(0, ·) = v0, and

h ˙u − v, σiL2(S

n)= 0 ∀σ ∈ Pq−1(In) ⊗ Vj,

h ˙v, ςiL2(S

n)+ h∇u, ∇ςiL2(Sn)= hf, ςiL2(Sn) ∀ς ∈ Pq−1(In) ⊗ Vj

(2.7) for n = 1, · · · , N . As we will see in Section 4 how u0 and v0 can be efficiently

approximated, we simply assume here and in the following that u0 and v0 are given.

Note that the test functions are of one degree lower in time, namely (q − 1), which takes the continuity of u and v at the temporal nodes into account, allowing for u and v to be computed successively on each space-time slab Sn as follows. Let un := u|In

and vn := v|In, then (2.7) is equivalent to the following variational problem: With

u0(t, ·) := u0, v0(t, ·) := v0∀t, find un, vn∈ Pq(In) ⊗ Vj for n = 1, · · · , N such that:

h ˙un− vn, σiL2(S

n)= 0 ∀σ ∈ Pq−1(In) ⊗ Vj,

h ˙vn, ςiL2(S

n)+h∇un,∇ςiL2(Sn)= hf, ςiL2(Sn) ∀ς ∈ Pq−1(In) ⊗ Vj,

un(tn−1,·) = un−1(tn−1,·),

vn(tn−1,·) = vn−1(tn−1,·).

(2.8)

The existence and uniqueness of the solution u and v of (2.7) are proven in [11]. In addition, u and v satisfy the following (simplified) error estimates:

Theorem 2.1 ([11]). Let U and V be the solution of (2.1), and let u and v be their approximations defined in (2.7). For t ∈ [0, T ], there holds

(i) ku(t) − U (t)kL2, kv(t) − V (t)kL2. hq+1+ 2−j(p+1),

(ii) ku(t) − U (t)kH1

0, kv(t) − V (t)kH10 . h

(5)

In addition, for 1 ≤ n ≤ N , there holds

(iii) ku(tn) − U (tn)kL2, kv(tn) − V (tn)kL2 . h2q+ 2−j(p+1),

(iv) ku(tn) − U (tn)kH1

0, kv(tn) − V (tn)kH01 . h

2q+ 2−jp.

Remark 2.2. The statements of Theorem 2.1 are meaningful only when the given data f , U0 and V0, and so the exact solution U and V , are sufficiently regular.

Further, the unspecified constants therein depend on the given data, the exact solution and the final time T (see [11] for details). In addition, in order to obtain the predicted rates of convergence in an optimal manner, hn and j should relate to each other in

one or other way. For example, the optimal L2 rate of convergence is obtained with

either hnh 2−

p+1 q+1j or h

n h 2−

p+1

2q j, as followed from parts (i) and (iii). 

Remark 2.3. Suppose that Vjis a C0Lagrange finite element space incorporating

Dirichlet boundary condition, consisting of continuous piecewise polynomials of degree p with respect to some mesh consisting of regular elements with diameters of order 2−j. Then it is well-known that, by interpolation between Sobolev spaces, Jackson

estimates (2.3) are satisfied. In addition, by the regularity of the elliptic problem (2.4), the estimates (2.5) are also satisfied (see e.g. [3, 2]). Further, not only estimates of type (2.3) and (2.5), but also several inverse estimates are satisfied by Pq(In), and so

by Vh, as well (see [11, 10] for details). All of these estimates are essential for the

derivation of Theorem 2.1. 

2.3. Matrix vector representation. Recall that h·, ·i will denote some canon-ical inner product (in this section, either h·, ·iL2(Sn), or h·, ·iL2(In), or h·, ·iL2(Ω)). For

ˆ

q≥ 0, let Ξqˆ

n = (ξn,0,· · · , ξn,ˆq)T be a basis for Pqˆ(In), then we write Pqˆ(In) ⊗ Vj ∋

wn =Pqi=0ˆ ξn,iwn,i with wn,i ∈ Vj. In addition, let Φj be a basis for Vj, then we

write Vj∋ wn,i= wn,iT Φj, and so, with wn:= (wTn,0,· · · , wTn,ˆq)T,

wn = wTn(Ξ ˆ q n⊗ Φj)

is the representation of wn with respect to the basis Ξqnˆ⊗ Φj for Pqˆ(In) ⊗ Vj, where

Ξqˆ

n⊗Φj= (ξn,0(Φj)T,· · · , ξn,ˆq(Φj)T). Note that, with Ijbeing the |Φj|×|Φj| identity

matrix where |Φj| is the cardinal of Φj (equal to the dimension of Vj), it holds that

wTn(Ξqˆ

n⊗ Φj) = wnT[Ξqnˆ⊗ Ij]Φj. (2.9)

In this paper, we write un = uTn(Ξqn⊗ Φj) and vn = vTn(Ξqn⊗ Φj). Employing

(2.9), (2.8) is represented as huT n( ˙Ξqn⊗ Φj) − vTn(Ξqn⊗ Φj), Ξq−1n ⊗ Φji = 0, hvT n( ˙Ξqn⊗ Φj), Ξq−1n ⊗ Φji + huTn(Ξqn⊗ ∇Φj), Ξq−1n ⊗ ∇Φji = fnT, uTn[Ξq n(tn−1) ⊗ Ij]Φj= uT[n]Φj, vTn[Ξq n(tn−1) ⊗ Ij]Φj= vT[n]Φj,

with the right hand side

fnT := hf, Ξq−1 n ⊗ Φji, and uT[1]Φj:= u0 and vT[1]Φj:= v0, whereas for n = 2, · · · , N uT[n]:= uT n−1[Ξqn−1(tn−1) ⊗ Ij] and vT[n]:= vn−1T [Ξqn−1(tn−1) ⊗ Ij],

(6)

or equivalently [hΞq−1 n , ˙Ξqni ⊗ hΦj,Φji]un− [hΞq−1n ,Ξqni ⊗ hΦj,Φji]vn= 0, [hΞq−1 n ,Ξqni ⊗ h∇Φj,∇Φji]un+ [hΞnq−1, ˙Ξqni ⊗ hΦj,Φji]vn= fn, [(Ξq n(tn−1))T ⊗ Ij]un = u[n], [(Ξq n(tn−1))T ⊗ Ij]vn= v[n], (2.10)

Note that the dimension of this coupled linear system is 2(q + 1)|Φj|. Further, since

(2.7) has a unique solution, (2.10) also has a unique solution, independently on the choice of bases.

3. Wavelet approximation and representation. In this section, we briefly discuss the concept of multiresolution approximation, that will be our spatial ap-proximation, and biorthogonal wavelets. In addition, based on the compressibility of operators in wavelet coordinates and the Riesz properties of wavelet bases, we prove some spectral properties of the representation of operators in wavelet coordinates.

3.1. Multiresolution approximation and wavelet bases. In the following, let Hs:= H0sfor 0 ≤ s ≤ 1 and Hs:= Hs∩ H01⊂ Hsfor s > 1, whereas Hs:= (H−s)′

for s < 0. We recall the following result:

Theorem 3.1 ([8]). Consider the following two multiresolution analyses V0⊂ V1⊂ · · · ⊂ L2, with closL2(∪j≥0Vj) = L2

˜

V0⊂ ˜V1⊂ · · · ⊂ L2, with closL2(∪j≥0j) = L2.

Suppose that

(Q) ∃ uniformly bounded biorthogonal projectors Qj: L2→ L2 such that

Im Qj= Vj, Im(I − Qj) = ˜V ⊥L2

j .

(J) both sequences satisfy Jackson estimates with parameters p + 1 > 0, ˜p+ 1 > 0 uniformly in j, i.e., inf vj∈Vj kv − vjkL2 .2−j(p+1)kvkHp+1 (v ∈ Hp+1), inf ˜ vj∈ ˜Vj k˜v− ˜vjkL2 .2−j( ˜p+1)k˜vkHp+1˜ (˜v∈ Hp+1˜ ).

(B) both sequences satisfy Bernstein estimates with parameters 0 < γ < p + 1, 0 < ˜γ <p˜+ 1 uniformly in j, i.e., for s ∈ [0, γ) and ˜s∈ [0, ˜γ), it holds that

kvjkHs .2 jskv jkL2 (vj ∈ Vj), k˜vjkHs˜.2 j˜sv jkL2 (˜vj ∈ ˜Vj).

Then, for s ∈ (−(˜p+ 1), γ) and ˜s∈ (−˜γ, p+ 1), with Q−1:= 0,

(R1)            k ∞ X j=0 wjk2Hs. ∞ X j=0 4jskwjk2L2 (wj∈ Im(Qj− Qj−1)) ∞ X j=0 4j˜sk(Qj− Qj−1)uk2L2.kuk2H˜ s (u ∈ Hs˜).

(7)

For s ∈ (−˜γ, γ), the mappings (wj)j 7→P∞j=−1wj and u 7→ ((Qj− Qj−1)u)j, which

are bounded in the sense of (R1), are each others inverse. Thus, for s ∈ (−˜γ, γ), (R2) kuk2 Hsh ∞ X j=0 4jsk(Q j− Qj−1)uk2L2 (u ∈ Hs).

Analogous results (R1∗) and (R2∗) are valid at the dual side, i.e., with interchanged roles of (Qj, d, γ) and (Q∗j, ˜d,γ).˜

Now, for j > 0 and Ij being some appropriate index set, let Ψj = (ψj,z)z∈Ij,

whose elements are called (biorthogonal) wavelets, be uniform L2-Riesz bases for the

detail spaces

Wj := Vj∩ ˜V ⊥L2

j−1 = Im(Qj− Qj−1).

Further, let Ψ0 be some arbitrary basis for W0 := V0. As a direct consequence of

Theorem 3.1, we obtain the following. Firstly, our spatial approximation spaces Vj

for j ≥ 0, that are nested closed subspaces of L2 = L2(Ω) incorporating Dirichlet

boundary condition (i.e. (vj)|∂Ω = 0 ∀vj∈ Vj), can be decomposed as

Vj = Vj−1⊕ Wj,

with V−1 := {0}. Repeating this two-level decomposition, we obtain the following

multilevel decomposition Vj= j M l=0 Wl.

Secondly, let the following wavelet collection

Ψj := j

[

l=0

Ψl= (ΨT0,· · · , ΨTj)T

be the corresponding wavelet basis, or multilevel basis, for Vj, then the following norm

equivalences hold uniformly (in j)

kcTΨjk2Hsh j X l=0 4lskclk2l2 (3.1) for s ∈ (−˜γ, γ) and c = (cT

0,· · · , cTj)T with cl ∈ IR|Ψl| and l = 0, · · · , j. In the rest

of this paper, c and d will denote vectors of wavelet coefficients, accordingly to this multilevel setting.

It follows from Bernstein estimates that if γ > 1 then the assumption Vj ⊂ H01=

H1 made in Subsection 2.2 is satisfied. In the next two subsections, we will briefly

discuss two types of approximation spaces and their wavelet bases suitable for our purpose. In particular, γ will be larger than 1 and the wavelets will have compact supports. The latter implies that an application of hΨj,Ψji

L2, or h∇Ψj,∇ΨjiL2, to

a vector of appropriate dimension can be carried out within linear complexity. For a general theory of multiresolution approximation and wavelets, see [5, 6].

Remark 3.2. Since (2.3) follows from (J) by interpolation between Sobolev spaces, both (2.3) and (J) are called Jackson estimates. When the spaces Vj

(8)

are spanned by compactly supported functions such that the diameters of their supports are of order 2−j uniformly, Jackson estimates with parameter p + 1 are typically valid

(cf. Remark 2.3). Likewise, when Vj are spanned by piecewise smooth functions from

Cm(Ω), Bernstein estimates with γ = m +3

2 are typically valid (m = −1 means that

they are not globally continue). Further, the norm equivalences (3.1) are also valid in the limit, i.e., kcTΨk2

Hs h P∞ l=04lskclk2l2 with Ψ := limj→∞Ψ j=S∞ l=0Ψl, and thus the collectionsS∞

l=02−lsΨl are in fact Riesz bases for Hs for s ∈ (−˜γ, γ). 

3.2. Finite element wavelets. In the finite element setting, the spaces Vj are

nothing else than C0Lagrange finite element spaces incorporating Dirichlet boundary

condition: With M0 being some fixed d-simplex (the reference element), we first fix

some conforming dyadic refinement of M0 into 2d congruent subsimplices. Now, let

M0be some (initial) collection of d-simplices M such that ∪M∈M0M is a conforming

triangulation, or mesh, of Ω. Starting from M0, we obtain a sequence of conforming

triangulation (Mj)j by defining Mj+1as the triangulation obtained by applying the

above fixed dyadic refinement to all M ∈ Mj (by means of affine bijections). The

spaces Vj are the collections of continuous piecewise polynomials of degree p, with

respect to the triangulations Mj, that vanish at ∂Ω. Employing standard finite

element techniques, compactly supported biorthogonal wavelet bases for Vj, with

˜

γ=32 and basically any desired ˜p+ 1, can be effectively constructed, see e.g. [8, 16]. Since the spaces Vj consist of continuous piecewise polynomials of degree p on

regular elements with diameters of order 2−j, the Jackson and Bernstein parameter

are p + 1 and γ = 3

2 respectively, as mentioned in the previous subsection. Although

this finite element approach has the advantage of inheriting the widely used and well documented finite element techniques, it is not feasible when d becomes larger than three. In this case, one can employ the wavelets discussed in the next subsection.

3.3. Tensor product based spline-wavelets. First, let us consider the simple case when Ω is the unit d-cube [0, 1]d. For d = 1, the spaces V

j are defined through a

basis, generated by dyadic dilations and translations of a single scaling function, that is a B-spline of degree p in our case, with suitable adaptations near the boundary so that Dirichlet boundary condition is satisfied. Compactly supported biorthogonal wavelets, with basically any desired ˜γ and ˜p+ 1, can then be constructed following traditional lines. For d > 1, wavelets on the d-cube are defined by tensor products of wavelets constructed on the interval [0, 1], and so the spaces Vjcan be defined. Now, for general

domains Ω, one might employ domain decomposition techniques as follows: First, Ω is decomposed into disjoint patches, each of them being some parametric image of the d-cube. Then, roughly speaking, wavelets on Ω are constructed by lifting wavelets on the d-cube to the patches and then ‘gluing’ them over the interfaces.

Since the spaces Vj have polynomials reproduction properties, meaning that

poly-nomials of degree p can be reproduced, the Jackson parameter is also p + 1, equal to that in the finite element case. But, since the scaling functions are splines, the Bernstein parameter is now γ = p +12, larger than that in the finite element case. For more information about tensor product based spline-wavelets, see e.g. [7].

3.4. Operators in wavelet coordinates. For r ∈ (−˜γ, γ), let Lr: Hr→ H−r

be a bounded elliptic operator. Here, boundedness means |hv, Lrui| . kvkHrkukHr

and ellipticity means hu, Lrui & kuk2Hr, It follows from the boundedness and ellipticity

of Lr that

(9)

In addition, we assume for simplicity that Lr,L∗r : Hγ → Hγ−2r are bounded, where

L∗

ris the adjoint of Lr, and thus Lr,L∗r: Hˆγ+r→ Hˆγ−rare bounded for any ˆγ= ˆγr∈

[0, γ − r]. For a discussion about sufficient conditions for this continuity property, see e.g. [18] and the references therein. In particular, when Lr is a differential operator

given by its associated bilinear form hv, Lrui =

X

|α|,|β|≤m

h∂βv, a

αβ∂αui,

with α and β being multi-indices, and ∂α and ∂β being the corresponding partial

differential operators, then a sufficient condition for this continuity property is that aαβ∈ C⌈(γ−r)⌉(Ω), with ⌈·⌉ being the ceiling function. By Cauchy-Schwarz inequality

and the norm equivalences (3.1), we then obtain for any ˆγ∈ [0, min{˜γ+ r, γ − r}) and wl∈ Wland wl′ ∈ Wl′ |hwl′,Lrwli| ≤ kwl′kH −ˆγ+rkLrwlkHˆγ−r .2 −ˆγ(l−l′) kwl′kHrkwlkHr, and analogously |hwl′,Lrwli| = |hL∗rwl′, wli| . 2−ˆγ(l ′−l) kwl′kHrkwlkHr, and so |hwl′,Lrwli| . 2−ˆγ|l−l ′| kwl′kHrkwlkHr. (3.2)

In fact, (3.2) is even valid for any ˆγ∈ [0, min{˜p+ 1 + r, γ − r}) as pointed out in [18]. Since h·, Lr·i h k · k2Hr, it follows from (3.1) that, for all j and c = (c

T

0,· · · , cTj)T,

cThΨj,L

rΨjic h cTD2rj c, (3.3)

with Dj being the diagonal matrix defined by Dj = diag((2lI[l])l=0,··· ,j) where I[l]

denotes the |Ψl| × |Ψl| identity matrix, i.e., Djc = (20cT0,· · · , 2jcTj)T. Note that

hΨj,L

rΨji is invertible and is the representation of Lrwith respect to Ψj. Also, (3.2)

can be written as |dT l′hΨl′,LrΨlicl| . 2−ˆγ|l−l ′| 2rl′ kdl′k2rlkclk, (3.4)

for any cl∈ IR|Ψl|, dl′ ∈ IR|Ψl′|and ˆγ∈ [0, min{˜p+ 1 + r, γ − r}).

In the following, with A and B being two matrices, we will write A h B if xTAx h xTBxfor all x of appropriate dimension. In particular, we write (3.3) as

hΨj,L

rΨji h D2rj , (3.5)

which is equivalent to hΨj,L

rΨji−1 h D−2rj . Further, it follows from (3.3) that, for

any d = (dT 0,· · · , dTj)T, sup d |dTD−r j hΨj,LrΨjiD−rj c| kdk ≥ cTD−r j hΨj,LrΨjiD−rj c kck h kck,

whereas it follows from the boundedness of Lrand (3.1) that

|dTD−rj hΨj,LrΨjiD−rj c| ≤ kD −r j d TΨjk HrkD −r j c TΨjk Hr h kdkkck, and so kD−rj hΨj,L

rΨjiD−rj ck h kck for all c, which will be written as

kD−rj hΨj,L

(10)

Lemma 3.3. With Lr being given as above and ˇr∈ IR, it holds that

hLrΨj,ΨjiD−2ˇj rhΨj,LrΨji h D2rD−2ˇrD2rj ,

provided that γ and ˜pare sufficiently large. Proof. Let Lj and L[l

,l]

j denote hΨj,LrΨji and hΨl′,LrΨli for l′, l = 0, · · · , j

respectively. Further, let L[j] denote the block diagonal matrix diag(L [l,l] j )l=0,··· ,j, i.e. L[j]c = ((L [0,0] j c0)T,· · · , (L [j,j]

j cj)T)T with c = (cT0,· · · , cTj)T as before. Note

that the counterparts of (3.5) and (3.6) within any level l are L[l,l]j h 22rlI[l] and kL[l,l]j k h 22rl respectively, with I[l]being the |Ψ

l| × |Ψl| identity matrix as above.

It is sufficient to prove that ∃ ǫ ∈ (0, 1) such that

kD−ˇjr[Lj− L[j]]ck ≤ ǫkD−ˇjrL[j]ck (3.7)

for all c. Indeed, since D−ˇr

j Ljc= D−ˇjrL[j]c+ D−ˇjr[Lj− L[j]]c, having proven (3.7), we infer that kDj−ˇrLjck ≤ kD−ˇjrL[j]ck + kDj−ˇr[Lj− L[j]]ck ≤ (1 + ǫ)kD−ˇjrL[j]ck and that kD−ˇr j Ljck ≥ kD−ˇjrL[j]ck − kDj−ˇr[Lj− L[j]]ck ≥ (1 − ǫ)kD−ˇjrL[j]ck, i.e., kD−ˇr j Ljck h kD−ˇjrL[j]ck.

In addition, it follows from kL[l,l]j k h 22rlthat

kD−ˇjrL[j]ck h kD−ˇjrD2rj ck.

The last two relations yield kD−ˇjrLjck h kD−ˇjrD2rj ck for all c, proving our assertion.

Now, we prove (3.7) as follows. Let d = (dT

0,· · · , dTj)T, then it holds that

|hd, Dj−ˇr[Lj− L[j]]ci| = |hDj−ˇrd,[Lj− L[j]]ci| = j X l=0 j X l6=l′=0 |h2−ˇrd l′, L[l ′,l] j cli|. (3.8)

By (3.4) and kL[l,l]j k h 22rl, it holds for any ˆγ= ˆγ

r∈ [0, min{˜p+ 1 + r, γ − r}) that j X l=0 j X l6=l′=0 |h2−ˇrd l′, L [l′,l] j cli| . j X l=0 j X l6=l′=0 2−ˆγ|l−l′| 2(r−ˇr)l′ kdl′k2rlkclk = j X l=0 j X l6=l′=0 2−ˆγ|l−l′|2−(r−ˇr)(l−l′)kdl′k2(2r−ˇr)lkclk . j X l=0 j X l6=l′=0 2−(ˆγ−|r−ˇr|)|l−l′|kdl′kk2−ˇrlL[l,l]j clk. (3.9)

Let A ∈ IR(j+1)×(j+1)and x, y ∈ IRj+1be defined by Ai,i′ = (1 − δi,i′)2−(ˆγ−|r−ˇr|)|i−i ′|

, and xi = k2−ˇr(i−1)L

[i−1,i−1]

j ci−1k and yi = kdi−1k, then the last double sum is

yTAx. Since yTAx≤ kAkkykkxk ≤ (kAk

1kAk∞)

1

(11)

j X l=0 j X l6=l′=0 |h2−ˇrd l′, L[l ′,l] j cli| . 2−(ˆγ−|r−ˇr|)( j X l′=0 kdl′k2) 1 2( j X l=0 k2−ˇrlL[l,l] j clk2) 1 2 = 2−(ˆγ−|r−ˇr|)kdkkD−ˇrL [j]ck. (3.10)

Since d is arbitrary, it follows from (3.8) and (3.10) that kD−ˇr

j [Lj− L[j]]ck . 2−(ˆγ−|r−ˇr|)kD−ˇjrL[j]ck, (3.11)

yielding (3.7) if γ and ˜p are sufficiently large, since ˆγ ∈ [0, min{˜p+ 1 + r, γ − r}) arbitrary.

In addition, let ˆr∈ (−˜γ, γ), and let Lrˆbe the corresponding operator satisfying

all of the assumptions above with interchanged roles of r and ˆr. Lemma 3.3 is an intermediate step to the following more general result.

Lemma 3.4. With Lr, Lrˆbeing given as above and ˇr∈ IR and α ≥ 0, it holds that

[hLrΨj,Ψji + αhLrˆΨj,Ψji]Dj−2ˇr[hΨj,LrΨji + αhΨj,LrˆΨji]

h [D2r+ αD2ˆr]D−2ˇr[D2rj + αD2ˆr]

independently on α, provided that γ and ˜pare sufficiently large.

Proof. The proof is analogously to and employs the proof of Lemma 3.3. Let Lj, L[j] and L[l

,l]

j be the stiffness, the block diagonal matrix and the block matrices

respectively, defined as in the proof of Lemma 3.3, corresponding to Lr. Likewise,

let ˆLj, ˆL[j] and ˆL [l′,l]

j be those matrices corresponding to Lˆr. Note that within any

level l we have [L[l,l]j + αˆL[l,l]j ] h (22rl+ α22ˆrl)I[l]independently on α, from which we

infer that k[L[l,l]j + αˆL[l,l]j ]clk & (22rl+ α22ˆrl)kclk independently on α. Combining this

with kL[l,l]j k h 22rl and kαˆL[l,l] j k h α22ˆrl, we obtain kL [l,l] j + αˆL [l,l] j k h 22rl+ α22ˆrl independently on α.

It is sufficient to prove that ∃ ǫ ∈ (0, 1) being independent on α such that kD−ˇjr[[Lj+ αˆLj] − [L[j]+ αˆL[j]]]ck ≤ ǫkD−ˇjr[L[j]+ αˆL[j]]ck (3.12)

for all c = (cT

0,· · · , cTj)T. Indeed, having proven (3.12), we infer that

kD−ˇjr[Lj+ αˆLj]ck h kD−ˇjr[L[j]+ αˆL[j]]ck

independently on α. In addition, it follows from kL[l,l]j + αˆL[l,l]j k h 22rl+ α22ˆrl that

kD−ˇjr[L[j]+ αˆL[j]]ck h kD−ˇjr[D 2r

j + αD2ˆr]ck

independently on α. The last two relations yield our assertion. Now, (3.12) is a consequence of (3.11), kL[l,l]j k h 22rl, kαˆL[l,l]

j k h α22ˆrl, and

kL[l,l]j + αˆL[l,l]j k h 22rl+ α22ˆrlas follows. For any ˆγ

r∈ [0, min{˜p+ 1 + r, γ − r}) and

ˆ

γˆr∈ [0, min{˜p+ 1 + ˆr, γ− ˆr}), it holds that

kD−ˇjr[[Lj+ αˆLj] − [L[j]+ αˆL[j]]]ck ≤ kDj−ˇr[Lj− L[j]]ck + kDj−ˇr[αˆLj− αˆL[j]]]ck .2− max{(ˆγr−|r−ˇr|),(ˆγˆr−|ˆr−ˇr|)}(kD−ˇr j L[j]ck+kαD−ˇjrLˆ[j]ck) .2− max{(ˆγr−|r−ˇr|),(ˆγˆr−|ˆr−ˇr|)}(kD−ˇr j D2rj ck+kαD −ˇr j D2ˆjrck) ≤ 2− max{(ˆγr−|r−ˇr|),(ˆγˆr−|ˆr−ˇr|)}2 kD−ˇr j [D2rj + αD2ˆjr]ck h 2− max{(ˆγr−|r−ˇr|),(ˆγˆr−|ˆr−ˇr|)}kD−ˇr j [L[j]+ αˆL[j]]ck

(12)

4. Decoupling and preconditioning. In this section, we show how the cou-pled linear system (2.10) can be decoucou-pled and preconditioned by choosing nodal bases, or Lagrange bases, for the temporal approximation and wavelet bases for the spatial approximation. Although we focus on nodal bases, we will briefly discuss how our approach would also work with other temporal bases.

4.1. Temporal discretization. In the following, we will use the notation de-fined in Section 2. In addition, for ˆq≥ 0, let Ξˆq = (ξ

0,· · · , ξ ˆ q ˆ

q)T be the nodal basis,

consisting of ˆq+ 1 Lagrange polynomials of degree ˆq, for Pˆq(0, 1) with respect to some

fixed partition 0 = τ0qˆ<· · · < τ ˆ q ˆ

q = 1 of (0, 1). For n = 1, · · · , N , with ϕn being the

affine bijections from (0, 1) to In defined by t = ϕn(τ ) = hnτ− tn−1, let Ξqnˆ be the

nodal basis for Pqˆ(In) defined by Ξqnˆ = Ξqˆ◦ ϕ−1n = (ξ ˆ q 0◦ ϕ−1n ,· · · , ξ ˆ q ˆ q◦ ϕ −1 n )T.

Since Ξq(0) = (1, 0, · · · , 0)T, we obtain from the last two subsystems in (2.10)

un,0= u[n],

vn,0= v[n].

(4.1)

Note that, uT

[1]Φj = u0 and vT[1]Φj = v0 as before whereas now, since Ξq(1) =

(0, · · · , 0, 1)T, u

[n] = un−1,q and v[n] = vn−1,q for n = 2, · · · , N . Further, by the

‘pull-back’ via the affine bijections ϕn above, it holds that

hΞq−1n ,Ξqni = hnhΞq−1, Ξqi and hΞq−1n , ˙Ξqni = hΞq−1, ˙Ξqi, (4.2)

with ˙Ξq := ((ξq

0)τ,· · · , (ξqq)τ)T. Let A denote the matrix obtained from hΞq−1, Ξqi

by deleting its first column a, i.e.,

a:= hΞq−1, ξq0i and A := hΞq−1,(ξq1,· · · , ξqq)Ti, (4.3)

and analogously, let

˙a := hΞq−1,(ξq0)τi and ˙A:= hΞq−1,((ξq1)τ,· · · , (ξqq)τ)Ti, (4.4)

then, by employing (4.1) and (4.2), we obtain the following linear system from the first two subsystems in (2.10)

" ˙ A⊗ Mj −hnA⊗ Mj hnA⊗ Sj A˙ ⊗ Mj # u0;n v0;n  = f1;n f2;n  , (4.5)

with the remaining unknowns

u0;n:= (un,1,· · · , un,q)T and v0;n:= (vn,1,· · · , vn,q)T,

the spatial mass- and stiffness matrix

Mj:= hΦj,Φji and Sj:= h∇Φj,∇Φji,

and the right hand side

f1;n:= − ˙a ⊗ Mju[n]+ hna⊗ Mjv[n] and f2;n:= −hna⊗ Sju[n]− ˙a ⊗ Mjv[n]+ fn.

Note that the dimension of the system (4.5) is 2q|Φj|. In the next subsection, we will

decompose this system into 2q systems of dimension |Φj| that yield, together with the

(13)

4.2. Decoupling. In order to decompose (4.5), we note that, besides Mj and

Sj, the two q × q matrices A and ˙Aare also invertible:

Lemma 4.1. For ˆq= q −1, q, let Ξqˆ= (ξˆq 0,· · · , ξ

ˆ q ˆ

q)T be the nodal bases, consisting

of r + 1 Lagrange polynomials of degree ˆq, for Pqˆ(0, 1) with respect to some fixed

partitions 0 = τ0qˆ<· · · < τqˆqˆ= 1 of (0, 1) as before. Then the two q × q matrices A and ˙Adefined in (4.3) and (4.4) respectively, i.e.,

A= hΞq−1,(ξq1,· · · , ξqq)Ti and ˙A= hΞq−1,((ξq1)τ,· · · , (ξqq)τ)Ti

are invertible.

Proof. First, suppose that the rank of A is smaller than q. Then ∃ x = (x1,· · · , xq)T 6= 0 such that Ax = 0. In other words, with p1 := xT(ξq1,· · · , ξqq)T ∈

Pq(0, 1), we have p16= 0 and hΞq−1, p1i = 0, or equivalently we have

0 6= p1⊥ Pq−1(0, 1)

on the one hand. On the other hand, since p1(τ0q) = p1(0) = 0, ∃ 0 6= p2∈ Pq−1(0, 1)

such that p1= τ p2, and so

hp1, p2i = hτ, p22i > 0,

which gives a contradiction. Thus A is invertible.

Next, note that the collection {(ξq0)τ,· · · , (ξqq)τ} spans Pq−1(0, 1). Since Ξq is a

partition of unity, i.e.,Pq

i=0ξ q i = 1, we have Pq i=0(ξ q i)τ= 0. In particular, we have (ξq0)τ = − q X i=1 (ξqi)τ,

from which we conclude that ((ξq1)τ,· · · , (ξqq)τ)T is a basis for Pq−1(0, 1). Thus, since

Ξq−1 is also basis for Pq−1(0, 1), ˙Ais invertible.

By this lemma, we employ a modified Schur complement of the system matrix in (4.5) to obtain the following decomposition of (4.5)

" ˙ A⊗ Mj −hnA⊗ Mj I⊗ Mj+ h2n[ ˙A−1A]2⊗ Sj 0 # u0;n v0;n  = ˜f1;n ˜f2;n ! ,

with ˜f1;nand ˜f2;n being certain functions of f1;nand f2;n. In other words, we obtain

the following q systems of dimension |Φj| for v0;n

[I ⊗ Mj]v0;n= [h−1n A

−1A˙ ⊗ M

j]u0;n− [h−1n A −1⊗ I

j]˜f1;n, (4.6)

with Ij denoting the |Φj| × |Φj| identity matrix as before whereas I = Iq×q denoting

the q × q identity matrix, and with u0;n being given by the following coupled system

of dimension q|Φj|

[I ⊗ Mj+ h2n[ ˙A−1A]2⊗ Sj]u0;n= ˜f2;n. (4.7)

In the following, we show how (4.7) can be decoupled. Let the real Schur decom-position of [ ˙A−1A]2 be given by

(14)

with Q being an orthogonal matrix and T being a real upper quasi-triangular matrix, by which we mean that T is a sum of a real upper triangular matrix T1 with zeros

on the diagonal and a real block diagonal matrix T2with block size of at most 2 × 2.

Note that any element λion the diagonal of T2is a real eigenvalue [ ˙A−1A]2and that

any 2 × 2 matrix Bi on the diagonal of T2 has two complex conjugate eigenvalues

λi,ℜ±iλi,ℑwith λi,ℑ6= 0, that are also two complex conjugate eigenvalues of [ ˙A−1A]2

(see e.g. [12]). We then obtain the following block upper quasi-triangular system from (4.7)

[I ⊗ Mj+ h2nT⊗ Sj]˜u0;n= [QT ⊗ Ij]˜f2;n, (4.9)

with ˜u0;n:= [QT ⊗ Ij]u0;n, i.e., the solution of (4.7) is given by

u0;n= [Q ⊗ Ij]˜u0;n.

Obviously, solving (4.9), and so (4.7), means successively solving either a system of dimension |Φj| of the form

[Mj+ λih2nSj]x = b (4.10)

or a coupled system of dimension 2|Φj| of the form

[I2×2⊗ Mj+ h2nBi⊗ Sj]¯x= ¯b, (4.11)

with I2×2 denoting the 2 × 2 identity matrix.

In order to decouple (4.11), we proceed as follows. First, let pi,ℜ+ ipi,ℑ be a

complex eigenvector corresponding to the complex eigenvalue λi,ℜ+ iλi,ℑ of Bi and

let Pi:= [pi,ℜ pi,ℑ], then Pi has full rank and

Bi= Pi  λi,ℜ λi,ℑ −λi,ℑ λi,ℜ  P−1i . Hence, we obtain the following system from (4.11)

Mj+ λi,ℜh2nSj λi,ℑh2nSj −λi,ℑh2nSj Mj+ λi,ℜh2nSj  x1 x2  = b1 b2  , (4.12) with (bT 1, bT2)T := [P −1 i ⊗ Ij]¯b and (xT1, xT2)T := [P −1

i ⊗ Ij]¯x, i.e., the solution of

(4.11) is given by ¯ x= [Pi⊗ Ij] x1 x2  .

Next, we employ a modified Schur complement to obtain from (4.12) the following equivalent system

"

Mj+ λi,ℜh2nSj λi,ℑh2nSj

Sj+ λ−2i,ℑh−4n [Mj+ λi,ℜh2nSj]S−1j [Mj+ λi,ℜh2nSj] 0

# x1 x2  = ˜b˜1 b2 ! , (4.13)

with ˜b1 and ˜b2being certain functions of b1 and b2. In other words, we obtain from

(4.12) the following two systems of dimension |Φj|

Sjx2= −

1 λi,ℑh2n

(15)

where

[Sj+ λ−2i,ℑh−4n [Mj+ λi,ℜh2nSj]S−1j [Mj+ λi,ℜh2nSj]]x1= ˜b2. (4.15)

In summary, solving (4.9), and so (4.7), means solving q systems of dimension |Φj| either of the form (4.10) or (4.14) or (4.15). Since (2.10) has a unique solution

independently on the choice of bases, any of those q systems has also a unique solution. In particular, we obtain for (4.10):

Lemma 4.2. Let λi∈ IR be an eigenvalue of [ ˙A−1A]2, then λi>0.

Proof. First, note that Mjand Sjare SPD. Let µ be an eigenvalue of S −1 2 j MjS −1 2 j ,

then µ > 0 and ∃ x 6= 0 such that

Mjx− µSjx= 0.

Now, suppose that λi <0. With h2n= − µ

λi >0, we then have

Mjx+ λih2nSjx= 0,

which contradicts the uniqueness of the solution of (4.10). In addition, since [ ˙A−1A]2

is invertible, we conclude that λi>0.

Since q ≪ |Φj| in practical applications, the temporal system matrices A and ˙A

and their inverses, and the Schur decomposition and the eigenvalues of [ ˙A−1A]2, and

the 2 × 2 matrices Pi and their inverses can be computed or approximated with high

precision, as a once and for all preprocessing computation. Thus, there is only one problem left, which will be treated in the next subsection: how to obtain iterative solutions of (4.6), (4.10), (4.14) and (4.15) efficiently.

Remark 4.3. For standard bases Ξqˆ

n for Pqˆ(In) consisting of monomials of

de-gree at most ˆq, since the two corresponding matrices A and ˙A are also invertible, one would obtain the same 2q subsystems of dimension |Φj| as given in (4.6), (4.10),

(4.14) and (4.15) (obviously with different right hand sides and unknowns). For or-thonormal bases Ξqˆ

n for Pqˆ(In) consisting of Legendre polynomials of degree at most

ˆ

q, by the invertibility of the redefined matrix A := hΞq−1,q 0,· · · , ξ

q

q−1)Ti and the

corresponding matrix ˙A= hΞq−1,((ξq

1)τ,· · · , (ξqq)τ)Ti as above, one would then

ob-tain 2(q + 1) subsystems of dimension |Φj| similar to the ones given in (4.6), (4.10),

(4.14) and (4.15). In either case, Lemma 4.2 is still valid. Further, since the ap-proximation spaces and the test spaces are different, we have in fact Petrov-Galerkin approximation, hence one could also choose (suitable) bases of different types for the

approximation spaces and for the test spaces. 

4.3. Preconditioning. In the following, we will drop the subscriptsL2 andl 2,

and we will take Φj = Ψj where Ψj = (ΨT0,· · · , ΨTj)T is some wavelet basis from

Section 3 with parameters p + 1 > γ > 1 and ˜p+ 1 > ˜γ > 0 such that γ and ˜p are sufficiently large. The mass matrix Mj and stiffness matrix Sj defined in Subsection

4.1 are thus

Mj = hΨj,Ψji and Sj= h∇Ψj,∇Ψji.

In terms of Subsection 3.4, Mj is the ‘stiffness’ matrix corresponding to the

iden-tity operator L0u= u and Sj is the stiffness matrix corresponding to the negative

Laplacian operator L1u= −∆u.

Note that the usual H1-seminorm | · |

H1 = k∇ · k is equivalent with the usual

H1-norm k · kH1 = (k∇ · k2+ k · k2) 1 2 on H1

(16)

so that cTS

jc= hcTΨj,−∆cTΨji h kcTΨjk2H1 for any c = (cT0,· · · , cTj)T as before.

Hence, from Section 3, in particular (3.5) with r ∈ {0, 1}, we conclude that Mj h Ij, Sj h D2j and S

−1

j h D

−2

j (4.16)

independently on j, with Ij being the |Ψj| × |Ψj| identity matrix as before and Dj

being the diagonal matrix defined in Subsection 3.4. In other words, Ij and D2j

are (SPD diagonal) preconditioners for Mj and Sj respectively, meaning that the

condition number of Mjand that of the two-sided preconditioned matrix D−1j SjD−1j

are bounded independently on j. Note that the mass- and stiffness matrix Mj and

Sj are not truly sparse. However, as mentioned in Section 3, the applications of

Mj or Sj to vectors of dimension |Ψj| can be carried out with linear complexity.

Since only matrix vector applications matter in the context of iterative methods, we conclude that iterative solutions of (4.6) and (4.14), as well as approximations to the Galerkin projections u0 and v0 of U0 and V0 from Section 2, of any desired

accuracy can be obtained by, e.g., the conjugated gradient method (CG) efficiently, in the sense that the computational complexity is linear, independently on j and hn.

In the following, we propose (SPD diagonal) preconditioners for (4.10) and (4.15), via relations similar to (4.16), such that the condition numbers of the corresponding (two-sided) preconditioned matrices are bounded independently on j and hn, and so

iterative solutions of (4.10) and (4.15) can also be obtained efficiently in this sense. Recall that it follows from Lemma 4.2 that the eigenvalue λiin (4.10) is positive, and

note that applications of S−1j in (4.15) can be iteratively approximated by, e.g., CG with the preconditioner D−2j .

Proposition 4.4. Let Mj, Sj, Ij and Dj be defined as above.

(i) For any ˆh >0, it holds that

[Mj+ ˆhSj] h [Ij+ ˆhD2j]

independently on j and ˆh.

(ii) For any a, b ∈ IR and ˆh >0, it holds that

Sj+ b2ˆh−4[Mj+ aˆh2Sj]S−1j [Mj+ aˆh2Sj]

h D2j+ b2ˆh−4[Ij+ aˆh2D2j]D −2

j [Ij+ aˆh2D2j]

independently on j and ˆh. Moreover, if a ≥ 0 then it also holds independently on a and b.

Proof.

(i) The assertion follows from (4.16).

(ii) If a ≥ 0, the assertion follows from (4.16) and Lemma 3.4. Otherwise, we write the assertion as

Sj+ [bˆh−2Mj− abSj]S−1j [bˆh−2Mj− abSj]

h D2j+ [bˆh−2Ij− abD2j]D−2j [bˆh−2Ij− abD2j]

(4.17)

with a > 0 and b > 0 (for b = 0, the assertion follows from (4.16)). Now, for any c = (cT

0,· · · , cTj)T, it follows from (4.16) that

cT[S

j+ [bˆh−2Mj− abSj]S−1j [bˆh−2Mj− abSj]]c

h kD−1j Sjck2+ kDj−1[bˆh−2Mj− abSj]ck2.

(17)

On the one hand, we have kD−1j Sjck2+ kDj−1[bˆh−2Mj− abSj]ck2 ≤ kD−1j Sjck2+ (kabD−1j Sjck + kbˆh−2D−1j Mjck)2 ≤ kD−1j Sjck2+ 2 kabD−1 j Sjck2+ 2 kbˆh−2D−1j Mjck2 . kD−1j Sjck2+ kbˆh−2D−1j Mjck2 h kDjck2+ kbˆh−2D−1j ck2, (4.19)

where we have used Lemma 3.4 in the last line. Note that the equivalence constant in the fourth line is 2(1 + a2b2), independent on ˆh. Analogously, we

have on the other hand

kD−1j Sjck2+ kDj−1[bˆh−2Mj− abSj]ck2

& kD−1j Sjck2+ kabD−1j Sjck2+ kD−1j [bˆh−2Mj− abSj]ck2

≥ kD−1j Sjck2+ 2−1(kabD−1

j Sjck + kD−1j [bˆh−2Mj− abSj]ck)2

≥ kD−1j Sjck2+ 2−1kbˆh−2D−1j Mjck2

h kDjck2+ kbˆh−2D−1j ck2

(4.20)

(the equivalence constant in the second line is (1 + a2b2)−1, independent on

ˆ

hagain). It follows from (4.18), (4.19) and (4.20) that cT[Sj+[bˆh−2Mj−abSj]S−1j [bˆh

−2M

j−abSj]]c h cT[D2j+b2ˆh−4D −2

j ]c (4.21)

independently on ˆh. With interchanged roles of (Mj, Sj) and (Ij, D2j) in

(4.19) and (4.20), we directly infer that cT[D2j+[bˆh−2Ij−abD2j]D −2 j [bˆh −2I j−abD2j]]c h cT[D 2 j+b 2ˆ h−4D−2j ]c. (4.22) The last two relations yield (4.17) independently on ˆh.

Remark 4.5. For a ≥ 0, b ≥ 0 and ˆh >0, it holds that cT[D2j+ [bˆh−2Ij+ abD2j]D−2j [bˆh

−2I

j+ abD2j]]c h cT[D2j+ b2ˆh−4D−2j ]c, (4.23)

with equivalence constants (1 + a2b2) and (1 + a2b2)−1. Hence, with

P+:= [D2j+ [bˆh−2Ij+ abD2j]D −2 j [bˆh−2Ij+ abD2j]], P−:= [D2j+ [bˆh−2Ij− abD2j]D −2 j [bˆh−2Ij− abD2j]], P0 := [D2j+ b2hˆ−2D −2 j ],

and similarly Q±:= [Sj+[bˆh−2Mj±abSj]S−1j [bˆh−2Mj±abSj]], it follows from (4.22),

(4.23) and the proof of Proposition 4.4 that P+, P− and P0 are preconditioners for

both Q+ and Q−. However, since P+, P−, Q+ and Q− share the same constants

(1 + a2b2) and (1 + a2b2)−1 in their equivalences with P

0, we may expect that P± are

quantitatively better than P0, being considered as preconditioners for Q±. In addition,

since P+ h Q+ independently on a and b, we may even expect that P− h Q−

independently on a and b also. Yet, we do not know how to prove it rigorously.  Remark 4.6. Let DMj and DSj be the diagonal of Mj and Sj respectively, then

DMj h Ij and DSj h D2j. Since DMj and DSj are the diagonal of Mj and Sj,

we may expect that DMj and DSj are quantitatively better than Ij and D2j, being

considered as preconditioners for Mj and Sj respectively. It is well-known that this is

indeed the case, and so we may expect better performances with DMj and DSj instead

of Ij and D2j in Proposition 4.4. In our numerical experiments, we have therefore

(18)

5. Numerical results. In this section, we present the numerical results of sev-eral one and two dimensional experiments. Since the predicted rates of convergence have been numerically observed in [11], we focus on preconditioning, by approximating the l2-condition numbers of the system matrices preconditioned accordingly to

Sub-section 4.3 and Remark 4.6. We only mention that the predicted rates of convergence have also been observed in our one dimensional experiments.

In our experiments, T = 5 and Ω = (0, 1)d for d ∈ {1, 2}. Further, our temporal

bases Ξq

n and Ξq−1n are nodal bases with q ∈ {1, 2, 4} and our spatial wavelet bases

Ψj are finite element wavelet bases with parameters p = ˜p∈ {1, 2} and γ = ˜γ = 32

constructed in [16]. Moreover, our spatial and temporal mesh are both uniform as follows. Let

M1:= [0, 1],

and

M21:= {(x1, x2) ∈ [0, 1]2: x2≥ x1} and M22:= {(x1, x2) ∈ [0, 1]2: x2≤ x1}.

For d = 1, our initial spatial mesh M0 consists of the two intervals resulted from a

dyadic refinement of M1, i.e., M0= {[0,12], [12,1]}. Analogously for d = 2, our initial

spatial mesh M0 consists of the eight triangles resulted from dyadic refinements of

M21 and M22. Thus, our spatial meshes Mj for j ≥ 0 consist of d2d(j+1) uniform

d-simplices for d ∈ {1, 2} (and so Vj = Vj(p) and Ψj can be defined, cf. Subsection

3.2). Finally, our temporal stepsizes hn are defined by either

hn= ¯h1(j) := T l T q2− p+1 q+1(j+1) m or hn= ¯h2(j) := T l T q2− p+1 2q (j+1) m ,

with ⌈·⌉ being the ceiling function (thus, we focus on L2 rates of convergence, cf.

Remark 2.2). Note that ¯h1(j) ≤ ¯h2(j). In both cases h = hn and N h = T .

Recall that the system matrices in (4.6) and (4.14) are Mj and Sj respectively,

and their preconditioners DMj and DSj are independent on the eigenvalues of the

matrix [ ˙A−1A]2 from Subsection 4.2. The condition numbers of the corresponding

preconditioned matrices (i.e., D−12

MjMjD −1 2 Mj and D −1 2 Sj SjD −1 2 Sj ), being dependent on p, will be denoted as κp,M,j and κp,S,j respectively.

In contrary, the system matrices in (4.10) and (4.15) depend on the eigenvalues of [ ˙A−1A]2, so do their preconditioners, as followed from Proposition 4.4. In particular,

the form of the preconditioners depend on whether the corresponding eigenvalues of [ ˙A−1A]2 are real or complex (note that two system matrices corresponding to a

pair of complex conjugated eigenvalues have the same preconditioner). For q = 1, [ ˙A−1A]2 has one (positive) real eigenvalue. For q = 2, it has one pair of complex

conjugated eigenvalues with positive real part. For q = 4 it has one pair of complex conjugated eigenvalues with positive real part and one pair of complex conjugated eigenvalues with negative real part. We did not numerically investigate the case q = 3, because [ ˙A−1A]2 has one (positive) real eigenvalue and one pair of complex

(19)

preconditioners are expected to behave similarly to the preconditioners in the two cases q ∈ {1, 2}. For q ∈ {1, 2}, the condition numbers of the preconditioned matrices, being dependent also on ¯hi(j) for i ∈ {1, 2}, will be denoted as

κ(p,q),i,j.

For q = 4, the condition numbers of the preconditioned matrices will be denoted as κ+(p,q),i,j and κ−(p,q),i,j,

corresponding to the pair of complex conjugated eigenvalues with positive and nega-tive real part respecnega-tively.

We approximated the condition numbers with the Lanczos method. In addition, we approximated the applications of S−1j for q ∈ {2, 4} with (preconditioned) CG as mentioned before. The numerical results are presented in Tables 5.1 and 5.2 for d = 1, and in Tables 5.3 and 5.4 for d = 2.

Table 5.1

κp,M,jand κp,S,j for d = 1 (with j = 10)

κp,M,10 κp,S,10

p= 1 2.2498e+00 1.5536e+01 p= 2 2.9034e+00 9.3342e+00

Table 5.2 κ(p,q),i,j, and κ+(p,q),i,jand κ

(p,q),i,j for d = 1 (with j = 10)

q= 1 q= 2 q= 4

κ(p,q),i,10 κ(p,q),i,10 κ+(p,q),i,10 κ−(p,q),i,10

p= 1 i= 1 5.1423e+00 1.4803e+01 1.4007e+01 1.5271e+01 i= 2 5.1423e+00 1.3523e+01 1.4771e+01 1.5491e+01 p= 2 i= 1 2.9033e+00 2.6361e+01 1.0945e+01 1.5670e+01 i= 2 2.9033e+00 1.8713e+01 7.9695e+00 8.6116e+00

Table 5.3

κp,M,jand κp,S,j for d = 2 (with j = 6)

κp,M,6 κp,S,6

p= 1 1.2214e+01 4.8632e+01 p= 2 1.2344e+01 5.9916e+01

Table 5.4 κ(p,q),i,j, and κ+(p,q),i,jand κ

(p,q),i,jfor d = 2 (with j = 6)

q= 1 q= 2 q= 4

κ(p,q),i,6 κ(p,q),i,6 κ+(p,q),i,6 κ−(p,q),i,6

p= 1 i= 1 2.0131e+01 3.7422e+01 4.9909e+01 5.5873e+01 i= 2 2.0131e+01 4.7867e+01 5.2218e+01 5.4894e+01 p= 2 i= 1 1.2414e+01 2.6170e+01 4.5416e+01 5.7078e+01 i= 2 1.2414e+01 3.9245e+01 5.5161e+01 6.0896e+01

(20)

From the numerical results, we draw the following conclusions:

• As discussed in [16], the preconditioned mass- and stiffness matrices above are, in a certain sense, quite close to the optimal ones. In comparison to the preconditioned mass- and stiffness matrices, our diagonal preconditioners show thus an overall satisfactory performance.

• The preconditioners correspond to pairs of complex conjugated eigenvalues with negative real part behave similarly to those correspond to pairs of com-plex conjugated eigenvalues with positive real part.

• For q = 1, the condition numbers of the preconditioned system matrices are proportional to those of the preconditioned mass matrices (which can also be theoretically verified).

• For q ∈ {2, 4}, the condition numbers of the preconditioned system matrices tend to be more proportional to those of the preconditioned stiffness matrices than to those of the preconditioned mass matrices.

REFERENCES

[1] M. Benzi and D. Bertaccini, Real-Valued Iterative Algorithms for Complex Symmetric Linear Systems, (2006). Submitted.

[2] S. Brenner and L. Scott, The Mathematical Theory of Finite Element Methods, Springer-Verlag, New York, 1994.

[3] P. Ciarlet, The Finite Element Method for Elliptic Problems, North-Holland Publishing Co., 1978.

[4] B. Cockburn, Discontinuous Galerkin methods, ZAMM Z. Angew. Math. Mech., 83 (2003), pp. 731–754.

[5] A. Cohen, Numerical analysis of wavelet methods, vol. 32 of Studies in Mathematics and its Applications, North-Holland Publishing Co., Amsterdam, 2003.

[6] W. Dahmen, Wavelet and multiscale methods for operator equations, in Acta numerica, 1997, vol. 6 of Acta Numer., Cambridge Univ. Press, Cambridge, 1997, pp. 55–228.

[7] W. Dahmen, A. Kunoth, and K. Urban, Biorthogonal spline wavelets on the interval – stability and moment conditions, Appl. Comput. Harmon. Anal., 6 (1999), pp. 132–196. [8] W. Dahmen and R. Stevenson, Element-by-element construction of wavelets satisfying

sta-bility and moment conditions, SIAM J. Numer. Anal., 37 (1999), pp. 319–352.

[9] R. Falk and G. Richter, Explicit Finite Element Methods for Symmetric Hyperbolic Equa-tions, SIAM J. Numer. Anal., 36 (1999), pp. 935–952.

[10] D. French and S. Jensen, Long-time behaviour of arbitrary order continuous time Galerkin schemes for some one-dimensional phase transition problems, IMA J. Numer. Math, 14 (1994), pp. 421–442.

[11] D. French and T. Peterson, A continuous space-time finite element method for the wave equation, Math. Comp., 65 (1996), pp. 491–506.

[12] G. Golub and C. Van Loan, Matrix Computations, The Johns Hopkins University Press, Baltimore, MD, third ed., 1996.

[13] C. Johnson, Discontinuous Galerkin finite element methods for second order hyperbolic prob-lems, Comput. Methods Appl. Mech. Eng., 107 (1993), pp. 117–129.

[14] P. Monk and G. Richter, A Discontinuous Galerkin Method for Linear Symmetric Hyperbolic Systems in Inhomogeneous Media, J. Sci. Comput., 22-23 (2005), pp. 443–477.

[15] M. Negreanu, A.-M. Matache, and C. Schwab, Wavelet filtering for exact controllability of the wave equation., SIAM J. Sci. Comput., 28 (2006), pp. 1851–1885.

[16] H. Nguyen and R. Stevenson, Finite element wavelets with improved quantitative properties, (2007). Submitted.

[17] M. Renardy and R. Rogers, An Introduction to Partial Differential Equations, Springer-Verlag, New York, 1993.

[18] R. Stevenson, On the compressibility of operators in wavelet coordinates, SIAM J. Math. Anal., 35 (2004), pp. 1110–1132 (electronic).

[19] T. von Petersdorff and C. Schwab, Numerical Solution of Parabolic Equations in High Dimensions, Mathematical Modelling and Numerical Analysis, 38 (2004), pp. 93–128.

References

Related documents

But instead of using simple differential equations the systems will be described stochastically, using the linear noise approximation of the master equation, in order to get a

Establishing a language for innovation to institutionalise a definition, use corporate culture as a control system, standardising internal processes, fostering a learning

xed hypotheses in general linear multivariate models with one or more Gaussian covariates. The method uses a reduced covariance matrix, scaled hypothesis sum of squares, and

related to tough-mindedness in both samples, and is mostly unrelated to uncertainty avoidance (except for regression analysis of the Swedish data, where conservation is related

For some constructions, we shall need nonsymmetric tree monomials, that is planar trees with decorations of internal vertices but no labelling of leaves (or, when it is needed

Hence given a labeled map and an LTL-formula entered via the skeleton, the framework generates respective cost maps that can be used to steer a path planning algorithm.. However

time integration using a time step ∆t &gt; 0 for the ve tor x an be made by any of the standard rst order

We construct a multiscale scheme based on the heterogeneous multiscale method, which can compute the correct coarse behavior of wave pulses traveling in the medium, at a