Conditioning of implicit Runge-Kutta integration for finite element approximation of linear diffusion equations on anisotropic meshes

TL;DR

This paper analyzes how the conditioning of implicit Runge-Kutta methods for finite element diffusion problems is affected by mesh anisotropy and proposes bounds that incorporate mesh and diffusion properties, with effective preconditioning strategies.

Contribution

It provides explicit bounds on the condition number for implicit RK methods on anisotropic meshes, considering different solution strategies and preconditioning effects, with geometric interpretations.

Findings

Condition number bounds depend on mesh size, nonuniformity, and diffusion matrix.

Diagonal preconditioning effectively mitigates mesh nonuniformity effects.

Numerical examples validate the theoretical bounds and strategies.

Abstract

The conditioning of implicit Runge-Kutta (RK) integration for linear finite element approximation of diffusion equations on general anisotropic meshes is investigated. Bounds are established for the condition number of the resulting linear system with and without diagonal preconditioning for the implicit Euler and general implicit RK methods. Two solution strategies are considered for the linear system resulting from general implicit RK integration: the simultaneous solution (the system is solved as a whole) and a successive solution which follows the commonly used implementation of implicit RK methods to first transform the system into smaller systems using the Jordan normal form of the RK matrix and then solve them successively. For the simultaneous solution in case of a positive semidefinite symmetric part of the RK coefficient matrix and for the successive solution it is shown that . If the smallest eigenvalue of the symmetric part of the RK coefficient matrix is negative and the simultaneous solution strategy is used, an upper bound on the time step is given so that the system matrix is positive definite. The obtained bounds for the condition number have explicit geometric interpretations and take the interplay between the diffusion matrix and the mesh geometry into full consideration. They show that there are three mesh-dependent factors that can affect the conditioning: the number of elements, the mesh nonuniformity measured in the Euclidean metric, and the mesh nonuniformity with respect to the inverse of the diffusion matrix. They also reveal that the preconditioning using the diagonal of the system matrix, the mass matrix, or the lumped mass matrix can effectively eliminate the effects of the mesh nonuniformity measured in the Euclidean metric. Numerical examples are given.

Tables

Table 1. Table 1 : Eigenvalues and singular values of the implicit RK coefficient matrix Γ Γ \Gamma and the minimum eigenvalue of ( Γ + Γ T ) / 2 Γ superscript Γ 𝑇 2 (\Gamma+\Gamma^{T})/2 .

Method	Order	Eigenvalues	${\sigma }_{\max }$	${\sigma }_{\min }$	${\lambda }_{\min }\left(\left(\Gamma +{\Gamma }^T\right)/2\right)$
Gauss	4	$0.25\pm 0.1443i$	0.6319	0.1319	0
Gauss	6	0.2153, $0.1423\pm 0.1358i$	0.6629	0.0635	-0.0563
Radau IA	3	$0.3333\pm 0.2357i$	0.5	0.3333	0.25
Radau IIA	5	0.2749, $0.1626\pm 0.1849i$	0.8023	0.0923	-0.0822
Lobatto IIIA	4	0, $0.25\pm 0.1443i$	0.7947	0	-0.0736

Equations

$$\{\begin{array}{l}\begin{array}{rlrlrl}& {\partial }_{t}u=\nabla \cdot \left(\mathbb{D}\nabla u\right),& & \mathbf{x}\in \Omega ,& & t>0,\\ & u(\mathbf{x},t)=0,& & \mathbf{x}\in \partial \Omega ,& & t>0,\\ & u(\mathbf{x},0)=u_0(\mathbf{x}),& & \mathbf{x}\in \Omega ,& & \end{array}\end{array}$$

$$\exists d_{\min },d_{\max }>0:d_{\min }I\le \mathbb{D}(\mathbf{x})\le d_{\max }I,\forall \mathbf{x}\in \Omega ,$$

$${\int }_{\Omega }\frac{\partial u_h}{\partial t}{v}_h\mathrm{d}\mathbf{x}=-{\int }_{\Omega }{(\nabla v_h)}^T\mathbb{D}\nabla u_h\mathrm{d}\mathbf{x},\forall v_h\in V^h,t>0,$$

$${\int }_{\Omega }{u}_h(\mathbf{x},0)v_h\mathrm{d}\mathbf{x}={\int }_{\Omega }{u}^0(\mathbf{x})v_h\mathrm{d}\mathbf{x},\forall v_h\in V^h.$$

$$u_h(\mathbf{x},t)=\sum _{j=1}^{N_{vi}}{u}_j(t){\varphi }_j(\mathbf{x}),$$

$$M\frac{d\mathbf{u}}{dt}=-A\mathbf{u},$$

$$M_{kj}={\int }_{\Omega }{\varphi }_k{\varphi }_j\mathrm{d}\mathbf{x}\text{and}{A}_{kj}={\int }_{\Omega }\nabla {\varphi }_k\cdot \mathbb{D}\nabla {\varphi }_j\mathrm{d}\mathbf{x},k,j=1,\dots ,N_{vi}.$$

$$M{\mathbf{v}}_k+\Delta tA\sum _{j=1}^s{\gamma }_{kj}{\mathbf{v}}_j=-A{\mathbf{u}}^n,k=1,\dots ,s,$$

$${\mathbf{u}}^{n+1}={\mathbf{u}}^n+\Delta t\sum _{k=1}^s{b}_k{\mathbf{v}}_k,$$

$$\mathbb{A}\equiv I_s\otimes M+\Delta t\Gamma \otimes A.$$

$$(M+\Delta tA){\mathbf{u}}^{n+1}=M{\mathbf{u}}^n.$$

$$P^{-\frac{1}{2}}(M+\Delta tA)P^{-\frac{1}{2}}{\mathbf{w}}^{n+1}=P^{-\frac{1}{2}}M{P}^{-\frac{1}{2}}{\mathbf{w}}^n,{\mathbf{w}}^n=P^{\frac{1}{2}}{\mathbf{u}}^n.$$

$$\frac{1}{2}{M}_D\le M\le \left(1+\frac{d}{2}\right)M_D\text{and}{M}_{jj}=\frac{2\left|{\omega }_j\right|}{(d+1)(d+2)}.$$

$$\frac{1}{d+2}{M}_{lump}\le M\le \frac{d+2}{2}{M}_{lump}.$$

$$A\le (d+1)A_D.$$

$$C_{\hat{\nabla }}\sum _{K\in {\omega }_j}\left|K\right|{\lambda }_{\min }\left({(F_K^{\prime })}^{-1}{\mathbb{D}}_K{(F_K^{\prime })}^{-T}\right)\le A_{jj}\le C_{\hat{\nabla }}\sum _{K\in {\omega }_j}\left|K\right|{\lambda }_{\max }\left({(F_K^{\prime })}^{-1}{\mathbb{D}}_K{(F_K^{\prime })}^{-T}\right),$$

$$h_{{\mathbb{D}}^{-1}}={\left(\frac{{\left|\Omega \right|}_{{\mathbb{D}}^{-1}}}{N}\right)}^{\frac{1}{d}},{\left|\Omega \right|}_{{\mathbb{D}}^{-1}}=\sum _{K\in {\mathcal{T}}_h}\left|K\right|\sqrt{\det ({\mathbb{D}}_K^{-1})}.$$

$$\hat{h}=\sqrt{2}{\left(\frac{d!}{\sqrt{d+1}}\right)}^{\frac{1}{d}},\hat{\rho }=\sqrt{\frac{2}{d(d+1)}}\hat{h},\hat{a}=\sqrt{\frac{d+1}{2d}}\hat{h}.$$

$$\frac{h_{K,{\mathbb{D}}^{-1}}^2}{{\hat{h}}^2}\le {\Vert {(F_K^{\prime })}^T{\mathbb{D}}_K^{-1}{F}_K^{\prime }\Vert }_2\le \frac{h_{K,{\mathbb{D}}^{-1}}^2}{{\hat{\rho }}^2},$$

$$\frac{{\hat{a}}^2}{a_{K,{\mathbb{D}}^{-1}}^2}\le {\Vert {(F_K^{\prime })}^{-1}{\mathbb{D}}_K{(F_K^{\prime })}^{-T}\Vert }_2\le \frac{d^2{\hat{a}}^2}{a_{K,{\mathbb{D}}^{-1}}^2}.$$

$${\hat{\rho }}^2{C}_{\hat{\nabla }}{h}_{{\mathbb{D}}^{-1}}^{-2}\sum _{K\in {\omega }_j}\left|K\right|{\left(\frac{h_{{\mathbb{D}}^{-1}}}{h_{K,{\mathbb{D}}^{-1}}}\right)}^2\le A_{jj}\le d^2{\hat{a}}^2{C}_{\hat{\nabla }}{h}_{{\mathbb{D}}^{-1}}^{-2}\sum _{K\in {\omega }_j}\left|K\right|{\left(\frac{h_{{\mathbb{D}}^{-1}}}{a_{K,{\mathbb{D}}^{-1}}}\right)}^2.$$

$$A_{jj}$$

$${\lambda }_{\min }({(F_K^{\prime })}^{-1}{\mathbb{D}}_K{(F_K^{\prime })}^{-T})={\Vert {(F_K^{\prime })}^T{\mathbb{D}}_K^{-1}{F}_K^{\prime }\Vert }_2^{-1}\ge {\hat{\rho }}^2{h}_{K,{\mathbb{D}}^{-1}}^{-2}.$$

$${\lambda }_{\min }(A)\gtrsim d_{\min }{N}^{-1}{\Psi }_E^{-1},$$

$${\Psi }_E=\{\begin{array}{ll}1,& d=1,\\ 1+\ln \left(\frac{|\bar{K}|}{|K_{\min }|}\right),& d=2,\\ {\left(\frac{1}{N}\sum _{K\in {\mathcal{T}}_h}{\left(\frac{|\bar{K}|}{|K|}\right)}^{\frac{d-2}{2}}\right)}^{\frac{2}{d}},& d\ge 3,\end{array}$$

$${\lambda }_{\max }(P^{-\frac{1}{2}}M{P}^{-\frac{1}{2}})\le \{\begin{array}{ll}\frac{|{\omega }_{\max }|}{d+1},& P=I_{N_{vi}},\\ 1+\frac{d}{2},& P=M_D,\text{leavevmode}{M}_{lump},\text{leavevmode}{M}_D+\Delta t{A}_D,\end{array}$$

$${\lambda }_{\min }(P^{-\frac{1}{2}}M{P}^{-\frac{1}{2}})\ge \{\begin{array}{ll}\frac{|{\omega }_{\min }|}{(d+1)(d+2)},& P=I_{N_{vi}},\\ 1/2,& P=M_D,\\ 1/(d+2),& P=M_{lump},\\ C{\left(1+\Delta t{h}_{{\mathbb{D}}^{-1}}^{-2}\underset{j}{\max }\sum _{K\in {\omega }_j}\frac{|K|}{|{\omega }_j|}{\left(\frac{h_{{\mathbb{D}}^{-1}}}{a_{K,{\mathbb{D}}^{-1}}}\right)}^2\right)}^{-1}& P=M_D+\Delta t{A}_D,\end{array}$$

$$P_{jj}\lesssim \left|{\omega }_j\right|+\Delta t{h}_{{\mathbb{D}}^{-1}}^{-2}\sum _{K\in {\omega }_j}\left|K\right|{\left(\frac{h_{{\mathbb{D}}^{-1}}}{a_{K,{\mathbb{D}}^{-1}}}\right)}^2.$$

$${\mathbf{v}}^T{P}^{-\frac{1}{2}}M{P}^{-\frac{1}{2}}\mathbf{v}\gtrsim \sum _j\frac{v_j^2\left|{\omega }_j\right|}{P_{jj}}\gtrsim \sum _j{v}_j^2{\left(1+\Delta t{h}_{{\mathbb{D}}^{-1}}^{-2}\sum _{K\in {\omega }_j}\frac{\left|K\right|}{\left|{\omega }_j\right|}{\left(\frac{h_{{\mathbb{D}}^{-1}}}{a_{K,{\mathbb{D}}^{-1}}}\right)}^2\right)}^{-1}.\text{qed}$$

$${\lambda }_{\max }(P^{-\frac{1}{2}}A{P}^{-\frac{1}{2}})\le \{\begin{array}{ll}C{h}_{{\mathbb{D}}^{-1}}^{-2}\underset{j}{\max }\sum _{K\in {\omega }_j}\left|K\right|{\left(\frac{h_{{\mathbb{D}}^{-1}}}{a_{K,{\mathbb{D}}^{-1}}}\right)}^2,& P=I_{N_{vi}},\\ C{h}_{{\mathbb{D}}^{-1}}^{-2}\underset{j}{\max }\sum _{K\in {\omega }_j}\frac{|K|}{|{\omega }_j|}{\left(\frac{h_{{\mathbb{D}}^{-1}}}{a_{K,{\mathbb{D}}^{-1}}}\right)}^2,& P=M,\text{leavevmode}{M}_D,\text{leavevmode}{M}_{lump},\\ (d+1)\Delta t{}^{-1},& P=M_D+\Delta t{A}_D,\end{array}$$

$${\lambda }_{\min }(P^{-\frac{1}{2}}A{P}^{-\frac{1}{2}})\gtrsim \{\begin{array}{ll}{d}_{\min }{N}^{-1}{\Psi }_E^{-1},& P=I_{N_{vi}},\\ {\lambda }_{\mathbb{D}},& P=M,\text{leavevmode}{M}_D,\text{leavevmode}{M}_{lump},\\ d_{\min }{\Psi }_{\mathbb{D}}^{-1},& P=M_D+\Delta t{A}_D,\end{array}$$

Full text

For the simultaneous solution in case of a positive semidefinite symmetric part of the RK coefficient matrix and for the successive solution it is shown that the conditioning of an implicit RK method behaves like that of the implicit Euler method. If the smallest eigenvalue of the symmetric part of the RK coefficient matrix is negative and the simultaneous solution strategy is used, an upper bound on the time step is given so that the system matrix is positive definite.

The obtained bounds for the condition number have explicit geometric interpretations and take the interplay between the diffusion matrix and the mesh geometry into full consideration. They show that there are three mesh-dependent factors that can affect the conditioning: the number of elements, the mesh nonuniformity measured in the Euclidean metric, and the mesh nonuniformity with respect to the inverse of the diffusion matrix. They also reveal that the preconditioning using the diagonal of the system matrix, the mass matrix, or the lumped mass matrix can effectively eliminate the effects of the mesh nonuniformity measured in the Euclidean metric. Illustrative numerical examples are given.

 

This is a preprint of a contibution published by Elsevier in J. Comput. Appl. Math. (2019) 112497.

© 2019. Licensed under CC-BY-NC-ND 4.0 (https://creativecommons.org/licenses/by-nc-nd/4.0).

The formal publication is available online at https://doi.org/10.1016/j.cam.2019.112497.

††journal: Journal of Computational and Applied Mathematics

1 Introduction

The nonuniformity of adaptive meshes has considerable effects on the conditioning of the discrete approximation of partial differential equations (PDEs) and their efficient solution. To study these effects, we investigate the implicit Runge–Kutta (RK) integration for the linear finite element (FE) approximation of linear diffusion equations on general simplicial anisotropic meshes for the initial-boundary value problem (IBVP)

[TABLE]

where $\Omega\subset\mathbb{R}^{d}$ ($d\geq 1$) is a bounded polygonal or polyhedral domain, $u_{0}$ is a given initial solution, and $\mathbb{D}$ is the diffusion matrix. We assume that $\mathbb{D}=\mathbb{D}(\boldsymbol{x})$ is time independent and symmetric and uniformly positive definite on $\Omega$, i.e.,

[TABLE]

where the less-than-or-equal sign means that the difference between the two matrices is positive semidefinite. We consider the Dirichlet boundary condition in this work but the analysis is applicable to other boundary condition types without major modifications.

Much effort has been made in the past to understand the effects of mesh nonuniformity on the conditioning of FE approximations. For example, Fried [8] obtains a bound on the condition number of the stiffness matrix for the linear FE approximation of the Laplace operator for general meshes. For the Laplace operator on isotropic adaptive grids, Bank and Scott [2] show that the condition number of the diagonally scaled stiffness matrix is essentially the same as for a regular mesh. Ainsworth, McLean, and Tran [1] and Graham and McLean [9] extend this result to the boundary element equations for locally quasi-uniform meshes and provide a bound in terms of patch volumes and aspect ratios. Du et al. [5] obtain a bound on the condition number of the stiffness matrix for a general diffusion operator on anisotropic meshes which reveals the relation between the condition number and some mesh quality measures. For the FE approximations of parabolic problems, Zhu and Du [21, 22] develop mesh dependent stability and condition number estimates for the explicit and implicit Euler methods. In the case of a lumped mass matrix ($M_{lump}$), Shewchuk [20] provides a bound on the largest eigenvalue of $M_{lump}^{-1}A$ in terms of the maximum eigenvalues of local element matrices, where $A$ is the stiffness matrix. The results mentioned above allow anisotropic adaptive meshes but do not fully take into account the interplay between the mesh geometry and the diffusion matrix. (An exception is the bound by Shewchuk which includes the effects of the diffusion coefficients.) Moreover, the existing analysis employs mesh restrictions in form of mesh regularity assumptions, e.g., local mesh uniformity, or parameters in final estimates that are related to mesh regularity, such as the maximum ratio of volumes of neighbouring elements or the maximum number of elements in a patch.

The objective of this work is to develop estimates on the conditioning of the resulting linear system that take the interplay between the mesh geometry and the diffusion matrix into full consideration, have explicit geometric interpretation, and make no prior assumptions on the mesh regularity. This is a continuation of our previous effort to develop bounds for the condition number of the stiffness matrix for the linear FE equations of a general diffusion operator on arbitrary anisotropic meshes [17, 18] and the largest permissible time steps for explicit RK schemes for both linear and high order FE approximations of the IBVP Eq. 1 [15, 16]. In particular, these bounds show [18] that the condition number of the stiffness matrix depends on three factors: the factor depending on the number of mesh elements and corresponding to the condition number of the linear FE equations for the Laplace operator on a uniform mesh, the nonuniformity of the mesh viewed in the metric defined by the inverse diffusion matrix, $\mathbb{D}^{-1}$, and the mesh nonuniformity measured in the Euclidean metric. Moreover, the Jacobi preconditioning, an optimal diagonal scaling for a symmetric positive definite sparse matrix [12, Corollary 7.6ff.], can effectively eliminate the effects of mesh nonuniformity and reduce those of the mesh nonuniformity with respect to $\mathbb{D}^{-1}$ [18]. Detailed characterizations of the condition number according to the mesh concentration distribution can be obtained using Green’s functions [8, 17].

We first consider the implicit Euler method (the simplest implicit RK method) and establish bounds for the condition number of the corresponding system matrix with and without diagonal scaling.

For general implicit RK methods, we consider two strategies for solving the resulting system: the simultaneous solution and a successive solution. We need to point out that for general implicit RK methods, the coefficient matrix is not necessarily symmetric and thus a condition number in the standard definition will not provide much information for iterative solvers. Motivated by the convergence estimates by Eisenstat et al. [6] for the generalized minimal residual method (GMRES) for nonsymmetric systems, we use here the ratio of the maximum singular value of the coefficient matrix and the minimum eigenvalue of its symmetric part to measure the conditioning.

For the simultaneous solution, the system is solved as a whole. Obtained estimates reveal that the conditioning of implicit RK methods with positive semidefinite symmetric part of the coefficient matrix behaves like that of the implicit Euler method (Theorem 5.1 in Section 5). If the smallest eigenvalue of the symmetric part of the RK coefficient matrix is negative, we provide an upper bound Eq. 45 on the possible time step so that the system matrix is positive definite. This condition can be serious and lead to $\Delta t=\mathcal{O}(h^{2})$.

For the successive solution, which follows the commonly used implementation of implicit RK methods [3, 4], the system is first transformed into a number of smaller systems using the Jordan normal form of the RK matrix and then solved successively. We show that the conditioning of the implicit RK integration behaves like that of the implicit Euler method and is determined by the conditioning of two types of matrices. The first one is similar to the implicit Euler method and corresponds to the real eigenvalues of the RK matrix, while the second, twice as large, corresponds to the complex eigenvalues of the RK matrix (cf. Eq. 52 in Section 5).

The paper is organized as follows. We first introduce the FE formulation and its implicit RK integration (Section 2) and provide preliminary estimates for the extremal eigenvalues of the mass and stiffness matrices (Section 3). The main results for the conditioning of the coefficient matrices are given in Sections 4 and 5, followed by numerical examples (Section 6) and conclusions (Section 7).

2 Linear FE approximation and implicit RK integration

Let $\{\mathcal{T}_{h}\}$ be a given family of simplicial meshes for the domain $\Omega$ and $N$, $N_{v}$, and $N_{vi}$ the number of mesh elements, vertices, and interior vertices, respectively. For convenience, we assume that the vertices are ordered such that the first $N_{vi}$ vertices are the interior ones. The element patch associated with the $j$-th vertex is denoted by $\omega_{j}$, $K$ denotes a given mesh element and $\hat{K}$ is the reference element which is assumed to have been taken as a unitary equilateral simplex. Element and patch volumes are denoted by ${\left\lvert K\right\rvert}$ and ${\left\lvert\omega_{j}\right\rvert}=\sum_{K\in\omega_{j}}{\left\lvert K\right\rvert}$. For each mesh element $K\in\mathcal{T}_{h}$, we denote the invertible affine mapping from $\hat{K}$ to $K$ and its Jacobian matrix by $F_{K}$ and $F_{K}^{\prime}$, respectively. Note that $F_{K}^{\prime}$ is a constant matrix and $\det(F_{K}^{\prime})={\left\lvert K\right\rvert}$.

Let $V^{h}\subset H_{0}^{1}(\Omega)$ be the linear FE space associated with $\mathcal{T}_{h}$. The piecewise linear FE solution $u_{h}(t)\in V^{h}$, $t>0$ for Eq. 1 is defined by

[TABLE]

subject to the initial condition

[TABLE]

It can be expressed as

[TABLE]

where $\phi_{j}$ is the linear basis function associated with the $j$-th vertex. Inserting this into Eq. 3 and taking $v_{h}=\phi_{k}$, $k=1,\dotsc,N_{vi}$, successively yields the system

[TABLE]

where $\boldsymbol{u}={(u_{1},\dotsc,u_{N_{vi}})}^{T}$ and $M$ and $A$ are the mass and stiffness matrices with

[TABLE]

For the time integration of Eq. 4 we consider a general implicit $s$-stage RK method with the Butcher tableau

[TABLE]

and assume that the eigenvalues of the RK matrix have nonnegative real parts. This requirement is satisfied by most implicit RK methods (e.g., see Table 1). Applying the method to Eq. 4 yields

[TABLE]

where $\boldsymbol{u}^{n}$ and $\boldsymbol{u}^{n+1}$ are the approximations of $\boldsymbol{u}(t_{n})$ and $\boldsymbol{u}(t_{n+1})$. The major cost of finding $\boldsymbol{u}^{n+1}$ in the above method is the solution of Eq. 6 for $\boldsymbol{v}_{1},\dotsc,\boldsymbol{v}_{s}$. Using the Kronecker matrix product $\otimes$ (e.g., see [19]), the $s\times s$ identity matrix $I_{s}$, and $\Gamma:={(\gamma_{kj})}_{k,j=1}^{s}$, the coefficient matrix of Eq. 6 can be expressed as

[TABLE]

The simplest implicit RK method is the implicit Euler method with $s=1$, $\Gamma=1$, $c_{1}=1$, and $b_{1}=1$ for which Eqs. 6 and 7 are reduced to

[TABLE]

In Section 4 we study the conditioning of the coefficient matrix $M+\Delta{}t{}\,A$ related to the efficient iterative solution of Eq. 9. The system Eq. 9 can also be solved using a symmetric and positive definite preconditioner $P$, which leads to

[TABLE]

The efficient iterative solution of Eq. 10 is related to the conditioning of $P^{-\frac{1}{2}}(M+\Delta{}t{}\,A)P^{-\frac{1}{2}}$. In this work, we consider as preconditioners the mass matrix $M$, its diagonal part $M_{D}$, the lumped mass matrix $M_{lump}$, and the diagonal part $M_{D}+\Delta{}t{}\,A_{D}$ of $M+\Delta{}t{}\,A$. The choice $P=M$ is of theoretical importance but impractical since $M^{-\frac{1}{2}}$ is expensive to compute. The other choices are simple and economic to implement.

The results for the Euler method will be used in Section 5 to study the conditioning of Eq. 8 and its diagonally preconditioned version for general implicit RK methods.

3 Preliminary estimates on the extreme eigenvalues

of the mass and stiffness matrices

Hereafter, $C$ denotes a generic constant which may have different values at different appearances and may depend on the dimension, the choice of the reference element, and the reference basis linear functions but is independent of the mesh and the IBVP coefficients. For notation simplicity, when using this generic constant $C$, we will sometime write $a\gtrsim b$ and $a\lesssim b$ meaning $a\geq Cb$ and $a\leq Cb$, respectively.

Lemma 3.1** ([18, proof of Theorem 3.1]).**

The mass matrix $M$ and its diagonal part $M_{D}$ are related by

[TABLE]

Lemma 3.2** ([16, Lemma 2.3]).**

The mass matrix $M$ and the lumped mass matrix $M_{lump}$ are related by

[TABLE]

Lemma 3.3** ([18, Lemma 4.1]).**

The stiffness matrix $A$ and its diagonal part $A_{D}$ satisfy

[TABLE]

Lemma 3.4** ([16, Lemma 2.5]).**

The diagonal entries of the stiffness matrix are bounded by

[TABLE]

where $\mathbb{D}_{K}=\frac{1}{{\left\lvert K\right\rvert}}\int_{K}\mathbb{D}(\boldsymbol{x})\;\mathrm{d}\boldsymbol{x}$ is the average of $\mathbb{D}$ on $K$ and $C_{\hat{\nabla}}=\frac{d+1}{d}{\left(\frac{d!}{\sqrt{d+1}}\right)}^{\frac{2}{d}}$.

The next two lemmas establish bounds for $A_{jj}$ with a more explicit geometric interpretation than Eq. 14. For this, we denote the diameter and the minimal height of $K$ in the metric $\mathbb{D}_{K}^{-1}$ by $h_{K,\mathbb{D}^{-1}}$ and $a_{K,\mathbb{D}^{-1}}$, respectively. The average element diameter is defined as

[TABLE]

Further, let $\hat{h}$, $\hat{\rho}$, and $\hat{a}$ be the diameter, the in-diameter, and the minimal height of the unitary equilateral $\hat{K}$, respectively, i.e.,

[TABLE]

Lemma 3.5** ([14, Lemmas 4.1 and 4.2]).**

It holds

[TABLE]

Lemma 3.6**.**

It holds

[TABLE]

Proof.

Combining Eqs. 17 and 14 yields

[TABLE]

which gives the right inequality of Eq. 18. On the other hand, Eq. 16 gives

[TABLE]

Combining this with Eq. 14 leads to the left inequality of Eq. 18. ∎

A mesh that is uniform in the metric $\mathbb{D}^{-1}$ (a $\mathbb{D}^{-1}$-uniform mesh) satisfies $h_{K,\mathbb{D}^{-1}}\sim a_{K,\mathbb{D}^{-1}}\sim h_{\mathbb{D}^{-1}}$. In such a case, Eq. 18 implies $A_{jj}\sim{\left\lvert\omega_{j}\right\rvert}h_{\mathbb{D}^{-1}}^{-2}\sim{\left\lvert\omega_{j}\right\rvert}N^{\frac{2}{d}}$.

Lemma 3.7** ([18, Lemma 5.1]).**

The smallest eigenvalue of the stiffness matrix is bounded by

[TABLE]

where

[TABLE]

${\left\lvert K_{\min}\right\rvert}$* and ${\lvert\bar{K}\rvert}=\frac{1}{N}{\lvert\Omega\rvert}$ are the minimal and the average element volumes and $d_{\min}$ is the global smallest eigenvalue of $\mathbb{D}$ (cf. Eq. 2).*

The factor $\Psi_{E}$ shows that the dependence of the lower bound on the mesh non-uniformity is mild although getting stronger in higher dimensions: $\Psi_{E}$ is mesh-independent in 1d, contains only $\ln\left({{\lvert\bar{K}\rvert}}/{{\left\lvert K_{\min}\right\rvert}}\right)$ in 2d, and involves the generalized $d/2$-mean of the ratio ${\left({{\lvert\bar{K}\rvert}}/{{\left\lvert K\right\rvert}}\right)}^{(d-2)/{d}}$ over all elements for $d\geq 3$.

The next lemmas summarize bounds for the extremal eigenvalues of $M$ and $A$.

Lemma 3.8**.**

The extremal eigenvalues of $P^{-\frac{1}{2}}MP^{-\frac{1}{2}}$ are bounded by

[TABLE]

and

[TABLE]

where ${\left\lvert\omega_{\max}\right\rvert}$ and ${\left\lvert\omega_{\min}\right\rvert}$ are the maximal and minimal patch volumes, respectively.

Proof.

The inequalities Eq. 21 and Eq. 22 for $P=I_{N_{vi}}$ (no preconditioning), $P=M_{D}$, and $M_{lump}$ follow from Lemmas 3.1 and 3.2. For $P=M_{D}+\Delta{}t{}\,A_{D}$, Eqs. 11 and 18 give

[TABLE]

Then, Eq. 22 follows from

[TABLE]

Lemma 3.9**.**

The extremal eigenvalues of $P^{-\frac{1}{2}}AP^{-\frac{1}{2}}$ are bounded by

[TABLE]

and

[TABLE]

where $\lambda_{\mathbb{D}}$ is the minimal eigenvalue of the operator $-\nabla\cdot(\mathbb{D}\nabla)$ and

[TABLE]

Proof.

For $P=I_{N_{vi}}$, Eq. 24 follows from Lemmas 3.3 and 3.6. For $P=M$, Eqs. 11 and 13 give

[TABLE]

which, together with Eqs. 11 and 18, implies Eq. 24. For $P=M_{D}$ and $P=M_{lump}$, Eq. 24 follows from Lemmas 3.1 and 3.2. For $P=M_{D}+\Delta{}t{}\,A_{D}$, Eq. 24 follows from Eq. 13.

The inequality Eq. 25 for $P=I_{N_{vi}}$ follows from Lemma 3.7. For $P=M$ and, similarly for $P=M_{D}$ and $P=M_{lump}$, Eq. 25 follows from the basic property of the conformal FE approximation of elliptic eigenvalue problems (e.g., [7, Theorem 1]), since the eigenvalue problem for $P^{-\frac{1}{2}}AP^{-\frac{1}{2}}$ is a FE approximation to the eigenvalue problem for the operator $-\nabla\cdot(\mathbb{D}\nabla)$.

For $P=M_{D}+\Delta{}t{}\,A_{D}$, using the strategy in the proof of [18, Lemma 5.1], we have

[TABLE]

with a parameter $q>2$ and some not-all-zero nonnegative numbers $s_{K}$, $K\in\mathcal{T}_{h}$ (to be determined later). Here, $\sum\limits_{K}$ denotes the summation over all elements in $\mathcal{T}_{h}$.

For notational simplicity, we denote $r_{K}=a_{K,\mathbb{D}^{-1}}^{-2}$ and rewrite Eq. 23 as

[TABLE]

For $d=1$, Eq. 29 implies

[TABLE]

and therefore

[TABLE]

For $d=2$, we choose $s_{K}={\left\lvert K\right\rvert}^{1-\frac{2}{q}}(1+\Delta{}t{}\,r_{K})$. Then, Eq. 29 implies

[TABLE]

and using Eq. 28 we obtain

[TABLE]

The denominator tends to infinity as $q\to\infty$ and to $2\max_{K}(1+\Delta{}t{}\,r_{K})$ as $q\to 2^{+}$. To find an optimal choice for $q$, we first estimate the denominator as

[TABLE]

where

[TABLE]

If $\psi_{\mathbb{D}}>e$, we choose $q=2\ln\psi_{\mathbb{D}}$ and obtain $q\psi_{\mathbb{D}}^{\frac{2}{q}}=2e\ln\psi_{\mathbb{D}}$. If $\psi_{\mathbb{D}}\leq e$, we use $q\to 2^{+}$ and obtain $q\psi_{\mathbb{D}}^{\frac{2}{q}}\leq 2e$. Combining these two cases yields

[TABLE]

and therefore

[TABLE]

For $d\geq 3$, we choose $s_{K}={\left\lvert K\right\rvert}^{\frac{2}{d}}(1+\Delta{}t{}\,r_{K})$ such that $P_{jj}^{-1}\sum\limits_{K\in\omega_{j}}r_{K}{\left\lvert K\right\rvert}^{\frac{2}{q}}\geq C$. From Eq. 28, we have

[TABLE]

which gives

[TABLE]

Further, Eq. 31 can be rewritten into Eq. 27 and

[TABLE]

Combining these with Eqs. 30 to 33 gives Eq. 25 with $P=M_{D}+\Delta{}t{}\,A_{D}$. ∎

4 Conditioning of the implicit Euler integration

In the following we estimate the conditioning of $M+\Delta{}t{}\,A$ (in the $l_{2}$-norm), which is related to the efficient iterative solution of Eq. 9.

4.1 Condition number of $M+\Delta{}t{}\,A$

Theorem 4.1**.**

The condition number of $M+\Delta{}t{}\,A$ is bounded by

[TABLE]

Proof.

Equations 11 and 13 yield

[TABLE]

From this and Eqs. 11 and 18 we obtain

[TABLE]

On the other hand, using $\lambda_{\min}(M+\Delta{}t{}\,A)\geq\lambda_{\min}(M)+\Delta{}t{}\,\lambda_{\min}(A)$ and Lemmas 3.8 and 3.9 we have

[TABLE]

Combining the above results yields Eq. 34. ∎

There are three factors influencing bound Eq. 34. The first factor is the number of the mesh elements $N$. The second factor is

[TABLE]

which reflects the mesh nonuniformity in the metric $\mathbb{D}^{-1}$ and is a constant for a $\mathbb{D}^{-1}$-uniform mesh, for which $a_{K,\mathbb{D}^{-1}}\sim h_{\mathbb{D}^{-1}}$ for all $K$. The third factor is the effect of the mesh nonuniformity (in the Euclidean metric) reflected by $\Psi_{E}$, ${{\left\lvert\omega_{j}\right\rvert}}/{{\left\lvert\omega_{\min}\right\rvert}}$, and $N{\left\lvert\omega_{\min}\right\rvert}$, which all become constants if the mesh is uniform.

The time step size $\Delta{}t{}$ plays the role of a homotopy parameter between the mass and the stiffness matrices:

[TABLE]

Thus, $\kappa(M)$ depends on the nonuniformity of the element patch volumes. On the other hand,

[TABLE]

which has been obtained previously [18, Theorem 5.2]. It depends on three factors as well: $N$ (through $h_{\mathbb{D}^{-1}}$), the mesh nonuniformity in the Euclidean metric (through $\Psi_{E}$), and the mesh nonuniformity with respect to $\mathbb{D}^{-1}$ (through the term in the maximum).

4.2 Condition number of $P^{-\frac{1}{2}}(M+\Delta{}t{}\,A)P^{-\frac{1}{2}}$

Theorem 4.2**.**

The condition number of $P^{-\frac{1}{2}}(M+\Delta{}t{}\,A)P^{-\frac{1}{2}}$ with $P=M$, $M_{D}$, or $M_{lump}$ is bounded by

[TABLE]

where $\lambda_{\mathbb{D}}$ is the minimal eigenvalue of $-\nabla\cdot(\mathbb{D}\nabla)$.

Proof.

Bound Eq. 38 is obtained by combining Lemmas 3.8 and 3.9 with the estimates

[TABLE]

The bounds on $\kappa(P^{-\frac{1}{2}}(M+\Delta{}t{}\,A)P^{-\frac{1}{2}})$ for $P=M$, $M_{D}$, and $M_{lump}$ are similar, while the last two choices lead to a simple diagonal scaling, which is easier to implement.

All three choices reduce the effects of the mesh nonuniformity: in comparison to Eq. 34, bound Eq. 38 does not depend on ${\left\lvert\omega_{\min}\right\rvert}$ or $\Psi_{E}$ directly and contains only the $\mathbb{D}^{-1}$-nonuniformity factor Eq. 35. This is intuitive, since the eigenvalues of $M^{-1}A$ approximate those of the underlying continuous operator, which are mesh-independent. However, $\kappa(M^{-1}A)$ is not necessarily smaller then $\kappa(A)$ and the overall effect depends on the magnitude of $\Delta{}t{}$: if $\Delta{}t{}$ is large, $\kappa(M^{-1}(M+\Delta{}t{}A))$ might not be better than $\kappa(M+\Delta{}t{}A)$. On the other hand, if $\Delta{}t{}$ is small, we can expect that $\kappa(M^{-1}(M+\Delta{}t{}A))<\kappa(M+\Delta{}t{}A)$. The numerical experiments in Section 6 support this argument (see Examples 6.2 and 3(b)).

Theorem 4.3**.**

The condition number of $P^{-\frac{1}{2}}(M+\Delta{}t{}\,A)P^{-\frac{1}{2}}$ with the Jacobi preconditioner $P=M_{D}+\Delta{}t{}\,A_{D}$ is bounded by

[TABLE]

where $\Psi_{\mathbb{D}}$ is given in Eq. 26.

Proof.

The proof is similar to that of Theorem 4.2. ∎

Bound Eq. 39 is comparable to Eq. 38 although the former is smaller than the latter in general, especially for large $\Delta{}t{}$, since the factor $\Psi_{\mathbb{D}}$ in Eq. 39 involves averaging over all elements, whereas Eq. 38 involves the maximum over patch averages.

In Eq. 39, $\Delta{}t{}$ is a homotopy parameter between the mass and stiffness matrices:

[TABLE]

On the other hand,

[TABLE]

where $\psi_{\mathbb{D}}$ becomes

[TABLE]

This bound is equivalent to the bound obtained in [18, Theorem 5.2].

Remark 4.1**.**

In comparison to the bound Eq. 34 for $M+\Delta{}t{}\,A$, bounds Eqs. 38 and 39 contain only two mesh-dependent factors: $N$ and the mesh nonuniformity in the metric $\mathbb{D}^{-1}$ (through the terms involving the ratio $h_{\mathbb{D}^{-1}}/a_{K,\mathbb{D}^{-1}}$). This shows that the effects of the mesh nonuniformity (in the Euclidean metric) on the condition number is effectively eliminated by the preconditioning.

Remark 4.2**.**

For a $\mathbb{D}^{-1}$-uniform mesh, $a_{K,\mathbb{D}^{-1}}\sim h_{\mathbb{D}^{-1}}$ and, hence, all terms involving the ratio $h_{\mathbb{D}^{-1}}/a_{K,\mathbb{D}^{-1}}$ will become a constant. Thus,

[TABLE]

which shows more clearly the role of $\Delta{}t{}$:

[TABLE]

5 Conditioning of general implicit RK integration

For a general implicit RK method, the matrix $\mathbb{A}$ in Eq. 8 is not necessarily normal and in this case a condition number in the standard definition does not provide much information for the convergence of iterative methods. When its symmetric part, $(\mathbb{A}+\mathbb{A}^{T})/2$, is positive definite, it is known [6] that the convergence of the generalized minimal residual method (GMRES) is

[TABLE]

where $\sigma_{\max}(\mathbb{A})$ is the largest singular value of $\mathbb{A}$, $\lambda_{\min}((\mathbb{A}+\mathbb{A}^{T})/2)$ is the minimal eigenvalue of the symmetric part, $r_{n}$ is the residual of the corresponding linear system at the $n$-th iterate, and $\|\cdot\|$ stands for the matrix or vector 2-norm. Thus, we can consider the “condition number”

[TABLE]

This definition reduces to the standard definition of the condition number (in 2-norm) for symmetric matrices.

We first recall that for an arbitrary square matrix $U$, its maximal singular value is given by

[TABLE]

From this, for any square matrices $U$ and $V$ we have

[TABLE]

The next lemma provides the eigenvalues and singular values of the Kronecker product of any square matrices $U$ and $V$.

Lemma 5.1** ([19, Theorems 13.10 and 13.12]).**

For any square matrices $U$ and $V$, the eigenvalues and singular values of $U\otimes V$ are $\lambda_{j}(U)\lambda_{k}(V)$ and $\sigma_{j}(U)\sigma_{k}(V)$, $j,k=1,2,\dotso$

The above lemma in particular implies

[TABLE]

and, for a positive semi-definite matrix $V$,

[TABLE]

5.1 General implicit RK methods: simultaneous solution

Although a successive strategy for implicit RK methods is more common (e.g., [10, p. 131]), the simultaneous solution can be appealing if an iterative method is used. It is also important to have a theoretical understanding of the overall system conditioning.

Theorem 5.1**.**

Let $P$ be a symmetric and positive definite preconditioner, $\Gamma$ an $s\times s$ implicit RK coefficient matrix, and $\mathbb{A}=I_{s}\otimes M+\Delta{}t\,\Gamma A$.

If $\lambda_{\min}((\Gamma+\Gamma^{T})/2)<0$ and $\Delta{}t$ is chosen sufficiently small such that

[TABLE]

then

[TABLE]

If $\lambda_{\min}((\Gamma+\Gamma^{T})/2)\geq 0$, then

[TABLE]

Proof.

First, notice that

[TABLE]

Then, from Eqs. 43 and 5.1 we have

[TABLE]

Moreover, the symmetric part of $(I_{s}\otimes P^{-\frac{1}{2}})(I_{s}\otimes M+\Delta{}t{}\,\Gamma\otimes A)(I_{s}\otimes P^{-\frac{1}{2}})$ is

[TABLE]

Then,

[TABLE]

In case of $\lambda_{\min}((\Gamma+\Gamma^{T})/2)\geq 0$ or if $\lambda_{\min}((\Gamma+\Gamma^{T})/2)<0$ and $\Delta{}t$ is chosen such that Eq. 45 holds, then the symmetric part of $(I_{s}\otimes P^{-\frac{1}{2}})(I_{s}\otimes M+\Delta{}t{}\,\Gamma\otimes A)(I_{s}\otimes P^{-\frac{1}{2}})$ is positive definite and Eqs. 47 and 46 follow from the above results and Eqs. 42 and 44. ∎

Remark 5.1**.**

In case of $P=I_{N_{vi}}$, $M$, $M_{D}$, $M_{lump}$, and $M_{D}+\Delta{}t{}\,A_{D}$, the estimates for the eigenvalues $\lambda_{\min}(P^{-\frac{1}{2}}MP^{-\frac{1}{2}})$, $\lambda_{\max}(P^{-\frac{1}{2}}MP^{-\frac{1}{2}})$, $\lambda_{\min}(P^{-\frac{1}{2}}AP^{-\frac{1}{2}})$, and $\lambda_{\max}(P^{-\frac{1}{2}}AP^{-\frac{1}{2}})$ are given by Lemmas 3.8 and 3.9. The quantities $\sigma_{\max}(\Gamma)$ and $\lambda_{\min}((\Gamma+\Gamma^{T})/2)$ are given in Table 1 for five implicit Runge–Kutta methods. One can see that $\sigma_{\max}(\Gamma)=\mathcal{O}(1)$ but $\lambda_{\min}((\Gamma+\Gamma^{T})/2)$ can be zero, positive, and negative.

For the 4th-order Gauss and 3rd-order Radau IA methods, $\lambda_{\min}((\Gamma+\Gamma^{T})/2)\geq 0$ and Eq. 45 poses no constraint on the choice of $\Delta{}t$. Equation 47 indicates that the conditioning of the system Eq. 8 resulting from implicit RK integration behaves more or less like that of $M+\Delta{}t{}\,A$ for the implicit Euler method, and the bound also involves three mesh-dependent factors and one of them can be eliminated effectively by diagonal preconditioning (see Remark 4.1).

For the other three methods in Table 1 with $\lambda_{\min}((\Gamma+\Gamma^{T})/2)<0$, Eq. 45 gives an upper bound on possible $\Delta{}t$. Note that this condition can be serious and lead to a condition like $\Delta t=\mathcal{O}(h^{2})$.

5.2 Diagonally implicit RK (DIRK) methods

The coefficient matrix $\Gamma={(\gamma_{kj})}_{k,j=1}^{s}$ of a DIRK method is a lower triangular matrix and the system Eq. 6 is solved by successively solving $s$ linear systems with $M+\Delta{}t{}\,\gamma_{jj}A$, $j=1,\dotsc,s$. The conditioning of $M+\Delta{}t{}\,\gamma_{jj}A$ is therefore similar to that of $M+\Delta{}t{}\,A$ and, hence, the analysis in Section 4 also applies to DIRK methods.

5.3 General implicit RK methods: successive solution

A successive solution procedure transforms the large system Eq. 6 into a number of smaller systems, which are then solved successively. We adopt the approach of Butcher [4] and Bickart [3] and carry out the transformation using the Jordan normal form of the RK matrix. To keep our analysis applicable to methods with a singular $\Gamma$, we use the Jordan normal form of $\Gamma$ instead of $\Gamma^{-1}$, which is the conventional choice (e.g., see [10, p. 131]).

Let the Jordan normal form of $\Gamma$ be

[TABLE]

where $T$ is a real invertible matrix and $J_{j}$, $j=1,\dotsc,p$, are the Jordan blocks, which either have the form

[TABLE]

or

[TABLE]

Recall that the eigenvalues of $\Gamma$ are assumed to have nonnegative real parts, i.e., $\mu_{j},\alpha_{j}\geq 0$.

Theorem 5.2**.**

Assume that the eigenvalues of the coefficient matrix $\Gamma$ for a given implicit RK method have nonnegative real parts and the system Eq. 6 is solved by using the Jordan normal form of $\Gamma$ to transform it into smaller systems. Then, the conditioning of implicit RK integration of Eq. 4 is determined by the conditioning of

[TABLE]

where $\mu_{j}$ and $\alpha_{j}\pm i\beta_{j}$ are the real and complex eigenvalues of $\Gamma$, respectively.

Moreover, for any symmetric and positive definite preconditioner $P$,

[TABLE]

Proof.

Using Eq. 49, we can rewrite Eq. 8 as

[TABLE]

where $n_{j}$ is the size of $J_{j}$ (with $\sum_{j=1}^{p}n_{j}=s$). Since the systems with the coefficient matrix $T\otimes I_{N_{vi}}$ or $T^{-1}\otimes I_{N_{vi}}$ can be solved directly and efficiently, we only need to consider the iterative solution of the systems associated with

[TABLE]

For a Jordan block in the form Eq. 50, the above matrix becomes

[TABLE]

which can be solved by successively solving $n_{j}$ systems with $M+\mu_{j}\Delta{}t{}\,A$ (backward substitution).

On the other hand, for a Jordan block in the form Eq. 51, the matrix Eq. 54 becomes

[TABLE]

which, again, can be solved by successively solving $n_{j}/2$ systems with $I_{2}\otimes M+\Delta{}t{}\,C_{j}\otimes A$. Hence, the conditioning in this case is determined by the matrices in Eq. 52.

The analysis of Section 4 can be used to estimate the condition number of $M+\mu_{j}\Delta{}t{}\,A$. Estimate Eq. 53 for the second matrix in Eq. 52 is obtained similarly as for Theorem 5.1, except that we now have $\sigma_{\max}(C_{j})=\sqrt{\alpha_{j}^{2}+\beta_{j}^{2}}$ and $\lambda_{\min}((C_{j}+C_{j}^{T})/2)=\alpha_{j}$. ∎

Note that Remark 5.1 applies to the current situation as well. However, estimate Eq. 53 can now be used for all methods in Table 1 without any restriction on $\Delta{}t$.

6 Numerical examples

In the following examples we consider the IBVP Eq. 1 in 2d ($d=2$) with $\Omega={(0,1)}^{2}$ and homogeneous Dirichlet boundary conditions. For the time integration we choose Radau IA method of order 3 with fixed time steps $\Delta{}t{}\,=10^{-1},\;10^{-3},\;10^{-5}$, and $\Delta{}t{}\,=N^{-\frac{1}{d}}\sim h$, where $h$ is the average mesh size. The Radau IA coefficient matrix is given by

[TABLE]

Example 6.1**.**

To compare the condition number $\tilde{\kappa}(I_{s}\otimes M+\Delta{}t{}\,\Gamma\otimes A)$ of Radau IA with the condition number $\kappa(M+\Delta{}t{}\,A)$ of the implicit Euler method we consider two diffusion matrices: isotropic

[TABLE]

and anisotropic

[TABLE]

with quasi-uniform (Fig. 1(a)) and $\mathbb{D}^{-1}$-uniform (Fig. 1(b)) meshes obtained using the mesh generator bamg [11]. A quasi-uniform mesh can be also seen as $\mathbb{D}^{-1}$-uniform for $\mathbb{D}=I$.

We expect the conditioning of Radau IA to be similar to that of the implicit Euler (see Eqs. 46 and 5.1). This is verified in Fig. 2 (left column), which shows that the conditioning for the Radau IA method has the same general behaviour as that for the implicit Euler scheme. With increasing number of mesh elements, the order of conditioning is $\mathcal{O}(N^{\frac{2}{d}})$ for a fixed $\Delta{}t{}$ and becomes $\mathcal{O}(N^{\frac{1}{d}})$ for $\Delta{}t{}\,\sim h\sim N^{-\frac{1}{d}}$ for both methods, which is in perfect agreement with the theoretical prediction in Remark 4.2.

Adapting towards the diffusion matrix improves the conditioning of the stiffness matrix [18]. This effect is observed in our test case as well: for anisotropic diffusion, the conditioning for diffusion-adapted meshes (Fig. 2(c)) is noticeably smaller than that with quasi-uniform meshes (Fig. 2(b)), especially as $N$ is getting larger.

To compare the exact values for the conditioning of Radau IA with the bound Eq. 46 in Theorem 5.1, we use the exact eigenvalues for the matrices $P^{-1/2}MP^{-1/2}$ and $P^{-1/2}AP^{-1/2}$. Figure 2 (right column) shows that the bound Eq. 47 is greater than the exact condition number but exhibits essentially the same behaviour as $N$ increases.

Example 6.2**.**

To illustrate the effect of the diagonal scaling, we consider $\mathbb{D}=I$ (Laplace operator), $P=M_{D}+\Delta{}t{}\,A_{D}$, $M_{D}$ and $M_{lump}$, and adaptive anisotropic meshes for the interpolation of the function

[TABLE]

generated with bamg [11] using an anisotropic adaptive metric [13]. This function simulates the interaction between a boundary layer along the $x$-axis and a shock wave along the line $x=y+0.5$ (see Fig. 1(c) for a mesh example).

Consistently with Remark 4.1, Fig. 3(a) shows that the diagonal scaling with $P=M_{D}+\Delta{}t{}\,A_{D}$ (Euler) and $I_{s}\otimes M_{D}+\Delta{}t{}\,\Gamma\otimes A_{D}$ (Radau IA) reduces the effects caused by the mesh nonuniformity: the conditioning of both Euler and Radau IA methods is reduced by a factor ranging from 3 ($\Delta{}t=10^{-1}$) to 8 ($\Delta{}t=10^{-3},\;10^{-5}$).

For the scaling with $P=M_{D}$, Fig. 3(b) shows that the conditioning is getting worse for large $\Delta{}t{}$ ($\Delta{}t=10^{-1}$) but is improving for small $\Delta{}t{}$ ($\Delta{}t=10^{-3}$, $\Delta{}t=10^{-5}$). To explain this, we recall that the diagonal entries of $M$ are proportional to the corresponding patch volumes (see Lemma 3.1). Hence, scaling with $\mathbb{M}_{D}$ improves the conditioning issues caused by the mesh volume-nonuniformity (the same applies to $M$ and $M_{lump}$). On the other hand, $\max_{j}A_{jj}\leq\lambda_{\max}(A)\leq(d+1)\max_{j}A_{jj}$ [18, Lemma 4.1] and Lemma 3.6 imply that the largest eigenvalue of $A$ depends on the element shape (see also [18, section 4.1] and [20, section 3]). Thus, $\kappa(M^{-1}A)$ is not necessarily smaller than $\kappa(A)$ since rescaling of $A$ with respect to the patch volumes does not necessarily improve the conditioning issues caused by the shape. Therefore, in general, the overall effect depends on the magnitude of $\Delta{}t{}$: $\kappa(I+\Delta{}t{}M^{-1}A)$ is not necessarily better than $\kappa(M+\Delta{}t{}A)$ for large $\Delta{}t{}$, but we can expect that $\kappa(I+\Delta{}t{}M^{-1}A)<\kappa(M+\Delta{}t{}A)$ for small $\Delta{}t{}$.

7 Conclusions

For the FE equations, the conditioning of implicit RK methods with positive semidefinite symmetric part of the coefficient matrix is comparable to the implicit Euler integration (Theorem 5.1).

If the smallest eigenvalue of the symmetric part of the coefficient matrix is negative, Eq. 45 gives an upper bound on possible $\Delta{}t$ so that the system matrix $I_{s}\otimes M+\Delta{}t{}\,\Gamma\otimes A$ is positive definite. This condition implies $\Delta{}t\leq C\frac{\lambda_{\min}(B)}{\lambda_{\max}(A)}$ and can lead to a condition $\Delta t=\mathcal{O}(h^{2})$.

For the successive solution procedure, the conditioning of the system matrix $I_{s}\otimes M+\Delta{}t{}\,\Gamma\otimes A$ of an implicit RK integration is determined by two types of smaller matrices (Theorem 5.2)

[TABLE]

The first matrix is similar to the implicit Euler method and corresponds to the real eigenvalues $\mu_{j}$ of the RK matrix $\Gamma$ while the second corresponds to the complex eigenvalues $\alpha_{j}\pm i\beta_{j}$ of $\Gamma$.

There are three mesh-dependent factors that affect the conditioning of general implicit RK methods: the number of mesh elements (average mesh size), the mesh nonuniformity (in the Euclidean metric), and the mesh nonuniformity with respect to the inverse of the diffusion matrix. Preconditioning by the diagonal part of the system matrix itself or the mass matrix or by the lumped mass matrix reduces the effects of the mesh nonuniformity in the Euclidean metric (Theorems 4.2 and 4.3). These results are consistent with previous studies for the boundary value problems [1, 2, 9, 18] and explicit integration of linear diffusion problems [15, 16].

Similar result can be shown for higher order FE approximations. The major result is based on the preliminary estimates for the mass and stiffness matrices (Lemmas 3.1, 3.2, 3.3, 3.4 and 3.7) which, in their essential part, hold for higher order as well (with different constants depending of the order). For the FE, a symmetric positive definite (possibly lumped) mass matrix is equivalent to a diagonal matrix of node patch volumes [15, Lemmas 2.4] so that estimates equivalent to Lemmas 3.1 and 3.2 hold for higher order as well (with different constants). For the stiffness matrix, estimates equivalent to Lemmas 3.3 and 3.4 are given in [15, Lemmas 2.2 and 2.3] and the derivation of the bound on $\lambda_{\min}(A)$ in Lemma 3.7 does not depend on the FE order (up to a constant) and therefore holds also for higher order FE.