Runge-Kutta Preconditioners

Martin J Gander

Michal Outrata

SPECTRAL ANALYSIS OF IMPLICIT S-STAGE BLOCK

Keywords: implicit Runge-Kutta methods, stage equations, preconditioned GMRES, convergence estimates, conformal maps MSC codes. 65L06, 65F10, 65E05

 for the stage equations of implicit Runge-Kutta methods for two stage methods. We simplify the formulas for the eigenvalues and eigenvectors of the preconditioned systems for a general s-stage method and use these to obtain convergence rate estimates for preconditioned GMRES for some common choices of the implicit Runge-Kutta methods. This analysis also allows us to qualitatively predict and explain the main observed features of the GMRES convergence behavior and we illustrate our analysis with numerical experiments.

1. Introduction. Runge-Kutta methods are a well-established family of onestep solvers for systems of ordinary differential equations (ODEs; see [START_REF] Wanner | Solving Ordinary Differential Equations I : Non-Stiff Problems[END_REF][START_REF] Wanner | Solving Ordinary Differential Equations II : Stiff and Differential-Algebraic Problems[END_REF] for an overview and further references). For implicit methods (IRK), their efficiency depends on the efficiency of a solver for the so-called stage equations -in general a system of ms non-linear equations, where m is the number of scalar ODEs in the system and s is the number of stages of the Runge-Kutta method. An important application arises from the space discretization of time-dependent partial differential equations (PDEs), resulting in a system of ODEs with very large m. If the spatial operator is linear, then the stage equations also form a system of linear algebraic equations and are often solved by an iterative solver, e.g., a Krylov method. In [START_REF] Rana | A New Block Preconditioner for Implicit Runge-Kutta Methods for Parabolic PDE Problems[END_REF], the authors introduced a family of preconditioners for GMRES for the stage equations, numerically showing that these preconditioners give an outstanding performance, especially under refinement of the spatial mesh, i.e., as m grows. Recently, there have also been other contributions in the direction of preconditioning the fully implicit Runge-Kutta stage equations for PDEs, see [START_REF] Southworth | Fast solution of fully implicit Runge-Kutta and discontinuous Galerkin in time for numerical PDEs, Part I: The linear setting[END_REF][START_REF] Southworth | Fast solution of fully implicit Runge-Kutta and discontinuous Galerkin in time for numerical PDEs, Part II: nonlinearities and DAEs[END_REF] but also [START_REF] Neytcheva | Numerical Solution Methods for Implicit Runge-Kutta Methods of Arbitrarily High Order[END_REF][START_REF] Munch | Stage-parallel fully implicit Runge-Kutta implementations with optimal multilevel preconditioners at the scaling limit[END_REF] and [START_REF] Clines | Efficient order-optimal preconditioners for implicit Runge-Kutta and Runge-Kutta-Nyström methods applicable to a large class of parabolic and hyperbolic PDEs[END_REF], introducing new ideas in terms of implementation as well as formulation and testing these numerically on a variety of test problems.

We focus on the setting considered in [START_REF] Rana | A New Block Preconditioner for Implicit Runge-Kutta Methods for Parabolic PDE Problems[END_REF], expand the 2-stage method analysis given in [START_REF] Gander | Spectral analysis of implicit 2-stage block Runge-Kutta preconditioners[END_REF], and consider the general s-stage case, giving a theoretical background for the performance and spectral properties observed. First, we recall some important preliminaries in Section section 2 so that we can deliver the analysis, based on the spectral analysis of the preconditioned system, in Section section 3. We support the analysis by considering more involved examples in Section section 4.

Model problem and preliminaries.

As our model problem we consider the heat equation on the unit square and a time interval (0, T end ), i.e., where ∆ is the Laplace operator, f, g, u 0 are given functions and Ω is the unit square Ω := (0, 1) × (0, 1). As in [START_REF] Gander | Spectral analysis of implicit 2-stage block Runge-Kutta preconditioners[END_REF] we discretize in space using a finite difference scheme on an equidistant grid with N + 1 rows and columns, and with mesh size h = 1/N as in Figure 1. The values at the interior grid points become unknown functions of time, which are governed by the system of ODEs (2.2)

∂ ∂t u i (t) = u i-N (t) + u i-1 (t) -4u i (t) + u i+1 (t) + u i+N (t) h 2 + b (ST ) i (t)
      , I =       1 . . . . . . 1      
, where L is of dimension n := (N -1) 2 and the blocks T, I are of dimension N -1.

We discretize [0, T end ] with M T end + 1 equidistant time points with time step τ =

T end /M T end , i.e.,

{0 = t 0 < • • • < t M T end = T end }, τ = T end M T end and t m = τ • m, m = 0, . . . , M T end .
This manuscript is for review purposes only.

Having a Butcher tableau

(2.5) c A b := c 1 a 1,1 . . . a
u m = u m-1 + τ s i=1 b i k m i ,
where the vectors k m 1 , . . . , k m s ∈ R n are the solutions of the linear system (2.7)

      I . . . I    - τ h 2    a 1,1 L . . . a 1,s L . . . . . . . . . a s,1 L . . . a s,s L       ≡Is⊗In-τ h 2 (A⊗L)=:M k m =    1 h 2 Lu m-1 + b (ST) (t m-1 + c 1 τ ) . . . 1 h 2 Lu m-1 + b (ST) (t m-1 + c s τ )    , with k m := k m 1 • • • k m s T ∈ R ns .
The symbol ⊗ stands for the Kronecker product (see [START_REF] Van Loan | The ubiquitous Kronecker product[END_REF] and references therein) and we note that (2.7) can be reformulated into a matrix equation, which is in general better suited for using a Krylov solver (see [START_REF] Palitta | Optimality properties of Galerkin and Petrov-Galerkin methods for linear matrix equations[END_REF]). Here we focus on the analysis of the results in [START_REF] Rana | A New Block Preconditioner for Implicit Runge-Kutta Methods for Parabolic PDE Problems[END_REF] and thus we do not address this any further but a study of the preconditioners from [START_REF] Rana | A New Block Preconditioner for Implicit Runge-Kutta Methods for Parabolic PDE Problems[END_REF] in the matrix equations setting seems worthwhile. Having p ≤ 2s as the order of convergence of the IRK method we assume that it is balanced with the spatial discretization error, i.e., that h 2 = C e τ p for some C e > 0.

The problem (2.7) with the sparse system matrix M can be very large for h (and τ ) small, suggesting an iterative solver such as GMRES, BiCG or GCR should be used, which in turn requires a preconditioner to attain efficiency. In [START_REF] Rana | A New Block Preconditioner for Implicit Runge-Kutta Methods for Parabolic PDE Problems[END_REF], the authors introduce the block preconditioners (2.8)

P d = I s ⊗ I n - τ h 2 diag(A) ⊗ L, P u = I s ⊗ I n - τ h 2 D A U A ⊗ L and P l = I s ⊗ I n - τ h 2 L A D A ⊗ L,
where L A , D A , U A are the LDU factors of the Butcher tableau matrix A. In addition, the authors also consider the block triangular preconditioners (2.9)

P GSL = I s ⊗ I n - τ h 2 A L ⊗ L and P GSU = I s ⊗ I n - τ h 2 A U ⊗ L,
where GSL/GSU stands for Gauss-Seidel lower/upper, and A L,U is the lower/upper triangular part of A, i.e.,

(A L ) ij = a ij if i ≥ j 0 otherwise , (A U ) ij = a ij if i ≤ j 0 otherwise .
Some of these -P d and P GSL -were considered already in [START_REF] Staff | Preconditioning of fully implicit Runge-Kutta schemes for parabolic PDEs[END_REF]. Notice that if a ii > 0 for all i = 1, . . . , s, then the preconditioners are invertible as L is symmetric, negativedefinite. More general conditions for non-singularity of the preconditioners can be also derived analogously to [24, Lemma 1].

Using GMRES for a linear system Cx = f with C being diagonalizable, i.e., C = SΛS -1 and Λ = diag(λ 1 , . . . , λ d ), a standard convergence bound for the residuals

r ℓ reads (2.10) ∥r ℓ ∥ ∥r 0 ∥ ≤ κ(S) min φ(0)=1 deg(φ)≤ℓ max 1≤i≤d |φ(λ i )|,
where κ(S) is the 2-norm condition number of the matrix S, see, e.g., [START_REF] Liesen | Krylov Subspace Methods: Principles and Analysis[END_REF]Section 5.7.2]. We highlight some aspects of the bound (2.10) that are often used to study GMRES convergence behavior.

Remark 2.1. As indicated above, the spectral information of the system matrix in GMRES (in our case of the preconditioned system) does not generally govern the convergence (see [START_REF] Greenbaum | Matrices that generate the same Krylov residual spaces[END_REF], [START_REF] Greenbaum | Any nonincreasing convergence curve is possible for GMRES[END_REF] and [START_REF] Arioli | Krylov sequences of maximal length and convergence of GMRES[END_REF] and also [START_REF] Liesen | Krylov Subspace Methods: Principles and Analysis[END_REF]Chapter 2 and 5.7] and the references therein). If the system matrix is normal, i.e., it is diagonalizable with S unitary, then the spectral information is enough to evaluate the ideal GMRES bound (2.10).

However, if C is non-normal, then a convincing argument needs to be put forward to validate linking spectral information with the convergence behavior of GMRES as the authors in [16, p. 303, Remark 1] point out.

Moreover, particular knowledge of the interaction of S and the initial residual r 0 can lead to a qualitative and quantitative improvement on (2.10), see, e.g., [START_REF] Liesen | GMRES convergence analysis for a convection-diffusion model problem[END_REF]. However, studying GMRES behavior with the bound (2.10), this interaction is completely lost.

In cases where (2.10) is justifiable, the next step is usually to bound from above the mixed 1 min-max problem in the right-hand side of (2.10) by replacing the discrete set over which we take the maximum, let us denote it by σ discr , by a non-discrete one, which we denote by σ non-discr , so that we have σ discr ⊂ σ non-discr . We highlight two important aspects of this step:

(

V := [v 1 , . . . , v n ] , (3.1)
Q :=    V . . . V    ∈ R sn×sn . and
Transforming M blockwise into the V basis gives M : with Λ = diag(λ 1 , . . . , λ n ). With the preconditioners proposed in (2.8-2.9) we write the spectrum of the preconditioned system as sp(M P -1 ) = sp(Q

= QM Q T , (3.2
T M P -1 Q) = sp(Q T M QQ T P -1 Q) = sp M P -1 ,
where P := Q T P Q stands for one of the right-preconditioners P d,GSU,u and an analogous formulation follows also for the left-preconditioners P l,GSL . As the preconditioners are defined blockwise as scalar multiplications of L and I, their blockwise transformation into the eigenbasis of L is a straight-forward calculation -replacing L with Λ (and keeping I). Next, such matrices -block matrices with each block being a square, diagonal matrix -can be permuted into classical block-diagonal matrices as the following lemma shows.

Lemma 3.1 (see [START_REF] Gander | Spectral analysis of implicit 2-stage block Runge-Kutta preconditioners[END_REF]Lemma 1]). Let C ∈ R ns×ns be a real matrix with block structure such that every block is a square diagonal matrix, i.e.,

(3.3) C =    Λ 11 . . . Λ 1s . . . . . . . . . Λ s1 . . . Λ ss    , with Λ ij = diag λ (ij) 1 , . . . , λ (ij) n ∀ij.
This manuscript is for review purposes only.

Then there exists a permutation matrix Π ∈ R ns×ns such that

(3.4) Π T CΠ =    C 1 . . . C n    with C ℓ =     λ (11) ℓ . . . λ (1s) ℓ . . . . . . . . . λ (s1) ℓ . . . λ (ss) ℓ     ∈ R s×s ,
for any ℓ = 1, . . . , n.

Hence, C is diagonalizable if and only if C ℓ is diagonalizable for all ℓ = 1, . . . , n,

and if C ℓ = V -1 ℓ D ℓ V ℓ is the eigendecomposition of C ℓ with D ℓ = diag(µ (1)
ℓ , . . . , µ

(s) ℓ ), then sp(C) = n ℓ=1 s i=1 µ (i) ℓ , and if (µ, v) is an eigenpair of some C ℓ , then µ, Π T (v ⊗ e ℓ ) is an eigenpair of C. As a result, if C is diagonalizable with C = V -1 DV , then κ(V ) = max ℓ=1,...,s κ(V ℓ ),
where κ(•) is the 2-norm condition number.

Remark 3.2. We note that an analogous lemma to Lemma 3.1 can also be formulated for non-normal matrices (replacing Q T by Q -1 ). Considering the Jordan canonical (or the Schur decomposition form) of C ℓ , Lemma 3.1 can be reformulated to obtain a block upper bi-diagonal (or block upper-triangular) matrix.

To shorten the notation we set

(3.5) θ k := τ h 2 λ k and Θ := τ h 2 Λ,
as these quantities appear always together in the computations. By a direct calculation (see [START_REF] Outrata | Schwarz methods, Schur complements, preconditioning and numerical linear algebra[END_REF]Appendix B.8]) we get the limit behavior of θ k as τ, h → 0,

(θ n , θ 1 ) → (- 8 C e , 0), (θ -1 1 , θ -1 n ) → -∞, - C e 8 , (LIM) p=1 (θ n , θ 1 ) → (-∞, 0), (θ -1 1 , θ -1 n ) → (-∞, 0). (LIM) p>1 (3.6)
Next we define the s-by-s matrices

M k :=       1 -a 11 θ k -a 12 θ k . . . -a 1s θ k -a 21 θ k 1 -a 22 θ k . . . . . . . . . . . . -a s1 θ k . . . . . . 1 -a ss θ k       and P ⋆ k :=       1 -α 11 θ k -α 12 θ k . . . -α 1s θ k -α 21 θ k 1 -α 22 θ k . . . . . . . . . . . . -α s1 θ k . . . . . . 1 -α ss θ k      
, where α ij are the entries of the replacement for A in M , e.g., taking ⋆ = d we have

α ij = a ij for i = j and α ij = 0 otherwise, while taking ⋆ = u we have α ij = (D A U A ) ij
where A = L A D A U A is the LDU factorization of A and so on. Using Lemma 3.1, we obtain the following result.

This manuscript is for review purposes only.

Proposition 3.3. Take M as in (2.7) and a preconditioner P from (2.8, 2.9).

Assuming P is invertible, the spectrum of M P -1 (or P -1 M ) is given as the union of the spectra of the matrices X k given by

(3.7) X ⋆ k := M k (P ⋆ k ) -1 (or P -1 k M k ), for k = 1, . . . , n. If all X ⋆ k are diagonalizable with (3.8) (S ⋆ k ) -1 X ⋆ k S ⋆ k = diag(ξ (k) 1 , . . . , ξ (k) s ),
then the condition number of the matrix of the eigenvectors of the preconditioned system is given by

κ (W ) • max k=1,...,n κ (S ⋆ k ) .
If the θ k have multiplicity at most m, then the eigenvalues of the preconditioned system have algebraic multiplicity at most ms. In particular, the preconditioned system can be non-diagonalizable but the longest Jordan vector chain has length at most ms.

Proof. Transforming M P -1 (or P -1 M ) into the basis of Q we use Lemma 3.1

for the matrix M P -1 (see (3.2)) and obtain the result. Now we are ready to generalize the results shown in [START_REF] Gander | Spectral analysis of implicit 2-stage block Runge-Kutta preconditioners[END_REF] for s = 2 to a general s-stage method.

Corollary 3.4 ([19, Proposition 7.5]). Under the assumptions of Proposition 3.3, we have for the right-preconditioner P d the formula (3.9)

X d k =        1 -a12θ k 1-a22θ k . . . -a1sθ k 1-assθ k -a21θ k 1-a11θ k 1 . . . . . . . . . . . . -a1sθ k 1-a11θ k . . . . . . 1        , with the characteristic polynomial p (s) k (λ) = (1 -λ) s + β s-2 (1 -λ) s-2 + β s-3 (1 -λ) s-3 + . . . + β 1 (1 -λ) + β 0 ,
where β j are continuous functions of θ k and a ii for i = 1, . . . , s. Hence, the eigenvalues become 1 -µ, where µ is a root of the parametrized polynomial

p(s) k (t) = t s + β s-2 t s-2 + β s-3 t s-3 + . . . + β 1 t + β 0 .
X GSU,u k =      1 0 . . . . . . . . . 0 ⋆ . . . ⋆ M k (P GSU,u k ) -1 2:s,2:s      , X GSL,l k =      1 ⋆ . . . . . . . . . ⋆ 0 . . . 0 (P GSL,l k ) -1 M k 2:s,2:s     
, and hence have one eigenvalue equal to one for each k. The entries replaced by ⋆ above do not affect the spectrum, only the eigenbasis. This manuscript is for review purposes only.

These results suggest 1 as a natural "central point" of the spectrum of the preconditioned system, generalizing the observations made for s = 2. We note that using these results we get both quantitative and qualitative insight into the spectra shown in [START_REF] Rana | A New Block Preconditioner for Implicit Runge-Kutta Methods for Parabolic PDE Problems[END_REF]], e.g., we see that for s = 3 the eigeninformation of M (P u ) -1 and (P l ) -1 M can still be obtained explicitly (see also [START_REF] Outrata | Schwarz methods, Schur complements, preconditioning and numerical linear algebra[END_REF]Section 7.4]) and on the other hand for s ≥ 6 there is no hope for these in general -but any bound on the eigeninformation of L can be used to obtain a bound on the eigeninformation of the preconditioned system by calculating with X k , see [START_REF] Gander | Spectral analysis of implicit 2-stage block Runge-Kutta preconditioners[END_REF]Section 4].

We show the spectra of the preconditioned systems and the corresponding GM-RES convergence behavior in Figure 2 and3, demonstrating observations and results from above. Notably, the bounds leave something to be desired, especially for P d where they are not descriptive at all. Moreover, increasing s seems to noticeably affect the quality of the preconditioners -see also [START_REF] Rana | A New Block Preconditioner for Implicit Runge-Kutta Methods for Parabolic PDE Problems[END_REF] for further numerical tests with various s and h. These numerical examples (as well as the ones in [START_REF] Clines | Efficient order-optimal preconditioners for implicit Runge-Kutta and Runge-Kutta-Nyström methods applicable to a large class of parabolic and hyperbolic PDEs[END_REF][START_REF] Gander | Spectral analysis of implicit 2-stage block Runge-Kutta preconditioners[END_REF]) are, as far as we can tell, representative of the general experience with these preconditioners. We highlight several key features illustrated in Figures 2 and3 that remained true in all of our experiments: 1. For s small, we have observed the staircase-like convergence behavior visible in the left upper-most plot in Figure 3 and this was most pronounced for the preconditioner P d .

2. We have usually not observed the desired superlinear convergence behavior, except for a speed-up after an initial stagnation (or slower speed convergence) phase.

3. In the vast majority of cases, the number of GMRES iterations to reach a certain tolerance grows only very moderately under mesh refinement and for P u , P l it remains almost constant. [START_REF] Driscoll | From potential theory to matrix iterations in six steps[END_REF]. In all of the experiments the spectra had the characteristic arc-like structure that we see in Figure 2.

Our goal is to explain all these features here as well as to investigate other bounds or estimates that would be more descriptive of the convergence behavior. This insight is of clear interest on its own but can be also used to further improve the used methods, e.g., looking at numerical optimization of the Butcher tableau in the spirit of [START_REF] Gander | Spectral analysis of implicit 2-stage block Runge-Kutta preconditioners[END_REF]Section 4]. We also note that the above results translate in a straight-forward fashion to the transformed system after we multiply (2.7) with (A -1 ⊗ I n ) from the left, obtaining

A -1 ⊗ I n - τ h 2 I s ⊗ L =:M transf k m = (A -1 ⊗ I n )    1 h 2 Lu m-1 + b (BC) (t m-1 + c i τ ) . . . 1 h 2 Lu m-1 + b (BC) (t m-1 + c i τ )    ,
and getting analogously the preconditioners,

R d = diag A -1 ⊗ I n - τ h 2 I s ⊗ L, R l = (D A -1 U A -1 ) ⊗ I n - τ h 2 I s ⊗ L and R u = (L A -1 D A -1 ) ⊗ I n - τ h 2 I s ⊗ L, R GSL = A -1 L ⊗ I n - τ h 2 I s ⊗ L and R GSU = A -1 U ⊗ I n - τ h 2 I s ⊗ L,
where A -1 has the LDU factorization

A -1 = L A -1 D A -1 U A -1 and A -1 L,U are
defined analogously to (2.9). These preconditioners were proposed in [START_REF] Neytcheva | Numerical Solution Methods for Implicit Runge-Kutta Methods of Arbitrarily High Order[END_REF] and then

This manuscript is for review purposes only. The spectra seemingly assemble in s "branches" in the first row and into s -1 "branches" in the other two with a central point at 1 + 0i. We set N = 50. used further in [START_REF] Munch | Stage-parallel fully implicit Runge-Kutta implementations with optimal multilevel preconditioners at the scaling limit[END_REF] but also [START_REF] Southworth | Fast solution of fully implicit Runge-Kutta and discontinuous Galerkin in time for numerical PDEs, Part I: The linear setting[END_REF][START_REF] Southworth | Fast solution of fully implicit Runge-Kutta and discontinuous Galerkin in time for numerical PDEs, Part II: nonlinearities and DAEs[END_REF]. For a general Butcher tableau, it is not possible to say whether the preconditioned transformed system gives a better performance than the original one. However, in [START_REF] Southworth | Fast solution of fully implicit Runge-Kutta and discontinuous Galerkin in time for numerical PDEs, Part I: The linear setting[END_REF][START_REF] Southworth | Fast solution of fully implicit Runge-Kutta and discontinuous Galerkin in time for numerical PDEs, Part II: nonlinearities and DAEs[END_REF] the authors propose different preconditioners and this analysis within this framework is going to be considered elsewhere. Also, we note that the extension of the above analysis for FEM discretization is a straightforward task -more details on both of these topics can be found in [START_REF] Outrata | Schwarz methods, Schur complements, preconditioning and numerical linear algebra[END_REF]Sections 7.6 and 7.7].

This manuscript is for review purposes only.

3.1. Spectral analysis. Next we turn to the spectral analysis, keeping in mind its limitation in the sense of Remark 2.1. For block-diagonal problems we obtain (3.11) ∥r ℓ ∥ ∥r 0 ∥ ≤ min

φ(0)=1 deg(φ)≤ℓ max j=1,...,n
∥φ (X j )∥ , which was studied in [START_REF] Faber | On Chebyshev polynomials of matrices[END_REF], where the authors showed that the extremal polynomials (i.e., the polynomial realizing the above bound) satisfies the equioscillation property but only every s iterations, where s is the size of the diagonal blocks. Relabeling the blocks in (3.11) we get

∥r ℓ ∥ ∥r 0 ∥ ≤ min φ(0)=1 deg(φ)≤ℓ max j=1,...,n ∥φ (X j )∥ = min φ(0)=1 deg(φ)≤ℓ max θj ∈sp( τ h 2 L) φ X θj .
Assuming each X θj is diagonalizable as in Proposition 3.3, we notice that {θ j } covers reasonably well the intervals I h,τ,... as h → 0 (see (3.6)) and, in the spirit of Section 2, the natural bound of (3.11) becomes

∥r ℓ ∥ ∥r 0 ∥ ≤ min φ(0)=1 deg(φ)≤ℓ max θ∈I h,τ,... ∥φ (X θ )∥ .
First, let us assume there is a uniform bound κ (S θ ) ≤ κ S for all θ ∈ I h,τ,... , which experimentally seems to be the case (see [START_REF] Outrata | Schwarz methods, Schur complements, preconditioning and numerical linear algebra[END_REF]) and can be confirmed analytically for s = 2, 3 (see [START_REF] Gander | Spectral analysis of implicit 2-stage block Runge-Kutta preconditioners[END_REF]) -this is an important and non-trivial assumption and a proper justi-

fication is an open problem. Next, we notice that the matrices X θ depend smoothly4 on θ and as a result so do their eigenproperties. In particular, the eigenvalues ξ (i) θ of X θ will -by definition -form an algebraic curve5 with s arcs (sometimes also called branches) some of which can be degenerate, e.g., reduced to just a point (incidentally, this is the case for at least one arc of the algebraic curve for any of the triangular preconditioners due to Corollary 3.5). Denoting the algebraic curve for the given Butcher tableau A and a choice of preconditioner P ⋆ by Γ, we obtain (3.12)

∥r ℓ ∥ ∥r 0 ∥ ≤ min φ(0)=1 deg(φ)≤ℓ max θ∈I h,τ,... κ (S θ ) max i=1,...,s φ ξ (i) θ ≤ κ S min φ(0)=1 deg(φ)≤ℓ max ξ∈Γ |φ (ξ)| .
Notice that if we replace in (3.12) the interval I h,τ,... with its limit I lim as h, τ → 0 (see (3.6)), we obtain a bound for all mesh sizes. Noticing that, in our case, the preconditioned system matrix has a limit as θ tends to either of the endpoints of I lim , it follows that the arcs of the corresponding algebraic curve correspond to the eigenvalues of these limit matrices. Hence, the effect of mesh refinement becomes sampling more points along Γ and stretching it towards these fixed endpoints (and possibly in increasing κ S ). This suggests that from a certain mesh size onward, the mesh refinement will have little effect on Γ and hence will not affect the min-max part of (3.12), shedding some light on why these preconditioners are quite robust under mesh refinement.

Remark 3.6. Note that the numerical experiments in [START_REF] Rana | A New Block Preconditioner for Implicit Runge-Kutta Methods for Parabolic PDE Problems[END_REF][START_REF] Clines | Efficient order-optimal preconditioners for implicit Runge-Kutta and Runge-Kutta-Nyström methods applicable to a large class of parabolic and hyperbolic PDEs[END_REF] as well as in [START_REF] Outrata | Schwarz methods, Schur complements, preconditioning and numerical linear algebra[END_REF] and in Section 4 clearly show that the spectra of the preconditioned systems cover reasonably well an algebraic curve. For two-stage methods, this behavior has been observed, proved and used to obtain descriptive GMRES bounds in [START_REF] Gander | Spectral analysis of implicit 2-stage block Runge-Kutta preconditioners[END_REF]. Moreover, for any algebraic curve Γ we have Γ = ∂ C Γ, which is convenient from the point of view of choosing σ non-discr , see Remark 2.1 and below.

We also emphasize that, in general, these preconditioners do not cluster eigenvalues (that is, any more than the θ ∈ I h,τ,... already are) but rather place them along a particular algebraic curve Γ ⊂ C. Hence, in general, we can reasonably expect linear convergence as opposed to superlinear, which can often be linked with clusters and numbers of outliers, in the sense of [16, Section 5.6.4].

Remark 3.6 also explains that the bound (2.11) is unlikely to be very descriptive or even usable. Indeed, the algebraic curves can reach into the left half-plane {Re(z) < 0}

(making the bound useless due to 0 being included in the bounding circle) or, in the more favorable case, the arcs of the algebraic curve are extremely unlikely to align with the circle so that the bound have some resemblance of being what we earlier called appropriate. Naturally, the bound on the right-hand side of (3.12) is constructed to remedy that but the key question becomes if this bound is also functional, namely if we can (approximately) evaluate it.

To this end, we follow the excellent paper [START_REF] Driscoll | From potential theory to matrix iterations in six steps[END_REF] on this topic and start by looking at the asymptotic convergence rate (justified by Remark 3.6 above). Considering (3.12)

we are led to look at the so-called logarithmic capacity of Γ, denoted by cap(Γ), which can be viewed as a measure of a compact set without isolated points in C.

In fact it is known to asymptotically correspond to the maximal modulus of the where the quantity on the left-hand side relates to the quantities we have seen in the GMRES bounds. There are two important caveats to using cap(Γ). The first one, which has been also highlighted as a caveat for using the analysis in [START_REF] Driscoll | From potential theory to matrix iterations in six steps[END_REF] overall, is the fact that that (3.13) only provides some information about the limit behavior as ℓ → +∞, whereas we are interested in the behavior for relatively small values of ℓ, say ℓ ≤ 50 or 100. To large extend this issue is addressed by Remark 3.6 that states that we expect a linear convergence rate throughout the iteration. The second one is the fact that (3.13) describes the limit scaling of the maximal modulus over all polynomials -it lacks the crucial scaling φ(0) = 1 of Krylov methods. This issue can be fixed by re-scaling (see [4, Section 2]), shifting our attention from the logarithmic capacity to Green's functions associated with Γ, as long as Γ is compact and without any isolated points.

Things simplify considerably if we assume that Γ is connected as then the normalized quantity

   min φ(0)=1 deg(φ)≤ℓ max z∈Γ |φ(z)|    1/ℓ
can be evaluated directly using conformal maps, in particular the Schwarz-Christoffel

This manuscript is for review purposes only.

maps. Without going into the details (the interested reader can find these in [4, Sections 2 and 3]), we obtain the asymptotic convergence factor estimate ρ est as

ρ est := lim ℓ→+∞    min φ(0)=1 deg(φ)≤ℓ max z∈Γ |φ(z)|    1/ℓ = 1 |Φ(0)| ,
where Φ(z) is the Schwarz-Chriostoffel map that maps the exterior of Γ to the exterior of the unit circle. In [4, Section 3, Theorem 2 and below], the authors put this as ". . . if Γ is connected, the estimated asymptotic convergence factor for a matrix iteration depends on how far the origin is from Γprovided that this distance is measured by level curves associated with the exterior conformal map."

We would like to emphasize the word estimate when talking about ρ est because we truly do not get a bound anymore -in fact we get an underestimate as highlighted also in [4, Section 5, equation (STEP1) and also Table 1]. However, we expect this estimate to be descriptive as explained above.

For not too complicated connected, compact sets the map Φ and its value at the origin can be calculated using the Schwarz-Christoffel MATLAB toolbox [START_REF] Driscoll | A MATLAB toolbox for Schwarz-Christoffel mapping[END_REF], but we immediately notice that in Figure 2 the set of eigenvalues along Γ is not connected and the actual algebraic curve Γ itself is also not available in an easy form, i.e., neither of these can be directly given as an input to the SC toolbox. We take the natural next step and approximate Γ by its linear interpolation based on the available eigenvalues

ξ (i)
θ . The linear interpolation gives us a good approximation of the arcs of Γ and we use the point 1 + 0i as the natural point to join them (also by linear interpolation) and denote the resulting set Γ h . Recalling the limit behavior in (3.6), we also see that Γ h will tend towards Γ as h → 0 for our model problem.

The calculation of ξ (i)

θ is independent for each k = 1, . . . , n but for large n the SC toolbox can suffer numerically when calculating with Γ h that is densely populated by the interpolation points -both in the sense of large computational complexity as well as in the sense of numerical issues (called over-crowding, see [START_REF] Driscoll | A MATLAB toolbox for Schwarz-Christoffel mapping[END_REF] but also [START_REF] Driscoll | Schwarz-Christoffel mapping[END_REF]Section 2.6]). Moreover, we usually have only rough estimates on the extremal eigenvalues θ min and θ max of L rather than its full spectrum. To this end, we recall the idea in [9,

Section 4] and instead of calculating Γ h we use the information about θ min,max and artificially sample a fixed number of "fake" points ϑ k between them, say q of them.

Then we replace θ k by ϑ k in the definition of Γ h , obtaining Γ q -an approximation of Γ h (and a further approximation of Γ) based on the linear interpolation given by the eigenvalues of the matrices X ϑ k . We illustrate these points in Figure 4.

Another key point is that using the SC toolbox6 -namely the functions extermap and evalinv -has difficulties (as far as we understand it) when the arcs of Γ q intersect, e.g., as is the case for s = 8 and the preconditioner P l , see Figure 2. Intuitively, this makes sense as the exterior of Γ q then has multiple components, making the original set-up more complicated (a theoretical treatment of such problems could be approached based on [START_REF] Embree | Green's functions for multiply connected domains via conformal mapping[END_REF]). We address this issue by taking the "envelope" of the arcs -if two arcs intersect, we follow the one staying outwards, e.g., in the case of s = 6

and the preconditioner P l we would exclude some portion of the arcs with smaller imaginary part (the densely populated portions) as these lie "inward" relative to the arcs with the larger imaginary part, see Figure 5. Finally, we illustrate the calculated Schwarz-Christoffel maps -or rather their contours -in Figure 5 together with the used inputs Γ q (with the exception of s = 6, 8 and the preconditioner P l , where we used the "envelopes") and also the asymptotic convergence factor estimate ρ est in Figure 6. First, we see that the results in Figure 6 fully support the arguments in Remark 3.6 for considering ρ est as the descriptive quantity for the convergence factor.

Including an estimate for κ S then gives also an estimate for GMRES convergence -not just its rate, see Section 4. Second, we note that for s = 8 and the preconditioner P u , the arcs turned so that the right-most arcs almost intersect themselves. This causes problems for the toolbox, which during the calculations raises a flag stating that the calculated map did not converge as expected. Although the predicted ρ est seems accurate, we see in Figure 5 that contours have ripples, confirming that the calculated results should be taken with caution. This can be fixed by a similar "envelope-like" approach we described for s = 6, 8 and the preconditioner P l , see Section 4, obtaining a further approximation. Although there are a few similar caveats concerning the implementation of the above ideas, we have always found that a simple solution (such as considering the envelope or pruning the fake points in order to alleviate the crowding) can be used to fix them and still give an appropriate insight into the GMRES convergence factor. As long as κ S does not completely dominate the ideal GMRES bound (2.10) this then translates to descriptive GMRES convergence estimates, see Section 4.

The above analysis also gives insight into the staircase-like behavior, which has been observed and explained for s = 2 and the preconditioner P d in [START_REF] Gander | Spectral analysis of implicit 2-stage block Runge-Kutta preconditioners[END_REF] working with the minimal residual polynomial φ MR ℓ (sometimes also called the GMRES polynomial; see [START_REF] Liesen | Krylov Subspace Methods: Principles and Analysis[END_REF]Section 5.7.1]). The arguments used in [START_REF] Gander | Spectral analysis of implicit 2-stage block Runge-Kutta preconditioners[END_REF] remain valid as long as the branches are not very close to each other 7 -as long as the branches are far apart, the maximum Fig. 5. In red: the curves Γq (first plots 1 to 7) and their "envelopes" (plots 8 and 9) for the Gauss Butcher tableau, taking q = 15. In black: the contours of the corresponding Schwarz-Christoffel map of the exterior of these curves (or envelopes) mapped to the exterior of the unit circle, see extermap in [START_REF] Driscoll | A MATLAB toolbox for Schwarz-Christoffel mapping[END_REF].

of the polynomial φ MR ℓ will decrease significantly more at the steps ℓ = s • j for j = 1, 2, . . . because only then each branch can get some attention. If the branches become close, then we do not expect this extra jump because keeping the absolute value of the polynomial small along one of the branches naturally translates into keeping the absolute value of the polynomial also small enough along another one. This is most pronounced in the first s iterations of GMRES, as we can see in Figure 6, where the convergence curves begin with a slower convergence phase -precisely s steps -for P d and P u , in contrast to the ones of P l , where the arcs intersect and are, in general, closer to each other. We illustrate this further in Figure 7 for the preconditioner P d for s = 4, 8 by looking at the polynomial φ MR ℓ and its roots (called harmonic Ritz values). We see that in the first row (4 branches, far apart) the possibility of "placing" one root along each of the branches was much more crucial (resulted in a more significant decrease of the modulus of the polynomial over the spectrum of the preconditioned system) than for the second row (8 branches with two complex conjugate pairs of branches that are close to each other). We note that an example of explanation (and prediction) of a complete staircase behavior of GMRES can be found in [4, Figure 9 and below].

other".

This manuscript is for review purposes only. GMRES convergence GMRES convergence estimate Fig. 6. The convergence behavior of preconditioned GMRES, using the Gauss Butcher tableau, together with the convergence estimate based on the calculated asymptotic convergence factor estimate ρest. Fig. 7. The level curves of the GMRES polynomials φ M R ℓ for the preconditioned system P l -1 M together with the spectrum of this system as well as the roots of φ M R ℓ (so-called harmonic Ritz values). We set N = 50.

We also want to comment on a similarity with the results in [START_REF] Liesen | Convergence of GMRES for tridiagonal Toeplitz matrices[END_REF][START_REF] Liesen | GMRES convergence analysis for a convection-diffusion model problem[END_REF]. There, the authors addressed the question of delay of convergence by using similar formulations to ours, also obtaining a GMRES problem reformulated as for a block-diagonal matrix using Kronecker-product-like techniques as in Lemma 3.1. In particular, in [15, Section 3.1] the authors use the equality

∥r ℓ ∥ = min φ(0)=1 deg(φ)≤ℓ φ       X 1 . . . X n       r 0 = min φ(0)=1 deg(φ)≤ℓ n j=1 φ (X j ) s (i) 0 ,
This manuscript is for review purposes only.

where s (i) 0 is the i-th subvector of length s of Q T Πr 0 , to obtain a lower bound

(3.14) ∥r ℓ ∥ 2 = min φ(0)=1 deg(φ)≤ℓ n j=1 φ (X j ) s (i) 0 2 ≥ n j=1 min φ(0)=1 deg(φ)≤ℓ φ (X j ) s (i) 0 2
on the GMRES convergence behavior, explaining the initial stagnation phase in an advection-diffusion problem. This way they bound the global minimization problem (corresponding to solving a problem with the block-diagonal matrix diag(X 1 , . . . , X n ))

by the sum of the local minimization problems (each given by the small s-by-s matrix X j ). By careful analysis of the interplay of the right-hand side (or initial residual)

and the diagonal blocks in [15, Section 3.1] (there the diagonal blocks are, moreover, tridiagonal and Toeplitz), the authors conclude ". . . the presence of at least one system with tridiagonal Toeplitz matrix T j = tridiag(γ j , λ j , µ j ) that is 'close to the Jordan block' (cf. [START_REF] Liesen | GMRES convergence analysis for a convection-diffusion model problem[END_REF]Section 3.3] but see also [START_REF] Liesen | Convergence of GMRES for tridiagonal Toeplitz matrices[END_REF]), and with l representing the index of the first significant entry of the corresponding right-hand side, prevents fast convergence of GMRES for the first N -l steps (N being the size of the blocks T j ) . . . . . . As explained in Section 3.1, the lower bound is useless for analyzing the convergence behavior after the step N -l, possibly even earlier. Hence the above approach cannot be used for quantifying any possible acceleration of convergence after the initial phase. "

We see that the approach is fundamentally different -both in the intended direction as well as in the results it can deliver -in spite of the fact that it works with the same technique.

We finalize this section with a remark on the field of values (sometimes also called the numerical range) and pseudospectra, which sometimes are extremely useful to understand and predict GMRES convergence behavior, especially if the eigenbasis of the system matrix is ill-conditioned, see, e.g., [START_REF] Embree | How descriptive are GMRES convergence bounds?[END_REF] and also [START_REF] Liesen | Krylov Subspace Methods: Principles and Analysis[END_REF]Section 5.7.3,pp. 296] and the references therein. where the X k are given as in (3.7) and the set addition is understood element-wise, i.e., ν(X

1 ) + ν(X 2 ) = {α 1 + α 2 | α 1 ∈ ν(X 1 ), α 2 ∈ ν(X 2 )}, or, more generally ν M P -1 ⊂ κ(Q) n i=1 ν (X k )
and analogously for ν P -1 M .

For the pseudospectrum we obtain an analogous formula, namely

σ δ M P -1 ⊂ κ(Q) n i=1 σ δ (X k )
and analogously for ν P -1 M .

This manuscript is for review purposes only.

In other words, the principle of working with the small matrices X k instead of the large matrix M P -1 naturally applies also to the other standard techniques for analyzing GMRES convergence behavior. However, adapting and using bounds based on field of values or the pseudospectrum of the preconditioned system for this set-up remains a topic for future research.

4. Numerical Examples. In this section we use the above analysis for more involved settings and, more importantly, also demonstrate the convergence estimates (instead of only the convergence rate estimates). We use a FEM discretization in space 8 for different geometries in Example 1 and 2, see Figure 8. We also fix the number of time steps to balance the spatial and time error (see the (L2) definition in Section 2), namely we fix

τ = h 2 p ,
where the 2 in the numerator is the order of the spatial error (since we use linear Lagrange polynomials in the FEM discretization) and p is the order of convergence of the Runge-Kutta method. We show for both examples the GMRES convergence together with the convergence estimates, namely

∥r ℓ ∥ ∥r 0 ∥ ⪅ min κ est S ρ ℓ est , 1 ,
where the estimate κ est S of κ S is computed from the eigenbasis condition numbers of the "fake sampled" matrices X ϑ k for k = 1, . . . , q. In our experience, the best results are obtained with q ≈ 15 -20, as increasing q further leads to crowding problems in the SC toolbox and eventually to problems with the convergence of the Schwarz-Christoffel map. We also found that spacing the fake points ϑ k logarithmically in the corresponding interval somewhat alleviates this issue and leads to more accurate predictions of the arcs of the given algebraic curve. We also recall that the seeming independence of the preconditioner quality on the spatial mesh size h was sufficiently documented elsewhere (see [START_REF] Rana | A New Block Preconditioner for Implicit Runge-Kutta Methods for Parabolic PDE Problems[END_REF][START_REF] Neytcheva | Numerical Solution Methods for Implicit Runge-Kutta Methods of Arbitrarily High Order[END_REF][START_REF] Clines | Efficient order-optimal preconditioners for implicit Runge-Kutta and Runge-Kutta-Nyström methods applicable to a large class of parabolic and hyperbolic PDEs[END_REF][START_REF] Gander | Spectral analysis of implicit 2-stage block Runge-Kutta preconditioners[END_REF][START_REF] Outrata | Schwarz methods, Schur complements, preconditioning and numerical linear algebra[END_REF]) and explained in Section 3 so that in our eyes, there is no need to address this direction here. Illustration of the solutions as well as further numerical experiments can be found in [START_REF] Outrata | Schwarz methods, Schur complements, preconditioning and numerical linear algebra[END_REF]Chapter 7].

Last but not least, we have not set a relative residual tolerance criterion for stopping GMRES, meaning that GMRES went on until either the relative residual was on the level of machine precision or the maximum number of iterations was reached. This is not a good choice from the point of view of the solution process efficiency but since our primary focus is on studying the preconditioners, we found this reasonable.

Example 1: Cookies in the oven . The first problem is a simulation of baking cookies in an electrical oven projected in 2D, an idea borrowed from [START_REF] Kressner | Low-rank tensor Krylov subspace methods for parametrized linear systems[END_REF]. The cookies have a worse heat conductivity than the surrounding air (piecewise constant in space and constant in time) and the setting demands various boundary conditions, resulting

∂u ∂t u = div (σ∇u) + f in Ω × (0, T ], ∂u ∂n u = 0 on Γ N × (0, T ], ∂u ∂n u + pu = 0 on Γ R × (0, T ], u = 0 at Ω × {0},
with Ω = (0, 4) × (0, 4) and the boundary of Ω is split into the Neumann and Robin parts Γ N , Γ R . We set the data as

Γ N = {x = 0} ∪ {y = 0} ∪ {y = 4}, Γ R = {x = 4}, p = 1, σ = 10 3 if (x, y) ∈ Cookie, 1 otherwise, f (x, y, t) = 3 if ∥(x, y) -(2, 2)∥ ≤ 1, 0 otherwise,
and show the GMRES convergence behavior with the estimates in Figure 9 as well as the sampling of the algebraic curves in Figure 10. Summary. Overall, we observe that the convergence factor estimates delivered very accurate results even for these more involved problems but the conditioning of the eigenbasis of the matrices X ϑ k notably deteriorated as we increased s. The fact that this does not show up in the GMRES convergence behavior suggests that more delicate bounds, such as mentioned in Remark 3.7 could give a more detailed insight

This manuscript is for review purposes only. Envelope/Polygon (i) (q = 15) Fig. 12. The algebraic curve polygon approximations that are used in the Schwarz-Christoffel MATLAB toolbox to calculate ρest for Example 2 -for some settings these correspond to the eigenvalues ξ (i) ϑ and in some these only enclose these. into the matter. However, in all cases the bounds lag behind the actual convergence behavior by 10-20 iterations, which is often still considered to be reasonably accurate.

We also showed the polygons used in the Schwarz-Christoffel toolbox -notice that in many of the plots we excluded part of the arcs, mainly because either (a) the arcs intersected and we took the envelope of the algebraic curve (usually for the preconditioner P l ) or (b) the points sampled along the arcs crowded sections of the arcs, which caused issues for the toolbox. In such cases we sparsified these regions by dropping some of these points. As a result, the Schwarz-Christoffel external map converged better and faster than for the problem in Section 3.1 and the contours were "ripple-free" for all of our problems, otherwise looking almost precisely as the ones in This manuscript is for review purposes only. 5. Concluding remarks. Our main goal has been to understand the block preconditioners considered in [START_REF] Rana | A New Block Preconditioner for Implicit Runge-Kutta Methods for Parabolic PDE Problems[END_REF][START_REF] Clines | Efficient order-optimal preconditioners for implicit Runge-Kutta and Runge-Kutta-Nyström methods applicable to a large class of parabolic and hyperbolic PDEs[END_REF][START_REF] Neytcheva | Numerical Solution Methods for Implicit Runge-Kutta Methods of Arbitrarily High Order[END_REF] in more detail and to try to explain their success and/or limitations. This goal was, in our eyes, mostly achieved but could be further improved in the sense of Remark 3.7 or by considering a more refined version of the bound (2.10), see [6, Section 2.1, equations (2.1) and (EV')] -this remains an area of interest for us for the future. Moreover, the above analysis can be directly used to try to optimize the Runge-Kutta methods, following the ideas in [START_REF] Rana | A New Block Preconditioner for Implicit Runge-Kutta Methods for Parabolic PDE Problems[END_REF][START_REF] Outrata | Schwarz methods, Schur complements, preconditioning and numerical linear algebra[END_REF][START_REF] Gander | Spectral analysis of implicit 2-stage block Runge-Kutta preconditioners[END_REF]. We also note that in practice, solving with either of the matrices P d,u,l,GSU,GSL,... is often done with some level of inaccuracy, e.g., using a multigrid method. The question of interaction of this inaccuracy with the overall GMRES convergence is an important one and to the best of our knowledge has been addressed only numerically in [START_REF] Outrata | Schwarz methods, Schur complements, preconditioning and numerical linear algebra[END_REF]Chapter 7]. We also note that adapting the above analysis to the framework presented in [START_REF] Southworth | Fast solution of fully implicit Runge-Kutta and discontinuous Galerkin in time for numerical PDEs, Part I: The linear setting[END_REF][START_REF] Southworth | Fast solution of fully implicit Runge-Kutta and discontinuous Galerkin in time for numerical PDEs, Part II: nonlinearities and DAEs[END_REF], or reformulating it from the vector equation to the matrix equation as suggested in [START_REF] Palitta | Optimality properties of Galerkin and Petrov-Galerkin methods for linear matrix equations[END_REF], and to study in detail the comparison of these approaches for the IRK setting are attractive directions for future research.

Fig. 1 .

 1 Fig. 1. Left: grid points for N + 1 = 4; right: lexicographical ordering of the unknowns for N + 1 = 4.

  , for i = N + 1, . . . , N (N -1) -1, where b

2

 2 collects the known values from the source terms, given by g and f , at the given point. Combining the unknowns in each grid column into one vector denoted by u k (t), i.e.,u k (t) := u N k+2 u N k+3 • • • u N (k+1)-1 T (t), u(t) := u T 1 (t) • • • u T N -1 (t) T ,and also analogously for b k (t) and b(t), we rewrite (2.Lu(t) + b (ST) (t)

Corollary 3 . 5 ([ 19 ,

 3519 Proposition 7.6]). Under the assumptions of Proposition 3.3, the block upper-triangular preconditioners P GSU,u give(3.10)

Fig. 2 .

 2 Fig. 2. The spectra of the preconditioned systems M P u,d -1 and P l -1 M for s = 4, 6, 8 and for three classical choices of fully implicit Runge-Kutta schemes -Gauss, RadauIIA and LobattoIIIC.The spectra seemingly assemble in s "branches" in the first row and into s -1 "branches" in the other two with a central point at 1 + 0i. We set N = 50.

Fig. 3 .

 3 Fig.3. The preconditioned GMRES convergence behavior for the preconditioned systems M P u,d -1 and P l -1 M for s = 4, 6, 8 and three classical choices of fully implicit Runge-Kutta schemes -Gauss, RadauIIA and LobattoIIIC -together with the GMRES bound (2.11) with c = 1 (we set the values to 1 if ρ ≥ 1). We set N = 50.

→

  cap(Γ), as ℓ → +∞,

Fig. 4 .

 4 Fig.4. The eigenvalues of the matrices X θ k (red) and X ϑ k (blue, for different values of q), using the preconditioner P d . Joining these together with line segments would yield the curves Γ h (red) and Γq (blue).

Remark 3 . 7 .

 37 Another commonly used bound for GMRES uses the field of values ν(C) or the δ-pseudospectrum σ δ (C) of the system matrix C. By a direct calculation we obtain, for our model problem, the field of values asν M P -1 = n i=1 ν (X k )and analogously for ν P -1 M ,

Fig. 8 .

 8 Fig. 8. The initial triangulations for Example 2 and 3 together with the boundary condition types and, for Example 2, also with highlighting the points with lower heat conductivity.

  in

Example 2 :

 2 The cabin heating . The second problem uses the 2D projection of an attic room of a cabin in the western Bohemia region, whose primary heating is the chimney (bottom-right corner, modeled with a Dirichlet boundary condition changing in time), with two windows (top and bottom) and a door (right), modeled with Robin boundary conditions with Robin parameters p w and p d , and a good insulation otherwise, modeled with a Neumann condition. We obtain the problem ∂u ∂t u = div (σ∇u) in Ω × (0, T ], ∂u ∂n u = 0 on Γ N × (0, T ], ∂u ∂n u + pu = 0 on Γ R × (0, T ], u = 0 at Ω × {0},This manuscript is for review purposes only.

Fig. 9 .Fig. 10 .

 910 Fig. 9. The GMRES convergence behavior with the convergence estimates based on ρest for Example 1 with n = 26985.

Fig. 11 .

 11 Fig. 11. The GMRES convergence behavior with the convergence estimates based on ρest for Example 2 with n = 26985.

Figure 5 .

 5 Figure 5.

This manuscript is for review purposes only.

In the sense of the minimum is over a non-discrete set while the maximum is over a discrete one.

We denote the boundary of a set S ⊂ C in C by ∂ C S.

Intuitively, we could expect that the bound will be appropriate only if σ discr covers the entirety of σ non-discr but because polynomials of complex variables are harmonic we can conclude that the maximum of the modulus of a polynomial over the set σ non-discr is attained along ∂ C σ non-discr and therefore only the relation of ∂ C σ non-discr and σ discr is important for the GMRES bound, see[START_REF] Driscoll | From potential theory to matrix iterations in six steps[END_REF] Section 2].This manuscript is for review purposes only.

That is, for our model problem. This changes if we consider, e.g., an indefinite spatial operator L instead of the negative-definite Laplacian.

We say that Γ is an algebraic curve provided there exists a bi-variate polynomial p(θ, t) such that Γ = {(θ, ξ) | p(θ, ξ) = 0}. Locally, this can also be viewed through the lens of perturbation theory, see[12, Chapter 2 Section 1.1].This manuscript is for review purposes only.

In our case, Γq qualifies as a degenerate polygon acceptable by the toolbox.This manuscript is for review purposes only.

In[START_REF] Gander | Spectral analysis of implicit 2-stage block Runge-Kutta preconditioners[END_REF], the branches are two line segments parallel to the imaginary axis that are, moreover, reasonably well separated along the real line, i.e., a natural case of being "not very close to each This manuscript is for review purposes only.

Wherever we talk about a FEM discretization, we use linear Lagrange polynomials on conforming triangular meshes. Those are refined by the standard quadrisection of a triangle, with additional post-smoothing of the mesh.This manuscript is for review purposes only.

Acknowledgements. Some of the ideas were stimulated by conversations with Mark Embree, Patrick Farell, Miroslav Tůma and Petr Tichý and we would like to thank them for their inspiring comments and suggestions.

Submitted to the editors DATE. Funding: This work was partially supported by the SNF grant number 178752 and by the FCS Swiss Excellence PhD Fellowship program of the Swiss Federation (ESKAS No. 2019.0384).


Keywords: force-velocity, intensity-time, critical power, modelling, locomotion

Application of a force-velocity-endurance model to cycling, rowing and running locomotion. ⋆ B. Morel Introduction. The production of force during brief maximal effort is limited by a force-velocity relationship, which is characterized by a negative linear function and determined by the theoretical maximum force (F 0 ) and velocity (V 0 ) [START_REF] Dorel | Torque and power-velocity relationships in cycling: relevance to track sprint performance in world-class cyclists[END_REF]. For longer efforts, the intensity-duration relationship has been described mathematically by an asymptotic decay function, with the critical intensity being the asymptote [START_REF] Poole | Critical power: an important fatigue threshold in exercise physiology[END_REF]. The force and velocity produced during an exercise can be modified by changing the cadence during cycling, for example, which alters the intensity-duration relationship [START_REF] Wright | The 3-minute all-out cycling test is sensitive to changes in cadence using the lode excalibur sport ergometer[END_REF]. As the effort is prolonged, changes in force and velocity capacities may occur differently. To account for the influence of the force-velocity and intensity relationships on each other, we developed a three-dimensional model that describes the force production capacity as a function of velocity and time, F (v, t). In addition, we designed a test to determine the F (v, t) relationship for an individual based on a 3-min all-out interspersed with a stop-start sprint (IFLET test) [START_REF] Morel | Force-velocity-endurance relationship: development of a mathematical model and validation of the iflet test 28[END_REF]. The aim of this study was to evaluate the parameters of the F (v, t) relationship for various types of human locomotion (cycling, rowing and running) and to compare the parameters obtained for different populations based on their training status.

Methods. The IFLET test was administered to 49 participants across cycling, rowing, and running. During cycling, 21 moderately active individuals who were not cyclists (NC), 19 subelite trained cyclists (SC), and nine elite cyclists (EC) participated. Twelve elite rowers participated for the rowing task and 16 U21 elite rugby players participated for running locomotion. The test involved a 3-minute all-out effort, with conditions changing every 30s to assess the force-velocity relationship at a specific point in time. During cycling, this was achieved by suddenly blocking the flywheel, while in rowing, a motor replaced the flywheel and was controlled to produce varying force-velocity conditions. During running, the test was conducted in shuttle mode, which required participants to start sprinting at zero speed. The force-velocity-time data recorded or computed during the various locomotions was fitted to the F (v, t) model's parameters to obtain the initial force (F 0i ), velocity (V 0i ), critical force (F 0c ), and velocity (V 0c ) capacities, as well as the time constant (τ ).

Results and discussion. The goodness of fit of the model from experimental data was excellent for all locomotion (all r2 > 0.93). Considering the group effect for cycling task locomotion, no difference was observed for V 0c (NC : 58.9 ± 12.3 %; SC : 59.3 ± 9.7 %; EC: 55.7 ± 6.3 %). However, all groups were statistically different in terms of F 0c (NC : 51.4 ± 11.2 %; SC : 64.2 ± 7.3 %; EC: 71.6 ± 10.4 %).

Conclusions and perspectives. This is the first time a model has been developed that simultaneously considers both velocity and time to describe the force capacity. This model accurately fits the experimental data obtained from the IFLET test, which is a 3-minute all-out sprint exercise interspersed with different locomotion types, such as running, cycling, or rowing. The applications of this approach are numerous and can be used in the evaluation of physical capacities as well as for performance enhancement through training or optimization of human-material interactions.

  1, * , M. Bowen 1 , S. Dorel 2 , Y. Bertron 1 , C. Delhaye 1 , L. Blervaque 1 , P-F. Busseneau 1 , L. Messonnier 1 , and P. Samozino 1

Interuniversity laboratory in Biology of Movement, Savoie-Mont Banc University, Chambéry, France

Laboratory Motricité, Interactions, Performance, Nantes University, France

⋆ This work was supported by French National Research Agency [ANR-22-CE14-0073]