Iain Smears

Robust and efficient preconditioners for the discontinuous Galerkin time-stepping method

Keywords: discontinuous Galerkin, time discretizations, parabolic PDE, preconditioning, conjugate gradient algorithm

The discontinuous Galerkin time-stepping method has many advantageous properties for solving parabolic equations. However, it requires the solution of a large nonsymmetric system at each timestep. This work develops a fully robust and efficient preconditioning strategy for solving these systems. Drawing on parabolic inf-sup theory, we first construct a left preconditioner that transforms the linear system to a symmetric positive definite problem to be solved by the preconditioned conjugate gradient algorithm. We then prove that the transformed system can be further preconditioned by an ideal block diagonal preconditioner, leading to a condition number bounded by 4 for any time-step size, any approximation order and any positive-definite self-adjoint spatial operators. Numerical experiments demonstrate the low condition numbers and fast convergence of the algorithm for both ideal and approximate preconditioners, and show the feasibility of the high-order solution of large problems.

Introduction

The discontinuous Galerkin (DG) time-stepping method is a single-step implicit scheme defined by a variational temporal discretization of parabolic evolution equations that generalizes the backward Euler method to higher-order approximations [START_REF] Delfour | Discontinuous Galerkin methods for ordinary differential equations[END_REF][START_REF] Eriksson | Time discretization of parabolic problems by the discontinuous Galerkin method[END_REF][START_REF] Jamet | Galerkin-type approximations which are discontinuous in time for parabolic equations in a variable domain[END_REF][START_REF] Hulme | Discrete Galerkin and related one-step methods for ordinary differential equations[END_REF]. For an introduction to this time discretization scheme in the context of abstract parabolic problems, we refer the reader to (Thomée, 2006). In certain cases, it coincides with the Radau IIA Implicit Runge-Kutta (IRK) schemes and the subdiagonal Padé approximations to the exponential function [START_REF] Axelsson | A class of A-stable methods[END_REF][START_REF] Hairer | Solving ordinary differential equations[END_REF][START_REF] Makridakis | A posteriori error analysis for higher order dissipative methods for evolution problems[END_REF]. It can be coupled with standard spatial discretization schemes, such as finite difference, finite element or spectral methods; in particular, when coupled to a spatial finite element method (FEM), it leads to a tensor-product space-time FEM.

The DG time-stepping method features many advantages that make it an attractive choice for solving parabolic problems. First, it permits arbitrarily high-order approximation to the solution, along with superconvergence at the time-step nodes [START_REF] Akrivis | Galerkin time-stepping methods for nonlinear parabolic equations[END_REF][START_REF] Chrysafinos | Error estimates for the discontinuous Galerkin methods for parabolic equations[END_REF]Makridakis & Babuška, 1997), see also (Schötzau & Schwab, 2000) for optimal-order a priori error estimates in natural norms with explicit dependence on the polynomial degree. Unlike linear multistep schemes, it is thus not constrained by the Dahlquist barrier theorem and it does not require an auxiliary scheme to compute the first few solution values. Furthermore, the DG time-stepping method allows fully variable time-step sizes, approximation orders, and even spatial mesh refinement/coarsening between time-steps; it is thus well-suited for adaptive algorithms driven by rigorous a posteriori error es-timates [START_REF] Akrivis | Optimal order a posteriori error estimates for a class of Runge-Kutta and Galerkin methods[END_REF][START_REF] Eriksson | Adaptive finite element methods for parabolic problems. I. A linear model problem[END_REF][START_REF] Eriksson | Adaptive finite element methods for parabolic problems. II. Optimal error estimates in L ∞ L 2 and L ∞ L ∞[END_REF][START_REF] Makridakis | A posteriori error analysis for higher order dissipative methods for evolution problems[END_REF]Schötzau & Wihler, 2010). The DG time-stepping method also permits the temporal version of hp-refinement, where one varies both the time-step size as well as the approximation order on each time-step, thereby yielding exponential convergence rates for a broad class of solutions with singularities induced by the initial datum or by the source term (Schötzau & Schwab, 2000). In these applications, it is common to use high temporal polynomial degrees in order to match the high-order spatial approximation, see for instance the experiments in (Schötzau & Schwab, 2000;[START_REF] Von Petersdorff | Numerical solution of parabolic equations in high dimensions[END_REF]Schötzau & Wihler, 2010;[START_REF] Werder | hp-discontinuous Galerkin time stepping for parabolic problems[END_REF].

In applying the DG time-stepping method in practice, we are faced with the challenge of solving a large nonsymmetric linear system at each time-step. In this work, we are interested in developing preconditioned iterative methods for solving these linear systems. We focus here on the DG timestepping method in the context of an abstract semi-discrete evolution problem of the form M U (t) + AU(t) = f (t) in (0, T ),

(1.1)

where the solution U : [0, T ] → V, with V a finite dimensional space, where M and A are symmetric positive definite matrices. Semi-discrete problems of the form (1.1) arise from the application of a wide range of spatial discretizations, including conforming and nonconfirming finite elements, finite differences, and spectral methods, to a broad class of parabolic problems, such as general self-adjoint second-order and fourth-order parabolic partial differential equations (PDE). We emphasize that our results are valid for general symmetric positive definite matrices M and A, although the reader may consider the heat equation as a concrete example, with M and A respectively representing the mass and stiffness matrices for a suitable approximation space V. The DG time-stepping method applied to the evolution problem (1.1) leads to a sequence of linear systems of the general block form where τ denotes the time-step size, the polynomial degree p defines the approximation order of the scheme, the solution coefficients u k ∈ V for k = 0, . . . , p, and, after mapping the time-step interval to the reference interval (-1, 1), we have

b jk := 1 -1 φ k φ j ds + φ k (-1)φ j (-1), c jk := 1 2 1 -1 φ k φ j ds, (1.3)
where {φ k } p k=0 is a chosen basis of P p the space of real-valued polynomials of degree at most p. For the case p = 0, the DG method reduces to the backward Euler method and the system is therefore symmetric. However, for p 1, the system matrix B of (1.2) is nonsymmetric with dimension dim V×(p+1), which is considerably larger than for linear multistep methods; even for moderate sizes of A and M, standard direct solution algorithms can be prohibitively expensive. Unlike the block triangular linear systems obtained from Diagonally IRK (DIRK) and Singly Diagonally IRK (SDIRK) schemes (see [START_REF] Hairer | Solving ordinary differential equations[END_REF] and the references therein), the system (1.2) does not immediately reveal any simple structure offering a straightforward solution algorithm. See also (Hairer & Wanner, 2010, p. 129-130) for a discussion of SIRK methods.

Since the DG time-stepping method is connected to the Radau IIA IRK scheme, it is interesting to relate our approach to the literature on solving the systems of IRK schemes. In particular, one of the earliest approaches [START_REF] Butcher | On the implementation of implicit Runge-Kutta methods[END_REF][START_REF] Bickart | An efficient solution process for implicit Runge-Kutta methods[END_REF] for solving the linear systems of general IRK schemes, such as (1.2), is based on transforming the system matrix to a block-diagonal matrix with blocks of the form M + τ µ i A, where the {µ i } p i=0 denote here the generalized eigenvalues of (b jk ) with respect to (c jk ), and where the transformation is given by the corresponding eigenvectors. This leads to p + 1 independent smaller linear systems that can be solved in parallel. It turns out that for the DG time-stepping method, the generalized eigenvalues µ i are complex numbers related to the roots of the denominator of a rational Padé approximation to the exponential function. Therefore, the resulting transformed system is complex-valued and non-Hermitian despite M and A being symmetric. In addition to the increased computational cost of complex arithmetic, this approach has an important shortcoming in terms of robustness and numerical stability for high-order approximations, due to ill-conditioning of the eigenvector transformations, as pointed out in (von Petersdorff & Schwab, 2004, Remark 5.4). This issue can be avoided by employing alternative factorizations, at the expense of the block-diagonal structure of the transformed problem: Schötzau & Schwab (2000) propose a factorization based on the Schur decomposition theorem, leading to a block-triangular complex transformed problem; the solution is then obtained by solving p + 1 complex non-Hermitian systems in sequence.

Preconditioned iterative methods offer an alternative approach to decoupling the system by complex transformations. In this direction, [START_REF] Richter | Efficient numerical realization of discontinuous Galerkin methods for temporal discretization of parabolic problems[END_REF] propose a linear iterative fixed point scheme based on an approximation of the block LU factorization of B, and they analyse the contraction rates of their method for p 3. An alternative approach is to apply directly a preconditioned Krylov subspace method to (1.2), as proposed in [START_REF] Mardal | Order-optimal preconditioners for implicit Runge-Kutta schemes applied to parabolic PDEs[END_REF] propose a block-diagonal preconditioner for IRK schemes to be used with GMRES. They show that the preconditioner is robust with respect to the time-step size τ, but their experiments show that it is not robust with respect to p. [START_REF] Weller | Efficient preconditioning of variational time discretization methods for parabolic partial differential equations[END_REF] develop a preconditioner specifically for the DG time-stepping method for p = 1 (as well as for the related continuous Galerkin method of same size) based on an approximate Schur complement preconditioner for one of the unknown coefficients u j in (1.2). It is thus apparent from these references that finding preconditioning strategies for the DG time-stepping method that are robust with respect to the polynomial degree p has been a challenging open problem.

In this work, we propose a robust and efficient preconditioned iterative method for the DG timestepping method. Instead of focusing exclusively on the block structure of (1.2), our approach draws upon the inf-sup analysis of the method and the underlying continuous analysis of parabolic PDEs, and we take advantage of the variational structure of the DG time-stepping method in an essential way. First, in section 3.1, we construct and apply a left-preconditioner P to the linear system Bu = f given by (1.2), resulting in a preconditioned system Lu = g with the key benefit that the transformed matrix L := P B is symmetric positive definite. We immediately point out that P = B , i.e. we are not forming the normal equations of the system. Instead, our construction of P is motivated by parabolic inf-sup theory, and we show that the matrix L represents the discrete version of the natural parabolic energy norm of the underlying evolution equation: for example, in the context of second-order parabolic PDEs, the matrix L is a discrete Gram matrix for the natural solution space H 1 (H -1 ) ∩ L 2 (H 1 ); we refer the reader to [START_REF] Wloka | Partial differential equations[END_REF] for an introduction to the continuous analysis of parabolic problems. The transformed symmetric positive definite system can therefore be solved by the preconditioned conjugate gradient (PCG) algorithm [START_REF] Hestenes | Methods of conjugate gradients for solving linear systems[END_REF]Málek & Strakoš, 2015;[START_REF] Wathen | Preconditioning[END_REF], which, in our case, minimizes the error in the physically relevant norm over a Krylov subspace. In order to obtain the fast and robust convergence of the PCG algorithm, in section 3.2, we construct a spectrally equivalent preconditioner H for L, such that the condition number of the preconditioned system satisfies

κ(H -1 L) 4, (1.4)
independently of all parameters τ, p, M and A. Therefore, the preconditioner is fully robust with respect to all problem and discretization parameters, including the polynomial degree. Furthermore, the preconditioners are efficient, firstly in the sense of guaranteeing the fast convergence of the PCG algorithm in the physically relevant norm, and secondly in the sense of computational cost, for the following reasons.

In section 4 we show that the preconditioners are well-suited for parallelization over the blocks and involve only simpler matrices of the form A and M + µA with real positive µ for which efficient solvers are often available. Furthermore, we show experimentally in section 5 that the ideal preconditioners can be approximated in practice by cheap spectrally equivalent approximations, such as a small number of multigrid V -cycles. We refer the reader to (Wathen, 2015, p. 367) and (Hiptmair, 2006, p. 705) for further discussions of the notion of efficiency of preconditioners. This paper is organized as follows: after introducing in detail the DG time-stepping method in section 2, we present the preconditioning strategy in section 3, where we construct the preconditioners P and H, and where we establish the condition number bound (1.4). Section 4 considers the efficient implementation of the method, and section 5 presents the results of numerical experiments testing the robustness and efficiency of the preconditioners.

Preliminaries

In this section, we introduce in detail the DG time-stepping method in the context of self-adjoint semidiscrete dissipative evolution equations. We also introduce two key ingredients in our approach. The first ingredient is the well-known temporal reconstruction operator commonly used in a posteriori analysis [START_REF] Makridakis | A posteriori error analysis for higher order dissipative methods for evolution problems[END_REF], which is the subject of section 2.3. The second ingredient is a spectral equivalence result for preconditioners from [START_REF] Pearson | A new approximation of the Schur complement in preconditioners for PDE-constrained optimization[END_REF], which we present in section 2.4.

Approximation space

Let V denote a finite dimensional real vector space, equipped with a given basis {e i } dim V i=1 . Let V be equipped with two inner products (•, •) M and (•, •) A , and let M and A be their matrix representations, given by M i j := (e i , e j ) M and A i j := (e i , e j ) A for all i, j = 1, . . . , dim V. The inner products (•, •) M and (•, •) A induce the norms • M and • A on V. Let P p denote the space of real-valued polynomials of degree at most p, and let V p be the space of V-valued polynomials of a single real variable with degree at most p. For example, if {φ j } p j=0 is a basis of P p , then every v ∈ V p is of the form

v : s → v(s) = p ∑ j=0 v j φ j (s),
with coefficients v j ∈ V for each j = 0, . . . , p. We gather these coefficients in the vector v = (v 0 , . . . , v p ) ∈ V p+1 . It follows that dim V p , the dimension of the space V p , is equal to (p + 1) × dim V. An equivalent point of view is to consider V p as the tensor-product space derived from P p and V.

REMARK 2.1 In this setting, it is natural to view functions in V p as mappings from time into V, and as a slight abuse of standard terminology, we will say that {φ j } p j=0 forms a basis of V p . Of course, this must be interpreted as the lengthier statement that {e i φ j } j=0,...,p i=1,...,dim V forms a basis of V p , where

{e i } dim V i=1 is a basis of V.
If L is a linear operator between V and either itself or its dual V * , then we can extend L to V p by applying it coefficient-wise:

Lv(s) := p ∑ j=0 (Lv j ) φ j (s) ∀ v = p ∑ j=0 v j φ j ∈ V p .
(2.1)

The inner products (•, •) M and (•, •) A also extend to V p in the natural way. Likewise, if π : P p → P p is a linear operator then we define πv := ∑ p j=0 v j (πφ j ). For instance, we define the time derivative

v of v ∈ V p by v := ∑ p j=0 v j φ j . Let (•, •) L 2 and • L 2 denote
respectively the L 2 -inner product and L 2 -norm over the interval (-1, 1). Let {L k } k 0 denote the set of Legendre polynomials, as defined for instance in (Gradshteyn & Ryzhik, 2015, Sec. 8.9). The Legendre polynomials are orthogonal in the L 2 -inner product: for all k, j 0,

(L k , L j ) L 2 = 2 2k + 1 δ k j . (2.2) Furthermore, L k (1) = 1 and L k (-1) = (-1) k for all k 0.

The DG time-stepping method

After mapping a given current time-step interval to the reference interval (-1, 1), the DG time-stepping method leads to the discrete problem of finding u ∈ V p such that

B(u, v) = F (v) ∀ v ∈ V p , (2.3)
where F is a bounded linear functional on V p , and the bilinear form B :

V p × V p → R is defined by B(u, v) := 1 -1 u , v M ds + (u(-1), v(-1)) M + τ 2 1 -1 (u, v) A ds, (2.4)
where τ denotes the time-step size and p denotes the polynomial degree of the approximation. In practice, the time-step size, the polynomial degree, or even the matrices M and A may vary between time-steps, although on any given time-step the linear system has the general form of (2.3). Given a basis {φ j } p j=0 for V p , the problem (2.3) can be represented by a linear system

Bu = f, (2.5)
where B is a (p+1)×(p+1) block matrix, where u ∈ V p+1 is the vector of coefficients of the expansion of u, and where

f = (f 0 , • • • , f p ) with each f j ∈ V * , j = 0, . . . , p, being the restriction of F to span φ j ⊗ V.
Therefore, the system (2.5) has the block structure shown in (1.2).

REMARK 2.2 In order to motivate our approach to preconditioning, it will be helpful to bear in mind the following point concerning the structure of the bilinear form B. It is well-known that the bilinear form B enjoys the following coercivity property

B(v, v) = 1 2 v(1) 2 M + 1 2 v(-1) 2 M + τ 2 1 -1 v 2 A ds ∀ v ∈ V p , (2.6)
which enables us to deduce the well-posedness of (2.3). Unfortunately, the norm defined by the righthand side of (2.6) does not include the time derivative, and thus coercivity and boundedness of B in this norm can only be obtained from an inverse inequality for the finite dimensional space V p , at the expense of introducing constants that are not robust with respect to τ and p. This point suggests that norm preconditioners for B that are based on the right-hand side of (2.6) are unlikely to be robust with respect to the discretization parameters.

As we shall see below, our approach is based on the inf-sup stability of the bilinear form B, which provides a sharper analysis of the structure of the problem than the coercivity result of (2.6). Some of our main tools are the reconstruction operator defined in section 2.3, and a corresponding suitable negative norm along with a key spectral equivalence result, which we recall in section 2.4.

Reconstruction operator

We introduce the reconstruction operator I : P p → P p+1 , defined by

I v := v -v(-1) (-1) p (L p -L p+1 ) 2 ∀ v ∈ P p .
(2.7)

As noted above, I naturally extends to an operator from V p to V p+1 . As explained in Remark 2.3 below, I is the reconstruction operator commonly used in a posteriori error analysis [START_REF] Makridakis | A posteriori error analysis for higher order dissipative methods for evolution problems[END_REF]. The key benefit of the reconstruction operator I for our purposes is that we may express the bilinear form B in the following equivalent form

B(u, v) = 1 -1 (I u) , v M + τ 2 (u, v) A ds ∀ u, v ∈ V p . (2.8)
We emphasize here that (I u) ∈ V p since I u ∈ V p+1 for u ∈ V p . To see how (2.8) is obtained from (2.7), we note that the properties of the Legendre polynomials imply that, for any v ∈ V p ,

I v(1) = v(1), I v(-1) = 0, 1 -1 (I v, w) M ds = 1 -1 (v, w) M ds ∀ w ∈ V p-1 .
(2.9) Substituting w for w , which belongs to V p-1 whenever w ∈ V p , in (2.9) and using integration by parts shows that 1 -1

(I v) , w M ds = 1 -1 (v , w) M ds + (v(-1), w(-1)) M ∀ w ∈ V p .
(2.10)

The equivalent form of B given in (2.8) then follows from (2.10).

REMARK 2.3 The operator I is the reconstruction operator commonly used in the a posteriori error analysis of the DG time-stepping method [START_REF] Makridakis | A posteriori error analysis for higher order dissipative methods for evolution problems[END_REF]. This operator is often defined by the interpolation conditions I v(s k ) = v(s k ) for all k = 1, . . . , p + 1, where -1 < s k s p+1 = 1 are the Gauss-Radau quadrature points, in addition to a further condition, chosen here as I v(-1) = 0. It turns out that the definition given above in (2.7) and this interpolatory definition are equivalent, since the Gauss-Radau points are the roots of the polynomial (1 -s)P

(1,0) p = L p -L p+1 , see [START_REF] Gautschi | Numerical analysis[END_REF]) and (Gradshteyn & Ryzhik, 2015, eq. 8.961.5). We note here that a straightforward consequence of the above properties is that, for any v ∈ P p , we have v ≡ 0 in (-1, 1) if and only if (I v) ≡ 0 in (-1, 1). We also note that we will not need the Gauss-Radau points for the implementation of the preconditioners in this work.

Negative norms and a result of Pearson and Wathen

In addition to the norms • M and • A , we will also use the negative norm

• MA -1 M defined by v MA -1 M := sup w∈V\{0} (v, w) M w A ∀ v ∈ V.
(2.11)

The negative norm • MA -1 M can be equivalently characterized by the identity

v 2 MA -1 M = v MA -1 Mv ∀ v ∈ V, (2.12)
where v is identified with its vector representation in the basis {e i } dim V i=1 . Indeed, (2.12) follows from the upper bound v 2

MA -1 M v MA -1
Mv, which is obtained by choosing w = A -1 Mv in (2.11), and from the lower bound v 2

MA -1 M v MA -1
Mv, which is obtained by applying the Cauchy-Schwarz inequality as follows:

w Mv = w A A -1 Mv √ w A w √ v MA -1
Mv, where we have simplified

(A -1 Mv) A(A -1 Mv) = v MA -1 Mv. The identity (2.12) thus shows that the norm • MA -1 M is in fact induced by an inner product (•, •) MA -1 M represented by the matrix MA -1 M.
The following result due to J. W. Pearson and A. J. Wathen [START_REF] Pearson | A new approximation of the Schur complement in preconditioners for PDE-constrained optimization[END_REF] will play a key part in the construction of our preconditioners.

LEMMA 2.1 Let A and M be arbitrary symmetric positive definite matrices and let µ 0 be a nonnegative real number. Then we have

1 2 v MA -1 M + µ A v v M + √ µ A A -1 M + √ µ A v 1 ∀ v ∈ V \ {0}.
(2.13)

Proof. For the original proof of this result, see (Pearson & Wathen, 2012, Thm 4). We provide here an alternative proof which highlights the negative norm structure of the matrix MA -1 M, as given in (2.11). First, the upper bound of (2.13) is immediate from

v 2 MA -1 M + µ v 2 A v 2 MA -1 M + 2 √ µ v 2 M + µ v 2 A . Now, for any v ∈ V, we have v 2 M v MA -1 M v A by (2.11). It follows from the inequality ab 1 2 a 2 + 1 2 b 2 for any a, b ∈ R that v (M + √ µ A) A -1 (M + √ µ A) v = v 2 MA -1 M + 2 √ µ v 2 M + µ v 2 A 2 v 2 MA -1 M + 2µ v 2 A = 2 v MA -1 M + µ A v.
(2.14)

Note that (2.14) is precisely the lower bound of (2.13).

Preconditioners

Left preconditioner

In this section, we construct a left preconditioner that will transform the linear system (2.5) to a symmetric positive definite system that represents a discrete parabolic energy norm. Our left preconditioner is defined simply in terms of a substitution for the test function v appearing in the bilinear form B(u, v).

We start by defining the operator P :

V p → V p by Pv := A -1 M(I v) + τ 2 v, (3.1)
where we recall the definition of the reconstruction operator I from (2.7) and its natural to extension to V p as explained in section 2.1. The fact that Pv ∈ V p for any v ∈ V p implies that the solution u of (2.3) also solves

L (u, v) = G (v) ∀ v ∈ V p , (3.2)
where the bilinear form L and linear functional G are obtained by substituting the test function Pv in place of v:

L (u, v) := B(u, Pv), G (v) := F (Pv) ∀ u, v ∈ V p . (3.3)
We note that the operator P can be viewed as defining a left preconditioner for the linear system (2.5) that represents (2.3). Indeed, let P denote the matrix representation of the linear operator P in the basis {φ j } p j=0 . Then, for any functions u and v ∈ V p , and their respective vector representations u and v ∈ V p+1 , we have B(u, Pv) = v P Bu and F (Pv) = v P f. Therefore, the linear system (3.2) is equivalent to Lu = g, (3.4)

with L = P B and g = P f denoting respectively the left-preconditioned matrix and right-hand side.

As the following theorem shows, L represents a discrete version of the natural energy norm for parabolic problems. Indeed, in applications to second-order parabolic PDEs, L is comparable to the inner product of L 2 (H 1 )∩H 1 (H -1 ), which is the natural solution space of the continuous problem [START_REF] Wloka | Partial differential equations[END_REF]. THEOREM 3.1 Let the bilinear form L : V p × V p → R and linear functional G : V p → R be defined by (3.3). Then, for any functions u and v ∈ V p , we have the identity

L (u, v) = 1 -1 (I u) , (I v) MA -1 M ds + τ 2 4 1 -1 (u, v) A ds + τ 2 (u(1), v(1)) M + τ 2 (u(-1), v(-1)) M (3.5)
Therefore, L is symmetric and positive definite on V p , and for any symmetric positive definite matrices A and M, any p 0 and any τ > 0, we have

L (v, v) v 2 D , |L (v, w)| 2 v D w D ∀ v, w ∈ V p , (3.6)
where

• D := D(•, •)
is the norm induced by the auxiliary bilinear form D defined by:

D(u, v) := 1 -1 (I u) , (I v) ) MA -1 M ds + τ 2 4 1 -1 (u, v) A ds. (3.7)
Thus, the function u ∈ V p solves (2.3) if and only if u solves (3.2).

Proof. First, it is straightforward to obtain the following identities which hold for any u, v ∈ V p :

(I u) , A -1 M(I v) M = ((I u) , (I v) ) MA -1 M , u, A -1 M(I v) A = (u, (I v) ) M . (3.8)
Next, we use the identities in (3.8) to simplify the different terms in L (u, v) = B(u, Pv), and eventually we obtain

L (u, v) = 1 -1 (I u) , (I v) ) MA -1 M ds + τ 2 4 1 -1 (u, v) A ds + τ 2 1 -1 (I u) , v M + u, (I v) M ds.
(3.9) By (2.10), we have 1 -1

(I u) , v M + u, (I v) M ds = 1 -1 d ds (u, v) M ds + 2 (u(-1), v(-1)) M = (u(1), v(1)) M + (u(-1), v(-1)) M ,
which, after substitution of the last terms in (3.9), implies the equivalent form for L given in (3.5).

Next, to show (3.6), we note that first that the lower bound L (v, v) v 2 D is immediate, whereas the upper bound |L (v, w)| 2 v D w D follows from the application of the Cauchy-Schwarz inequality to (3.9). Finally, it follows that the problem (3.2) has a unique solution, momentarily denoted ũ ∈ V p . Moreover, it is well-known that (2.3) has a unique solution u ∈ V p . Since the definition of L and G in (3.3) shows that u is a solution of (3.2), we deduce from uniqueness that ũ = u. Therefore, the problems (2.3) and (3.2) are equivalent. REMARK 3.1 Theorem 3.1 can be viewed as part of the inf-sup analysis of the DG time-stepping method: defining the energy norm

• L := L (•, •) and the norm • X := 1 -1 • 2
A ds, we have

u L = sup v∈V p \{0} B(u, v) v X , |B(u, v)| u L v X ∀ u, v ∈ V p , (3.10)
where the first equality is attained by choosing the optimal test function v = Pu, since Pu X = u L . This observation highlights the fact that the operator P represents the operator that gives the optimal test function in the inf-sup analysis of the DG time-stepping method. It is in this sense that our preconditioning strategy is directly motivated by the inf-sup theory of the DG time-stepping method.

Spectrally equivalent norm preconditioner

As shown by Theorem 3.1, the bilinear form L is symmetric and positive definite. Therefore, the linear system (3.4) can be solved iteratively by the preconditioned conjugate gradient (PCG) algorithm [START_REF] Hestenes | Methods of conjugate gradients for solving linear systems[END_REF]Málek & Strakoš, 2015;[START_REF] Wathen | Preconditioning[END_REF]. In this section, we construct a spectrally equivalent and easily applicable preconditioner for L, thereby leading to the robust and fast convergence of the PCG algorithm. Recall that the bilinear form L is spectrally equivalent to the bilinear form D defined by (3.7). The first step in our construction of a preconditioner is to choose an advantageous temporal basis for V p that block-diagonalizes the matrix D that represents D. In a second step, we construct an easily applicable preconditioner, denoted by H, under this basis, and we apply Lemma 2.1 to show robust spectral equivalence between H and L.

3.2.1

Definition of the basis. It follows from Remark 2.3 that the symmetric bilinear form (u, v) → 1 -1 (I u) (I v) ds is positive definite and thus defines an inner-product on P p . Therefore, there exists a set of linearly independent polynomial eigenfunctions {ϕ j } p j=0 ⊂ P p and corresponding real positive eigenvalues {λ j } p j=0 ⊂ R >0 such that

λ j 1 -1 (I ϕ j ) (I v) ds = 1 -1 ϕ j v ds ∀ v ∈ P p , j = 0, . . . , p. (3.11)
We note that the eigenvalues and eigenfunctions generally depend on p, and are not necessarily hierarchical. The eigenfunctions are chosen to be orthonormalized:

1 -1 (I ϕ j ) (I ϕ k ) ds = δ jk ∀ 0 j, k p. (3.12)
We provide a practical method for computing this basis along with the corresponding eigenvalues in section 4.1. The following result characterizes the distribution of the eigenvalues.

THEOREM 3.2 Let p 0 be a nonnegative integer, and let λ 0 • • • λ p be the eigenvalues of the generalized eigenvalue problem (3.11). Then, there exists a positive constant C, independent of j and p, such that λ j C( j + 1) -2 ∀ j = 0, . . . , p, ∀ p 0.

(3.13)

Moreover, there exists a constant c, independent of p, such that

λ 0 • • • λ p c (p + 1) -4 . (3.14)
We leave the proof of Theorem 3.2 to Appendix A since it is based on some results presented in the later section 4.1. Figure 2 shows that the orders of the bounds of Theorem 3.2 are sharp. The eigenfunctions for p = 4 are shown in Figure 1.

Construction of the norm preconditioner.

It is advantageous to use the basis {ϕ j } p j=0 in computations since in this basis, the dominant terms in the bilinear form L , namely those belonging to the bilinear form D, are represented by a block-diagonal matrix. Indeed, let D denote the matrix representation of the bilinear form D, defined in (3.7), under the basis {ϕ j } p j=0 of V p , where we recall the terminology from Remark 2.1. Then, in this basis, the matrix D has a simple block-diagonal structure:

D = diag MA -1 M + τ 2 λ j 4 A p j=0
.

(3.15)

In other words, D is block-diagonal, with blocks MA -1 M + λ j τ 2 4 A for j = 0, . . . , p. Keeping in mind Lemma 2.1, it is then natural to define the norm preconditioner H by

H := diag M + τ λ j 2 A A -1 M + τ λ j 2 A p j=0 . (3.16)
The matrix H is symmetric positive definite, and its inverse is trivially given by

H -1 = diag    M + τ λ j 2 A -1 A M + τ λ j 2 A -1    p j=0
.

We propose to use H as a preconditioner for the PCG algorithm applied to (3.4). Each PCG iteration requires the application of H -1 , which comprises two applications of a solver for a weighted backward Euler step and one multiplication by A per block.

We now give the central result of this work, which shows that L and H are spectrally equivalent with fully robust bounds. THEOREM 3.3 Let p 0, and let {ϕ j } p j=0 ⊂ P p and {λ j } p j=0 ⊂ R >0 be the eigenfunctions and eigenvalues of (3.11). Let L be the matrix representation of L in the basis {ϕ j } p j=0 . Then, for any symmetric positive definite matrices A and M, any τ > 0 and any p 0, we have

1 2 v L v v H v 2 ∀ v ∈ V p+1 \ {0}.
(3.17)

Proof. Theorem 3.1, in particular (3.6), shows that

1 v L v v D v 2 ∀ v ∈ V p+1 \ {0}. (3.18)
Lemma 2.1, applied block-wise with µ = λ j τ 2 /4, implies that

1 2 v D v v H v 1 ∀ v ∈ V p+1 \ {0}. (3.19)
Therefore, (3.18) and (3.19) imply (3.17). Theorem 3.3 immediately implies that the condition number κ of the preconditioned system, defined as the ratio of extremal eigenvalues of H -1 L, satisfies κ 4.

(3.20)

Therefore, we may expect very fast convergence from the PCG algorithm for Lu = g using the preconditioner H. Indeed, for an initial guess u 0 ∈ V p , the iterates {u k } k 0 of the PCG algorithm satisfy the well-known bound [START_REF] Wathen | Preconditioning[END_REF] u

-u k L u -u 0 L 2 √ κ -1 √ κ + 1 k 2 3 k ∀ k 0, (3.21)
where we recall the energy norm • L = L (•, •). The fact that the condition number κ 4 shows that the left preconditioner P and the norm preconditioner H are very effective at preconditioning the original system matrix B, despite B being nonsymmetric and poorly conditioned. This is a key aspect of the efficiency of the proposed preconditioners. Furthermore, we highlight that (3.20) and (3.21) are valid for any time-step size τ, any polynomial degree p, and any symmetric positive definite matrices M and A. Thus the preconditioners are fully robust with respect to all discretization and problem parameters.

REMARK 3.2 The fact that the energy norm • L appears in the bound (3.21) is advantageous in practice. Since the endpoint value of the solution at s = 1 serves as initial datum for the next time-step, it is beneficial to be able to control the accuracy to which it is computed. For example, (3.5) shows that the energy norm controls the value at s = 1 in the M-norm, which is the natural norm for this quantity of interest. We note that the factor of τ appearing alongside v(1) 2 M in the energy norm can be scaled out as it appears in both the numerator and denominator of (3.21). Therefore (3.21) is a guaranteed convergence rate for this quantity of interest in the physically relevant norm.

Implementation

Computation of eigenfunctions and eigenvalues

The algorithm requires the computation of the eigenfunctions {ϕ j } p j=0 and eigenvalues {λ j } p j=0 defined in (3.11), which involves solving a symmetric positive definite eigenvalue problem of dimension p + 1.

-1 -0.5 0 0.5 1 -0.5

0 0.5 1 ϕ 0 ϕ 1 ϕ 2 ϕ 3 ϕ 4 FIG. 1. Eigenfunctions {ϕ j } p
j=0 defined by (3.11) and (3.12), computed for p = 4 by the method of section 4.1 and ordered by decreasing eigenvalue λ j λ j+1 . Note that the j-th eigenfunction ϕ j need not be of degree at most j. Furthermore, the set of eigenfunctions generally depends on p, although they are independent of A, M and τ. This poses little difficulty, as in practice p is usually small, especially in comparison to dim V. Since this step only depends on p and is fully independent of A, M and τ, it can be pre-computed to very high accuracy. We now show how to assemble this eigenvalue problem in a form where it can be solved numerically. Since the cases p = 0 and p = 1 can be easily computed by hand, we present a general method for p 2. In the following, we will use the following identity for the Legendre polynomials {L k } p k=0 , see (Gradshteyn & Ryzhik, 2015, Sec. 8.914):

L k = L k+1 -L k-1 2k + 1 ∀ k 0. (4.1)
LEMMA 4.1 Let p 2 be a nonnegative integer. Define the polynomials

ψ 0 := 1 √ 2 (L 1 + L 0 ) ψ p := L p -L p-1 √ 4p + 2 , ψ k := L k+1 -L k-1 √ 4k + 2 , k = 1, . . . , p -1. (4.2)
Then, there holds

(I ψ k ) = k + 1 2 L k ∀ k = 0, . . . , p. (4.3)
Therefore, we have, for each 0 k, j p,

1 -1 (I ψ k ) (I ψ j ) ds = δ k j . (4.4)
Proof. Consider the case k = 0: ψ 0 (-1) = 0, and L 1 = L 0 . Therefore I ψ 0 = ψ 0 , and thus (I ψ 0 ) = L 0 / √ 2, which is (4.3) for k = 0. Now consider the general case k = 1, . . . , p -1. Since L j (-1) = (-1) j for all j 0, we have ψ k (-1) = 0 and thus I ψ k = ψ k . Now, the identity (4.1) implies that

(I ψ k ) = 2k + 1 √ 4k + 2 L k = k + 1 2 L k , (4.5)
thus verifying (4.3) for k = 1, . . . , p -1. For the case k = p, we use the fact that L p (-1) -L p-1 (-1) = 2(-1) p to compute

I (L p -L p-1 ) = L p -L p-1 - 2(-1) p (-1) p 2 (L p -L p+1 ) = L p+1 -L p-1 . (4.6)
Therefore, similarly to (4.5), the identity (4.3) for k = p follows from (4.1). Define the matrix Z ∈ R (p+1)×(p+1) such that ψ k = ∑ p j=0 Z k j L j as in (4.2). Define the diagonal matrix D := diag{2/(2 j + 1)} p j=0 , which represents the L 2 -inner product in the Legendre polynomial basis. Define the matrix T ∈ R (p+1)×(p+1) by

T k j := 1 -1 ψ k ψ j ds, k, j = 0, . . . , p.
(4.7)

Note that T = ZDZ . It follows from (4.4) that the eigenvalue problem (3.11) can be expressed as: find V ∈ R (p+1)×(p+1) and {λ j } p j=0 such that

TV = V diag{λ j } p j=0 , VV = V V = Id. (4.8)
The eigenvalue problem (4.8) can thus be solved numerically by standard eigenvalue solvers. The eigenfunctions {ϕ j } p j=0 can then be recovered as

ϕ j (s) = p ∑ k=0 V k j ψ k (s) = p ∑ k=0 Q k j L k (s), Q := Z V. (4.9)

Implementation of preconditioners

In order to compute the action of the left-preconditioner P , we use the matrix representation of the mapping v → (I v) in the basis given by {ϕ j } p j=0 . Thus we need to find the matrix K such that (I ϕ k ) = ∑ p j=0 K k j ϕ j for all k = 0, . . . , p. We show below that this is easily computed. Inverting (4.1) by induction yields

L k = k-1 ∑ j=0 1 -(-1) k-j j + 1 2 L j ∀ k 0. (4.10)
Therefore, we have

(I L k ) = p ∑ j=0 K * k j L j ∀ k = 0, . . . , p, (4.11)
where, after some calculation, it is found that the matrix K * ∈ R (p+1)×(p+1) can be defined in terms of K ∈ R (p+1)×(p+1) , x ∈ R p+1 , and y ∈ R p+1 by

K * := K -x y , Kk j := 1 -(-1) k-j j + 1 2 if j k, 0 otherwise,
x j := (-1) j , y k := (-1) 1-k k + 1 2 . LEMMA 4.2 Let p 2 be an integer, and let {ϕ j } p j=0 be the eigenfunctions defined by (4.9). Then, for each k = 0, . . . , p, there holds

(I ϕ k ) = p ∑ j=0 K k j ϕ j , K := V D -1/2 K * D 1/2 V, (4.12)
where V is as in (4.8), the matrix K * is as in (4.11), and where D = diag{2/(2 j + 1)} p j=0 . Proof. The orthogonality of the matrix V in (4.8) implies that ψ k = ∑ p m=0 V km ϕ m . Since the basis {ψ k } p k=0 is orthonormal in the inner product of (4.4), there holds

L j = p ∑ r=0 1 -1 (I L j ) (I ψ r ) ds ψ r = p ∑ r=0 p ∑ m=0 K * jm (r + 1 2 ) 1/2 1 -1 L m L r ds ψ r = p ∑ r=0 K * jr r + 1 2 -1/2 ψ r = p ∑ r=0 (K * D 1/2 ) jr ψ r = p ∑ m=0 (K * D 1/2 V ) jm ϕ m ,
where we have used (4.3) and (4.11) in the first line. Therefore, we compute

(I ϕ k ) = p ∑ j=0 V jk (I ψ j ) = p ∑ j=0 V jk ( j + 1 2 ) 1/2 L j = p ∑ j=0 p ∑ m=0 V jk (D -1/2 K * D 1/2 V ) jm ϕ m ,
which completes the proof. We now show how the matrix K is used in applying the preconditioner

P . For v = ∑ p k=0 v k ϕ k , there holds Pv = p ∑ j=0 A -1 Mw j + τ 2 v j ϕ j , w j := p ∑ k=0 K k j v k , F (Pv) = p ∑ j=0 f j A -1 Mw j + τ 2 v j = p ∑ k=0 MA -1 p ∑ j=0 K k j f j + τ 2 f k v k .
Therefore, g := P f, where g = (g 0 , . . . , g p ), can be computed componentwise by

g k = MA -1 p ∑ j=0 K k j f j + τ 2 f k , k = 0, . . . , p. (4.13)
The action of P requires the solution of p + 1 independent systems with matrix A, which can be performed in parallel. If we ignore communication costs, then the cost of computing g is independent of the polynomial degree p on a parallel machine with sufficiently many computing nodes. After application of the preconditioner P , the linear system has the form

Lu = g,
where L = P B, which can be solved by the PCG algorithm with preconditioner H, as suggested in section 3. The PCG algorithm requires the action of the matrices L and H -1 at each iteration. There are two ways to implement the action of L, the first being the application of B followed by P . The downside of this approach is that it leads to greater communication costs on a distributed memory parallel machine where each node holds in memory only a few coefficients of u: the matrices B and P are generally block-dense, so all vector components must be exchanged between all computational nodes for both steps. The second approach is to use (3.5), which shows that L can be expressed as a diagonal matrix plus a "rank-two" term: in the basis {ϕ j } p j=0 , we have

L = diag MA -1 M + τ 2 λ j 4 A p j=0 + τ 2 q 1 q 1 ⊗ M + τ 2 q -1 q -1 ⊗ M, (4.14)
where q ±1 := (q 0 (±1), . . . , q p (±1)) are the vectors of endpoint values of the ϕ j . Therefore, we may compute Lu, with u = (u 0 , . . . , u p ), as follows:

w j := Mu j , j = 0, . . . , p, z ±1 := p ∑ j=0 ϕ j (±1)w j , (4.15a)
(Lu) j = MA -1 w j + τ 2 λ j 4 Au j + τ 2 ϕ j (1)z 1 + τ 2 ϕ j (-1)z -1 . (4.15b)
This approach requires the same number of matrix-vector products as the first approach outlined above, but the communication costs are greatly reduced. Indeed, it is seen that the above procedure requires the parallel computation of the w j , which are then gathered and reduced by one or two nodes tasked with computing z ±1 . The results are then broadcast back to all nodes, after which all subsequent computations can be performed in parallel.

Finally, the application of H -1 is trivially parallel, as it requires only that each node solves two linear systems involving the matrix M + τ λ j /2A, and one application of A. Theorem 3.2 shows that the eigenvalues remain uniformly bounded from above, and that the smaller eigenvalues approach zero as p increases. Therefore, for large p and j, the matrix M + τ λ j /2A becomes increasingly similar to M, implying that these systems will become increasingly easy to solve in many applications.

Numerical experiments

Condition numbers

In this experiment, we study the dependence of the condition numbers κ of the preconditioned system H -1 L on the problem and discretization parameters. In order to compute accurately the condition numbers, we consider first a one-dimensional problem. Let M and A be the mass and stiffness matrices obtained by applying piecewise-linear continuous finite elements on a uniform subdivision of the domain Ω = (0, 1) into elements of size h = 2 -k , k = 5, . . . , 10. We note that M and A thus depend on h, although this is left implicit in our notation. For this experiment, we use direct solvers to implement the action of A -1 and (M + τ λ j /2A) -1 .

Table 1 shows the condition numbers κ as a function of τ for a wide range of parameters. For this experiment, we set p = 2 and h = 2 -5 . It is found that for either very large or very small τ, κ → 1. This is easily explained by the fact that

lim τ→∞ v L v v H v = lim τ→0 v L v v H v = 1. (5.1)
Therefore the condition number κ → 1 as τ → 0 or τ → ∞. Thus the maximal condition number is found for intermediate values of the time-step size τ, although in all cases it satisfies the theoretical bound κ 4 as shown in (3.20). Table 2 shows that the condition number has little to no dependence on mesh refinement, i.e. variation of M and A. As p becomes very large, the condition number approaches an asymptotic value of around 2.686 4, as shown by Table 3; our goal in testing our preconditioners here with very high polynomial degrees is merely to ascertain the asymptotic behaviour of the condition number. The proposed preconditioners are thus fully robust with respect to the parameters, in agreement with the theoretical bound κ 4 from (3.20). Furthermore, this experiment suggests that the efficiency of the preconditioners can exceed theoretical expectations, as shown by condition numbers κ 2.686 throughout these tests.

Iteration counts and multigrid preconditioning

The condition number bound (3.20) holds for the preconditioner H, which assumes that exact solvers are used to apply the inverses of the matrices M + τ λ j /2A. However, in practice, it is desirable to use τ 10 -6 10 -5 10 -4 10 -3 10 -2 10 -1 1 10 κ 1.011 1.103 1.749 2.031 2.028 2.019 1.693 1.089 Table 1. Dependence of the condition number κ of H -1 L on a wide range of values for the time-step size τ. Observe that for either very large or very small τ, the condition number is asymptotically 1, in agreement with (5.1). This explains the observation that the maximum condition number is observed for intermediate values. 3. Asymptotic behaviour of the condition number with respect to the polynomial degree p = 2 m , m = 3, . . . , 8, with fixed h = 2 -5 and τ = 0.1. The asymptotic value of κ appears to be κ ≈ 2.686. This suggests that the theoretical bound κ 4 of (3.20) may not be quantitatively sharp in all cases, although it is nevertheless very close. a cheap approximation, such as a small number of iterations from an iterative solver for M + τ λ j /2A, which leads to an approximation of the ideal preconditioner H -1 .

κ h = 2 -5 h = 2 -6 h = 2 -7 h = 2 -8 h = 2 -9 h = 2 -10 p = 1 1.
In this experiment, we study the effect of this approximation, in particular when the inverse of M + τ λ j /2A is approximated by a small number of multigrid V-cycles. Let M and A be defined as the mass and stiffness matrices of the piecewise-linear finite element space defined on a uniform triangulation of size h = 2 -k , k = 6, . . . , 10 in the two-dimensional domain Ω = (0, 1) 2 . The finest mesh thus leads to more than one million degrees of freedom (DoF) for V. We consider the system Lu = g, with a chosen exact solution u. We approximate solvers for M + τ λ j /2A by applying n V-cycles, n ∈ {1, 2, 3}, with symmetric Gauss-Seidel smoothers. Thus each application of the preconditioner requires 2n V-cycles per block. To obtain a fair comparison of all the preconditioners, a relative tolerance of 10 -6 of the true error in the energy norm was used to determine convergence of PCG, and zero initial guesses were used for all computations.

Table 4 shows the PCG iteration counts for both ideal and approximate preconditioners; in this ex- is explained by the fact that the error on the 9-th iteration was very close to, but above, the convergence criterion.

Mesh size DoF Direct 1 V-cycle 2 V-cycles 3 V-cycles h = 2 -
periment, we set p = 2 and τ = 0.1, and thus there are more than 3 million degrees of freedom for V p on the finest mesh. Fast and robust convergence of PCG is observed in all cases, showing the effectiveness of these preconditioners. For this experiment, the method using 1 V-cycle appears as the most efficient, as the reduction in cost per iteration outweighs the additional PCG iteration required. In Table 5, we study the dependence of the iteration counts on the polynomial degree. It is found that the use multigrid preconditioners retains the robustness of the preconditioner with respect to the approximation order p. This is indeed the expected result, since Theorem 3.2 shows that as p is increased, for large j the matrices M + τ λ j /2A are spectrally closer to M and thus the efficiency of the multigrid V-cycle increases.

Our main goal in this experiment is to check the robustness of our preconditioners with respect to the polynomial degree, rather than to actually propose employing high-order approximations in time coupled with low-order approximations in space. Nevertheless, such high-order temporal approximations are encountered in the context of hp-version methods (Schötzau & Schwab, 2000;[START_REF] Von Petersdorff | Numerical solution of parabolic equations in high dimensions[END_REF]Schötzau & Wihler, 2010;[START_REF] Werder | hp-discontinuous Galerkin time stepping for parabolic problems[END_REF]. We have also computed the results of this experiment using smaller time-step sizes τ, and we observed a decrease in the iteration counts, as predicted by (5.1). The conclusion drawn from these results is that in practice, standard approximate solvers can be used while retaining the key properties of the ideal preconditioner H, namely the fast convergence of the PCG algorithm and the robustness with respect to the discretization parameters.

Parabolic problems with inexact solvers

The action of L as given in (4.15) requires the action of A -1 . In many practical applications, it is desirable to approximate this step by using an inexact solver, such as a fixed number of multigrid V-cycles, leading to an approximation L ≈ L. This raises the question of whether this can be performed without affecting firstly the accuracy of the solution, and secondly the performance of the preconditioning strategy. This experiment provides evidence that inexact solvers can indeed be used without compromising these important objectives. Consider the piecewise-linear simplicial conforming finite element approximation of the heat equation in the unit square Ω = (0, 1) 2 , final time T = 0.1, imposed with initial datum u 0 (x, y) = x(1x) sin(πy) and with homogeneous Dirichlet lateral boundary conditions. This initial datum is chosen as it leads to a solution with decreased temporal regularity, see (Schötzau & Schwab, 2000). To compare near-exact and inexact solvers, we consider two approaches: 2.546 × 10 -2 3.078 × 10 -3 3.078 × 10 -3 1.303 × 10 -8 1/2 1.475 × 10 -2 3.934 × 10 -4 3.934 × 10 -4 2.129 × 10 -8 1/4 8.008 × 10 -3 5.444 × 10 -5 5.441 × 10 -5 3.219 × 10 -8 1/8 4.178 × 10 -3 8.718 × 10 -6 8.641 × 10 -6 8.126 × 10 -8 Table 6. Final time errors u(T ) -u τ (T ) L 2 (Ω ) obtained by the DG time-stepping method for the problem of section 5.3, for p = 0 (the backward Euler method) and for p = 1 with direct solvers (D) and inexact multigrid (MG) solvers. The last column gives u τ (T ) -ûτ (T ) L 2 (Ω ) , where u τ , respectively ûτ , denotes the solution obtained by (D), respectively (MG).

Iterations τ = T τ = T /2 τ = T /4 τ = T /8 (D) 4 3.5 3 3 (MG) 6 5 5 5

Table 7. Average number of PCG iterations per time-step for the problem of section 5.3, with p = 1 and with either direct (D) or inexact multigrid (MG) solvers.

(D) direct solvers are used in the application of P , H -1 and L, (MG) 5 multigrid V-cycles are used to approximate the application of A -1 in L and P , and 1 multigrid V-cycle is used to approximate (M + τ λ j /2A) -1 in H -1 , as in section 5.2.

We point out that the method (MG) does not use any direct solvers throughout the entire computation, and that each PCG iteration costs 7 V-cycles per block. Table 6 shows the final time errors u(T ) -u τ (T ) L 2 (Ω ) , for varying τ, where u denotes the exact solution of the PDE, and where u τ denotes the discrete solution. Here, we used h = 2 -8 and p = 1. For comparison, we also show the errors attained for p = 0, i.e. the backward Euler (BE) method. The inexact solvers for (MG) retain the accuracy of the method, as the difference in solutions between (D) and (MG) is several orders of magnitude smaller than the difference to the exact solution. For p = 1, the expected third order super-convergence rate is observed.

Table 7 shows the average number of PCG iterations per time-step required by both approaches in order to obtain a residual tolerance of 10 -6 , as well as the final time error between the approximate solution u τ for (D) and ûτ for (MG). This shows that the number of PCG remains robust with respect to the approximation of L entailed by (MG). This experiment shows that the use of inexact solvers can retain both the accuracy of the DG time-stepping method as well as the efficiency of the proposed preconditioners.

Conclusion

We have developed efficient and robust preconditioners enabling the fast solution of the DG timestepping method by the preconditioned conjugate gradient algorithm. The analysis and numerical experiments show that the ideal and approximate preconditioners are robust with respect to all discretization parameters and lead to low condition numbers for the preconditioned system. Thus the high-order solution of large problems by the DG time-stepping method is tractable under the proposed approach.

  + c 00 τ A • • • b 0p M + c 0p τ A . . . . . . . . . b p0 M + c p0 τ A • • • b pp M + c pp τ A

j

  FIG.2. Eigenvalues {λ j } p j=0 computed by the method of section 4.1 and arranged in decreasing order, for p = 5, 10, 100 and 100. All plots are on a logarithmic scale, except for p = 5 which is given on a semilogarithmic scale. It appears that the bounds of Theorem 3.2 are sharp.

Table 2 .

 2 Dependence of the condition number κ for p = 1, . . . , 6 and h = 2 -5 , . . . , 2 -10 , with τ = 0.1. For fixed p, the condition number is insensitive to the mesh size h, i.e. variation of M and A. The condition number slowly increases with p towards its asymptotic value of approximately 2.686, see Table3.

	p = 2 p = 3 p = 4 p = 5 p = 6	318 2.019 2.243 2.353 2.416 2.493	1.319 2.019 2.243 2.353 2.417 2.493	1.319 2.019 2.243 2.353 2.417 2.493	1.319 2.019 2.243 2.353 2.417 2.493	1.319 2.019 2.243 2.353 2.417 2.493	1.319 2.019 2.243 2.353 2.417 2.493
		p κ 2.558 2.643 2.674 2.684 2.686 2.686 8 16 32 64 128 256

Table

Table 4 .

 4 Number of PCG iterations required for convergence of the PCG algorithm in the experiment of section 5.2, for various mesh sizes and corresponding number of degrees of freedom. Fast convergence of PCG is observed in all cases, with the iteration counts being mesh-independent. The method using 1 V-cycle appears as the most efficient in this experiment.

	11 907 48 387 195 075 783 363 h = 2 -10 3 139 587 6 h = 2 -7 h = 2 -8 h = 2 -9	7 7 7 7 7	8 8 8 8 8	7 7 7 7 7	7 7 7 7 7

Table 5 .

 5 Dependence of the PCG iteration counts on the polynomial degrees p = 4, . . . , 14 in the experiment of section 5.2, with approximate solvers using 1 multigrid V-cycle. The number of PCG iterations remains uniformly bounded, showing full robustness with respect to the polynomial degree, even up to 15 million DoF. The additional iteration for the case p = 12, h = 2 -10

† Corresponding author. Email:iain.smears@inria.fr

Acknowledgements

The author wishes to express his thanks to Paul Houston, Lorenz John, Christian Kreuzer, Endre Süli and Martin Vohralík for many helpful discussions. The author was supported by an EPSRC Doctoral Prize at the Mathematical Institute, University of Oxford, during the period in which part of this work was completed.

A. Analysis of eigenvalues

In this section, we present the proof of Theorem 3.2. Without loss of generality, it is sufficient to consider only p 2 throughout this section. We will make use of Weyl's Theorem from eigenvalue perturbation theory [START_REF] Demmel | Applied numerical linear algebra[END_REF].

THEOREM A.1 (Weyl) For a positive integer n, let T and T denote symmetric n × n matrices. Let α 1 • • • α n denote the eigenvalues of T and let α1 • • • αn denote the eigenvalues of T . Then, for each j = 1, . . . , n, we have

where • 2 denotes the matrix 2-norm.

We are now ready to prove Theorem 3.2. Proof of Theorem 3.2. Recall that the matrix T is defined by (4.7). We start by noting that the matrix T is pentadiagonal, i.e. T k j = 0 if |k -j| > 2. Moreover, it is easy to show from the properties of Legendre polynomials that there exists a constant C independent of p such that

The Cauchy-Schwarz inequality implies that |T jk | T kk T j j , and thus all entries of T are uniformly bounded. Since T has at most 5 nonzero entries per row, and all entries are uniformly bounded, the Gershgorin discs of T are bounded independently of p. The Gershgorin Disc Theorem therefore implies that there is a constant C, independent of p, such that λ 0 , the maximal eigenvalue of T , satisfies λ 0 C. This corresponds to (3.13) for the case j = 0. We consider now j 1, and without loss of generality, we may assume that p 2. For 0 q p -1, we define the orthogonal projector π q : P p → P q by

where the polynomials ψ j are defined by (4.2), and v = (v 0 , . . . , v p ) ∈ R p+1 . Define the matrix Tq ∈ R (p+1)×(p+1) by ( Tq ) k j := (π q ψ k , π q ψ j ) L 2 . Note that ( Tq ) k j is zero if either k or j is greater than q, and equals T k j otherwise. Thus, the matrix Tq has the general form

where T denotes the (q + 1) × (q + 1) principal submatrix of T . The main step of the proof is to use (A.2) repeatedly to find an upper bound for

For arbitrary u and v ∈ R p+1 \ {0}, define the polynomials u := ∑ p j=0 u j ψ j and v = ∑ p j=0 v j ψ j . Then, we have

It follows from the definition of π q in (A.3) that

Note that this implies the following a priori estimate for the orthogonal projector:

Thus (A.6) and the Cauchy-Schwarz inequality imply that

For k q and j q + 1, we have (ψ k , ψ j ) = 0 if k < q -1 and j > q + 2, and thus

The Cauchy-Schwarz inequality thus implies that

(A.9)

An identical argument shows that

Combining (A.8), (A.9) and (A.10) implies that there exists a constant C, independent of p and q, such that T -Tq 2 C(q + 1) -2 . (A.11)

Therefore, Weyl's Theorem implies that the eigenvalues λ j of T and λ j of Tq satisfy

where the constant C is independent of j, q and p. However, it is clear from (A.4) that λ j = 0 for j q + 1, and thus there exists a constant C independent of p and q such that

The left-hand side of this inequality is independent of q while the right-hand side is valid of all q j -1. Therefore, the choice q = j -1 yields (3.13) for j = 1, . . . , p.

The lower bound (3.14) follows from inverse inequalities. Indeed, λ p satisfies

.

To obtain (3.14), it is therefore enough to show the inverse inequality

For any v ∈ P p , standard inverse inequalities [START_REF] Schwab | p-and hp-finite element methods[END_REF] imply that


Josiena Gotzsch
email: jose.gotzsch@grenoble-em.com

Jean-Jacques Chanaron
email: jean-jacques.chanaron@grenoble-em.com

David Birchall
email: david.birchall@henleymc.ac.uk

Product development with a focus on attractive product expression: an analysis of case studies

Creating outstanding products is vital for a company's endurance in competitive markets. A mix of functionality, ergonomics, aesthetics, symbols and price aspects all play a role in making a product desirable.

Some products carry a personal meaning for its user. Others communicate its user's identity or the company's brand image. This study concentrates on communicative and meaningful aspects in a product's design. It examines how the creation of a communicative design occurs during the new product development process. The present research has an exploratory nature. For the field research, the case study method was chosen and a guideline for semi-directed interviews developed. This interview guideline was used to analyse

multiple product development projects in two distinct companies. In two projects, specific attention was given to communicative aspects in the product's design. These two projects are examined in this paper.

In the two selected projects, the product development teams carefully studied the users' preferences for aesthetics and product messages. A user-centred approach was used in both development processes. The choice to purposely improve the communicative value of the product's design was on the one hand influenced by the limited possibilities to create other advantages, such as improving the product's functionality or reducing its cost price. A weakness in the competitor's design strategy allowed the successful improvement of the product's communication on the other hand.

Product development with a focus on attractive product expression: a case study analysis 1 Introduction

Products can transmit messages that are meaningful to users, for example, by communicating status or by evoking personal memories. Different research approaches are possible to examine product development activities leading to a communicative product design. A first research approach consists of categorising the different messages that a product allows to communicate. A second research approach might analyse the methods used to identify the users' preferences for communicative product aspects. In the third place, research can focus on examining how to incorporate these requirements for communicative value in a product's design during the development process. This article gives attention to both the second and third research theme. It examines how two companies identified the preferences for product's appearance. It also investigates how these companies incorporated these user requirements in the product's design.

Literature

Terminology

In the literature, multiple terms are used to describe the communicative aspect of a product's design. Some terms directly refer to the product. This phenomenon is then called "product language" [START_REF] Bürdek | Design : Geschiedenis, theorie en praktijk van de produktontwikkeling, Uitgeverij ten Hagen & Stam, The Netherlands[END_REF][START_REF] Steffen | Design als Produktsprache[END_REF], "product meaning" [START_REF] Muller | Vormgeven: ordening en betekenisgeving, Uitgevery Lemma[END_REF], "product soul" [START_REF] Durgee | Product soul[END_REF] or "product charisma" [START_REF] Gotzsch | Designed to be Loved: the Art of Product Charisma[END_REF].

Another part of these terms refer to the reaction that a product elicits. The term "product emotion" is proposed by [START_REF] Desmet | Designing Emotions[END_REF] and "emotional domain", "soft functions" by [START_REF] Lebbon | The Emotional Domain in Product Design[END_REF]. Meaningful expression might also give the product an "added emotional value" [START_REF] Desmet | Designing Emotions[END_REF], provide the user with a "product experience" [START_REF] Marzano | Suffusing Design Through the Organisation[END_REF] or make the product "desirable" [START_REF] Sanders | Converging Perspectives: Product Development Research for the 1990s[END_REF].

Propositions for research

The design management and new product development literature do not specifically examine the creation of communicative design. However, when analysing the literature, several aspects appear relevant for the incorporation of messages in a product's design. These research propositions are summarised in table 1 and discussed in the following text.

proposition 1 Competition is a major driver for the creation of a communicative product appearance proposition 2 A positive attitude towards design in the organisation is important for the creation of a communicative product appearance.

proposition 3 A cooperative, multifunctional team with appropriate and sufficient design skills is important for the creation of a communicative product appearance.

proposition 4 User focus in the pre-development stage to understand the users' preferences for the product's aesthetics and significance is particularly important for the creation of a communicative product appearance.

proposition 5 As part of a quality design process, the identified preferences for "product expression" have to be translated into a communicative product appearance.

Table 1: The research propositions

In the first place, competition can act as a driver (see table1, proposition 1) for the creation of a communicative design. Cooper and Kleinschmidt (1995) established the importance of the superiority of a product as key success factor. [START_REF] Cote-Colisson | Les PMI françaises et le design[END_REF] and [START_REF] Potter | The Benefits and Costs of Investment in design, Design Innovation Group[END_REF] specifically validate the role of design in differentiating a product from competition. [START_REF] Marzano | Suffusing Design Through the Organisation[END_REF] and [START_REF] Lebbon | The Emotional Domain in Product Design[END_REF] confirm that with more products available, consumers specifically attribute importance to the product's emotional quality.

In other words, the presence of competition pushes a company towards the creation products that differentiate from competition. Incorporating communicative qualities into a product's design allows such a differentiation.

In the second place, the literature indicates that the design management process is helped by an open-minded attitude towards design at all levels in the organisation [START_REF] Bruce | Marketing and Design Management[END_REF][START_REF] Cote-Colisson | Les PMI françaises et le design[END_REF] and an innovative culture with openness to change and encouragement for innovation [START_REF] Bruce | Marketing and Design Management[END_REF]). An positive attitude towards design

(see table 1, proposition 2)
is, therefore, expected to facilitate specific aspects of the design process, such as the creation of products with communicative value Additional factors found to facilitate the design management process are experience with design and care for product details [START_REF] Bruce | Marketing and Design Management[END_REF], importance given to aesthetics [START_REF] Freeze | Through the Backdoor: The Strategic Power of Case Studies in Design Management[END_REF], a strong company or brand identity and a drive for continuous product improvement [START_REF] Hart | The Effects of Managerial Attitudes to Design on Company Performance[END_REF].

The commitment of people involved in the project and design professionalism [START_REF] Bruce | Marketing and Design Management[END_REF], small, close-knit, multi-disciplinary teams [START_REF] Hart | The Effects of Managerial Attitudes to Design on Company Performance[END_REF][START_REF] Bruce | Marketing and Design Management[END_REF] with regular and good communication, mutual respect and harmony were found as key factors for the quality of the development project. The design skills within a cooperative multifunctional team (see table1, proposition 3), therefore, seem to be important for the creation of products with communicative value.

Cooper and Kleinschmidt (1995) highlight the importance of pre-development activities in product development. A "customer focused process" and "homework" that is well done before the development process, are crucial. The importance of an early user-orientation is also pointed out by [START_REF] Hart | The Effects of Managerial Attitudes to Design on Company Performance[END_REF] and [START_REF] Bruce | Marketing and Design Management[END_REF]. To develop communicative product design, it therefore seems vital, to understand the users' preferences for product appearance (see table1, proposition 4).

Decision-making qualities are important for a high quality development process (Cooper and Kleinschmidt, 1995). This process must be protected by high positions in the organisation [START_REF] Bürdek | Design : Geschiedenis, theorie en praktijk van de produktontwikkeling, Uitgeverij ten Hagen & Stam, The Netherlands[END_REF]. Companies that create products with communicative value, probably intend to create the best possible product for their users. It is, therefore, assumed that they use a quality design process to incorporate the users' preferences for product expression into the design (see table1, proposition 5).

The above research propositions were specifically examined in the field research phase of which the methodology is described in the following section.

Methodology

This study investigates how the creation of a communicative design takes place. It has an exploratory character. The case study method was selected as an appropriate approach and a structured interview guide with open and semidirected questions was developed. Two pilot studies took place to test this interview guide. This allowed confirming the duration of the interviews and verifying details such as the use of the taperecorder to register the interviews. Six in-depth case studies followed. The type of products that were selected for the case studies were consumer products for use in a "home environment", such as products used in the living room, kitchen or bathroom. From the above category, products were chosen with which a consumer normally is in regular physical contact, because this makes the affective contact with these products more important. The companies are operating in competitive markets.

Interview structure

For each project, several designers, design managers or product / marketing managers were interviewed. Before the interviews, publications about the projects were collected. The case studies were written, based on the information from the interviews (kept on notes and audio tapes) and completed with the additional publications. The case studies describe the product development process and focus on conditions and activities in the product development process contributing to the creation of product expression. The case study description depended on the information that was obtained from the multiple sources. Each case description has, however, the following structure:

project origin, product expression, product development process and future of the project.

The interviewees were asked to verify the written summary. A double check was made by asking other participants in each of the six projects to review the text.

In two of the six projects, it appeared that the companies specifically aimed to integrate an attractive user-centred expression in the design. These are the two projects that are described and analysed in this article.

Results

The two selected projects originated in very different companies. One is a multinational company with an important in-house design team. The other is a small independent design consultancy.

Company A has its headquarters in the Netherlands and is, with approximately 10,000 employees, part of a much larger multinational group. Company A commercialises small domestic consumer products and works with an internal design team. This internal design team is part of the company's design department with 500 employees. Between 1997 and 1999, Company A developed a range of small kitchen appliances (see figure 1).

Influence of competition

Proposition 1: Competition as a major driver

Competition pushed the multinational Company A towards product differentiation. Company A's market is very competitive and the success of this range of kitchen appliances was of strategic importance. Design is a major sales argument for small kitchen appliances within this price segment.

Improving the functionality of these "mature" products is difficult. Improving the products' appearance, however, still offers opportunities. An extensive international trend and market research took place to identify the visual preferences of its user group. The users of these small kitchen appliances were described as "quality-oriented traditionals". They were portrayed as mature, high-end, quality-oriented consumers that appreciate design, but that are quite suspicious when confronted with radical changes. The products had to be specifically appealing to this group of users. During the product development process, colours and shapes that refer to traditional kitchen materials such as ceramics and metals were selected. The products were designed to communicate good quality. They had to look modern, but not revolutionary.

A search for product differentiation also motivated the second company in this case study, the independent design consultancy (company B), to propose a "lively" toothbrush. This product has a small base "with feet" and was specifically designed to please children (see figure 2, left -front). At the time, children's toothbrushes used to look similar to those for adults. This very different "playful" toothbrush design was a commercial success, since it increased the client's European market share from 8% to 42% in 1991.

Competition simply ignored this opportunity at the time.

The design team continued to make the design of the three toothbrushes attractive for children. A second toothbrush "with arms" replaced the first design in 1994 (see figure 2, left -back). The brush itself remained unchanged, but its support had feet and arms integrated into its design. A third toothbrush with a "suction disc" that adheres to tiles was launched in 1997 (see figure 2, right).

To be competitive, both companies needed a desirable product. They used a user-centred approach to create a product appearance, intended to be specifically appealing to their clearly identified users. Creating an appearance dedicated to such a user group proved to be a good strategy, especially because competition underestimated or ignored this possibility at that moment.

Company culture and design skills

Proposition 2: The importance of a positive attitude towards design

Within Company A, design has a high standing. The marketing department works closely with the internal designers. Company B's client, on the other hand, had little experience of working with product designers, but its marketing manager was open to design solutions.

The creation of a meaningful design, therefore, proves possible in a company that does not place design on a strategic level (yet), but that has an innovative open-minded attitude towards design. Design experience, however, might become crucial when competition is also focusing on this design aspect.

Proposition 3: Design skills in a cooperative, multifunctional team

Company A has a product division specifically responsible for coffee-makers and another one responsible for food-preparation products. Both divisions were in financial difficulties around 1997. A new Business Unit Manager encouraged changes in the organisation in an effort to overcome these problems. This stimulated the marketing managers of the coffee-makers and the foodpreparation products to develop a common marketing strategy. The design department also used a new and innovative approach by working with trend watchers in New York, Paris and Tokyo. The objective was to obtain a clear vision of trends in colours and styling that could appeal to the user group. This mixture of innovative marketing people, quality design skills, and new research in marketing and design appeared decisive in this project.

In Company's B toothbrush case, the role of the design consultancy team and their client's marketing manager was important. A small team with creative ideas and a limited budget worked on the project.

In a multifunctional team, the cooperation between marketing and design executives appears a real condition for a successful project. The professional quality of the team members is important to define and obtain a precise "feeling" of what is needed in the market and to incorporate this into the product's design. These characteristics, however, seem to be fundamental for every product development project and not only for the creation of a communicative design.

Product development activities

In the projects, different methods were used to identify the product appearance that users might prefer. Company A did extensive market and trend research before the development phase. Company B used spontaneous idea generation and involved users (in this case children) during the creation process.

Proposition 4: understanding the users' preferences

Company A analysed its market and obtained profound knowledge of its target consumer before the development process. Styling characteristics that would appeal to these users today and in 3 years time were examined in different countries by trend analysts. For company A, it was new to place so much emphasis on market and trend research. Company B used a much "lighter" user-centred approach. During a brainstorming session, the idea was developed to create the product in a style that appealed to children.

A user-orientation is necessary to understand preferences and sensibilities of the product users. Serious competition appears to create the need for profound and extensive user-focused research.

Proposition 5: the importance of a quality product development process

How are the user's visual preferences translated into product design? Company A changed the organisation of its project teams, to enable one dedicated team to develop the entire product range for a specific user group. Before this reorganisation, teams were specialised in developing one specific type of product. A marketing manager was made responsible for the coherence in quality and appearance of the total product range. The project team felt confident during the product development phase. They sensed that the initial market and trend analysis studies were well done. Hardly any product testing was used.

Company B "only" used the team's creativity to create the first toothbrush.

Three years later, competition also started commercialising attractive toothbrushes for children. As a response to growing competition, a budget for research was made available and children were involved in the project for idea development.

The methods used by the companies differ. Company

Discussion

Conclusions from the cases

It appears that there is evidence to support the claim that competition acts as a major driver to develop products with communicative value. A product appearance having meaning to its user seems an effective approach to add value to the product. The aspects in a product's design, that a company might choose to pay specific attention to during the product development process, depend on multiple aspects. This choice depends on many variables; the users' preferences, the degree of competition in the market, the competitor's design strategy, and on other possibilities to integrate advantages, such as functional improvements, technical excellence and price aspects.

A positive attitude to design (proposition 2) and design skills in a cooperative multifunctional team (proposition 3) appear very important for a successful product creation process, but are not specific for the development of communicate design.

However, to develop a product appearance that is attractive to its users, it appears crucial to have a creative focus on the users' preferences during the pre-development phase and the creation process (proposition 4 and 5).

Research in this pre-development phase might make use of different techniques such as brainstorming, market and trend research, depending on the budget and the degree of competition. In the creation process, user involvement can be a source of additional ideas for the creation of a meaningful design.

Limitations and further research

The conclusions are based on two successful product development projects.

These two projects were chosen from the six case studies, because these two specifically focus on creating an appearance having meaning for the product user. The other case studies aimed at differentiating the product by using other methods, such as improving the product's functionality or concept. This number of cases is too limited to generalize the findings. The two companies in the case studies are also very different, because one case concerns the in-house design practice of a multinational company (company A) and the other involves a small independent design consultancy (company B), working for a multinational client. It will be a next step in research to examine how other companies handle the development process of communicative meaningful product designs. It is an important process, because it is an additional, underestimated way for product differentiation that allows bringing commercial success in competitive markets.

  Key words: User-centred Product, Design Management Process, New Product Development Process, Development Product Expression, Product Meaning, Product Language, Product Value Biographical notes: Josiena Gotzsch is Associate Professor in the Department of Technology Management at the Grenoble Ecole de Management (GEM) and Program Director at the Grenoble Graduate School of Business (GGSB) for the international Undergraduate Program, France. Josiena Gotzsch obtained a Doctorate of Business Administration (DBA) at Henley Management College, UK and a Master of Science (MSc) in Industrial Design Engineering from the Delft University of Technology, the Netherlands. She teaches Industrial Design, the Product Development Process and Technology Management. Research interests include product design, product communication (aesthetics and symbolic value), product attachment, design management and product development. Professor J.J. Chanaron is currently Research Director within the French National Centre for Scientific Research (CNRS) and Chief Scientific Advisor at the Grenoble Graduate School of Management. Jean-Jacques has published extensively via books, articles in refereed journals and conference papers in Industrial Economics, Economics of Innovation and Technology Management since 1973 when he received his PhD at the University of Grenoble. He also holds a HDR in Economics since 1994. He is Associated Professor and Researcher with Henley Management College, Manchester University and Newcastle University in the UK as well as Tongji University in Shanghai, China. He is a well-recognized expert in the automotive industry. He is consultant to International Organizations (EU, OECD, ILO, UNIDO), professional organizations (CCFA, FIEV, JAMA, CLEPA), OEMs (PSA, Renault, Toyota, Nissan, DaimlerChrysler, VW) and numerous component manufacturers. He is a member of the French Society of Automotive Engineers (SIA) and the GERPISA International Network of Researchers on the Auto Industry. In April 2004, he has been granted the IAMOT award for research excellence in Technology and Innovation Management. David Birchall Professor David Birchall is Director of the School for Management Knowledge and Learning and Henley Learning Advisory Service at Henley Management College, UK. Professor Birchall teaches and researches aspects of work in the digital economy, including innovation management. He currently directs projects funded by the European Union, UK Government Department and several commercial organisations. His latest book Capabilities for Strategic Advantage: Leading Through Technological Innovation -co-authored with George Tovstiga -was published by Palgrave in May 2005

Figure 1 :

 1 Figure 1: Products of the kitchen appliances range

Figure 2 :

 2 Figure 2: Three generations of toothbrushes (left) and the toothbrush with suction cup (right)photograph: B. Moyen -Lucy in the Sky, Annecy

  A committed a project team to the product range. The organisation moved from an organisation based on a product type to one based on the entire product range dedicated to specific users. The team used the information from the pre-development research. Company B did no pre-development research for the design of the first successful toothbrush. They used a simple creative user-centred approach within the development team. They imagined what would appeal to children and their parents. But when a few years later competition intensified, users (children) became involved in the creation process.

Table 2 :

 2 Structure of interview guideThe interview started with general information (see table 2, A) about the company and the interviewee, such as the interviewee's background, role in the project and experience in design. The next theme concerns the company culture (table 2, B) and the position of design in the company. The following questions focused on a specific project in the company (table 2, C). The history of this project and its product development process were examined in open questions. Detailed questions about the pre-development and product creation activities followed (table 2, D, E, F, G, H, I). The questions in the interview guide were directly related to the research propositions. This is indicated as p1 to p5 (proposition 1 to proposition 5) in table 2.

			Theme	Main question	Details
	A		General information	About the company	Name, size, products, …
				About the interviewee	Role in project,
					background in design,
					…
	B	p2	Company culture	Design in your company	Position, origin, attitude
				Innovative culture
	C	p1	Product / project choice Short product history	When developed,
				Strong & weak points	competition…
	D		Development process	Open question
			Expression used
	E	p4	Pre-development	Method to identify need for	Source, objective, who,
			process	expression	…
	F	p5	Creation process	What took place	Who, how, decision
					process, …
	G	p3	R&D resources	People
				Time (budget)
	H		Facilitating &	Which aspects facilitated
			constraining aspects	or constraint the process
	I		Results	Consumer reaction	Expected / surprise
				Product personality:	Product messages
				details	Market position,
				Commercial results	financial, other
				Company learning process	results, …