Stability of explicit Runge-Kutta methods for high order finite element approximation of linear parabolic equations

TL;DR

This paper analyzes the stability of explicit Runge-Kutta methods in high order finite element approximations of linear parabolic equations, providing bounds on the system matrix eigenvalues that influence time step size.

Contribution

It introduces new bounds on the eigenvalues of the system matrix, linking stability to mesh geometry, basis functions, and diffusion matrix properties.

Findings

Bounds on eigenvalues depend on matrix entry ratios and mesh geometry.

The tightness of bounds varies with problem parameters and mesh configuration.

Results offer insights into stability constraints for explicit schemes on nonuniform meshes.

Abstract

We study the stability of explicit Runge-Kutta methods for high order Lagrangian finite element approximation of linear parabolic equations and establish bounds on the largest eigenvalue of the system matrix which determines the largest permissible time step. A bound expressed in terms of the ratio of the diagonal entries of the stiffness and mass matrices is shown to be tight within a small factor which depends only on the dimension and the choice of the reference element and basis functions but is independent of the mesh or the coefficients of the initial-boundary value problem under consideration. Another bound, which is less tight and expressed in terms of mesh geometry, depends only on the number of mesh elements and the alignment of the mesh with the diffusion matrix. The results provide an insight into how the interplay between the mesh geometry and the diffusion matrix affects the stability of explicit integration schemes when applied to a high order finite element approximation of linear parabolic equations on general nonuniform meshes.

Equations

$$\{\begin{array}{l}\begin{array}{rlrlr}& u_t=\nabla \cdot \left(\mathbb{D}\nabla u\right),& & \mathbf{x}\in \Omega ,& t\in \left(0,T\right],\\ & u(\mathbf{x},t)=0,& & \mathbf{x}\in {\Gamma }_D,& t\in \left(0,T\right],\\ & \mathbb{D}\nabla u(\mathbf{x},t)\cdot \mathbf{n}=0,& & \mathbf{x}\in {\Gamma }_N,& t\in \left(0,T\right],\\ & u(\mathbf{x},0)=u^0(\mathbf{x}),& & \mathbf{x}\in \Omega ,& \end{array}\end{array}$$

$$\{\begin{array}{l}\begin{array}{rlrl}{\Vert u(\cdot ,t)\Vert }_{L^2(\Omega )}& \le {\Vert u^0\Vert }_{L^2(\Omega )},& & t\in \left(0,T\right],\\ |||u(\cdot ,t)|||& \le |||u^0|||,& & t\in \left(0,T\right],\end{array}\end{array}$$

$${\int }_{\Omega }{\partial }_t{u}^h{v}^h\mathrm{d}\mathbf{x}=-{\int }_{\Omega }\nabla v^h\cdot \mathbb{D}\nabla u^h\mathrm{d}\mathbf{x},\forall v^h\in V_D^h,$$

$${\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_D^h.$$

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

$$M{\mathbf{U}}_t=-A\mathbf{U},\mathbf{U}(0)={\mathbf{U}}_0,$$

$$M_{ij}={\int }_{\Omega }{\varphi }_i{\varphi }_j\mathrm{d}\mathbf{x}\text{and}{A}_{ij}={\int }_{\Omega }\nabla {\varphi }_i\cdot \mathbb{D}\nabla {\varphi }_j\mathrm{d}\mathbf{x}$$

$${\Vert u_n^h\Vert }_{L^2(\Omega )}\le \sqrt{\frac{c_2}{c_1}}{\Vert u_0^h\Vert }_{L^2(\Omega )}\text{and}|||u_n^h|||\le |||u_0^h|||,$$

$$\underset{i}{\max }\left|R\left(-\tau {\lambda }_i\left({\tilde{M}}^{-1}A\right)\right)\right|\le 1.$$

$$A\le \eta A_D.$$

$${\mathbf{u}}^{T}S\mathbf{v}+{\mathbf{v}}^{T}S\mathbf{u}\le {\mathbf{u}}^{T}S\mathbf{u}+{\mathbf{v}}^{T}S\mathbf{v}.$$

$${\mathbf{u}}^{T}A\mathbf{u}$$

$$\le \sum _i\eta {\int }_{\Omega }{\left(u_i\nabla {\varphi }_i\right)}^T\mathbb{D}\left(u_i\nabla {\varphi }_i\right)\mathrm{d}\mathbf{x}$$

$$=\eta \sum _i{u}_i^2{\int }_{\Omega }\nabla {\varphi }_i^T\mathbb{D}\nabla {\varphi }_i\mathrm{d}\mathbf{x}$$

$$={\mathbf{u}}^T\eta A_D\mathbf{u}.\text{qed}$$

$$C_{H^1}=\underset{i}{\max }{|{\hat{\varphi }}_i|}_{H^1(\hat{K})}^2.$$

$$A_{ii}\le C_{H^1}\sum _{K\in {\omega }_i}\left|K\right|\underset{\mathbf{x}\in K}{\max }{\Vert {(F_K^{\prime })}^{-1}\mathbb{D}{(F_K^{\prime })}^{-T}\Vert }_2.$$

$$A_{ii}={\int }_{\Omega }\nabla {\varphi }_i^T\mathbb{D}\nabla {\varphi }_i\mathrm{d}\mathbf{x}=\sum _{K\in {\omega }_i}{\int }_K\nabla {\varphi }_i^T\mathbb{D}\nabla {\varphi }_i\mathrm{d}\mathbf{x}.$$

$$A_{ii}$$

$$\le \sum _{K\in {\omega }_i}\left|K\right|{\Vert \hat{\nabla }{\hat{\varphi }}_i\Vert }_{L^2(\hat{K})}^2\underset{\mathbf{x}\in K}{\max }{\Vert {(F_K^{\prime })}^{-1}\mathbb{D}{(F_K^{\prime })}^{-T}\Vert }_2$$

$$\le C_{H^1}\sum _{K\in {\omega }_i}\left|K\right|\underset{\mathbf{x}\in K}{\max }{\Vert {(F_K^{\prime })}^{-1}\mathbb{D}{(F_K^{\prime })}^{-T}\Vert }_2.\text{qed}$$

$$\mathcal{W}=\mathrm{diag}\left(\left|{\omega }_1\right|,\dots ,\left|{\omega }_{N_{\varphi }}\right|\right).$$

$${\hat{\lambda }}_{\tilde{M}}\mathcal{W}\le \tilde{M}\le {\hat{\Lambda }}_{\tilde{M}}\mathcal{W}.$$

$${\mathbf{u}}^T\tilde{M}\mathbf{u}$$

$$\le \sum _K\left|K\right|{\hat{\Lambda }}_{\tilde{M}}{\Vert {\mathbf{u}}_K\Vert }_2^2={\hat{\Lambda }}_{\tilde{M}}\sum _i{u}_i^2\sum _{K\in {\omega }_i}\left|K\right|$$

$$={\hat{\Lambda }}_{\tilde{M}}\sum _i{u}_i^2\left|{\omega }_i\right|={\hat{\Lambda }}_{\tilde{M}}{\mathbf{u}}^T\mathcal{W}\mathbf{u}.$$

$$\frac{{\hat{\lambda }}_{{\tilde{M}}_1}}{{\hat{\Lambda }}_{{\tilde{M}}_2}}{\tilde{M}}_2\le {\tilde{M}}_1\le \frac{{\hat{\Lambda }}_{{\tilde{M}}_1}}{{\hat{\lambda }}_{{\tilde{M}}_2}}{\tilde{M}}_2.$$

$${\Vert u_n^h\Vert }_{L^2(\Omega )}\le \sqrt{\kappa (M_{\hat{K}})\kappa ({\tilde{M}}_{\hat{K}})}{\Vert u_0^h\Vert }_{L^2(\Omega )}$$

$$|||u_n^h|||\le |||u_0^h|||.$$

$$\frac{1}{\kappa ({\tilde{M}}_{\hat{K}})}{\tilde{M}}_D\le \tilde{M}\le \kappa ({\tilde{M}}_{\hat{K}}){\tilde{M}}_D.$$

$${\hat{\lambda }}_{\tilde{M}}{\mathcal{W}}_{ii}\le {\tilde{M}}_{ii}\le {\hat{\Lambda }}_{\tilde{M}}{\mathcal{W}}_{ii},$$

$$u_i{\hat{\lambda }}_{\tilde{M}}{\mathcal{W}}_{ii}{u}_i\le u_i{\tilde{M}}_{ii}{u}_i\le u_i{\hat{\Lambda }}_{\tilde{M}}\mathcal{W}{u}_i\text{for any ui}.$$

$${\hat{\lambda }}_{\tilde{M}}\mathcal{W}\le {\tilde{M}}_D\le {\hat{\Lambda }}_{\tilde{M}}\mathcal{W}.$$

$$\max_{i}\frac{A_{ii}}{\tilde{M}_{ii}}\leq\lambda_{\max}\bigl{(}\tilde{M}^{-1}A\bigr{)}\leq\eta\,\kappa(\tilde{M}_{\hat{K}})\max_{i}\frac{A_{ii}}{\tilde{M}_{ii}},$$

Full text

This is a preprint of a contribution to A. Abdulle et al. (eds.), Numerical Mathematics and Advanced Applications — ENUMATH 2013, Lecture Notes in Computational Science and Engineering, vol. 103.

© Springer International Publishing Switzerland 2015. The final version is available at https://doi.org/10.1007/978-3-319-10705-9_16.

††journal: Numerical Mathematics and Advanced Applications — ENUMATH 2013

1 Introduction

We consider the initial-boundary value problem (IBVP)

[TABLE]

where $\Omega\subset\mathbb{R}^{d}$ ($d\geq 1$) is a bounded polygonal or polyhedral domain, $\Gamma_{D}\cup\Gamma_{N}=\partial\Omega$, $\operatorname{meas}_{d-1}\Gamma_{D}>0$, $u^{0}$ is a given function, and $\mathbb{D}=\mathbb{D}(\bm{x})$ is the diffusion matrix, which is assumed to be time-independent, symmetric and uniformly positive-definite on $\Omega$. If u^{0}\in H^{1}_{D}(\Omega)=\set{v\in H^{1}(\Omega):\text{v=0$on$\Gamma_{D}}} and $u$ is sufficiently smooth, then the solution of the IBVP satisfies the stability estimates

[TABLE]

where ${|||u|||}={\left\lVert\mathbb{D}^{1/2}\nabla u\right\rVert}_{L^{2}(\Omega)}$ is the energy norm. We are interested in the stability conditions so that the numerical approximation preserves these stability estimates.

The stability of explicit Runge-Kutta methods depends on the largest eigenvalue of the corresponding system matrix, which, in turn, depends on the mesh and the coefficients of the IBVP. For our model problem this means that we need to estimate the largest eigenvalue of $M^{-1}A$, where $M$ and $A$ are the mass and stiffness matrices for the finite element discretization of the IBVP 1 [4, Theorem 3.1]. For the Laplace operator on a uniform mesh it is well known that $\lambda_{\max}(M^{-1}A)\sim N^{2/d}$, where $N$ is the number of mesh elements. For general meshes and diffusion coefficients, estimates have been derived recently in Huang et al. [4] and Zhu and Du [6, 7] (see also [1, 2, 3, 5] for estimates on $M$ and $A$). All of these works allow anisotropic diffusion coefficients and anisotropic meshes, while the former employs a more accurate measure for the interplay between the mesh geometry and the diffusion matrix and gives a sharper estimate on $\lambda_{\max}(M^{-1}A)$ than the latter. On the other hand, [4] considers only linear finite elements whereas the estimates in [7] are valid for both linear and higher order finite elements.

The purpose of this paper is to extend the result of [4] to high order Lagrangian finite elements as well as provide a mathematical understanding of how the interplay between the mesh geometry and the diffusion matrix affects the stability condition. We show that the main result of [4, Theorem 3.3] holds for high order finite elements as well. The analysis is based on bounds on the mass and stiffness matrices. We follow the approach in [4, 5] and derive simple but accurate bounds for the case of high order Lagrangian finite elements on simplicial meshes (Lemmas 2, 3, 4 and 5). We also consider the more general case of surrogate mass matrices $\tilde{M}$. The main result (Theorem 1) shows that $\lambda_{\max}(\tilde{M}^{-1}A)$ is proportional to the maximum ratio between the corresponding diagonal entries of the stiffness and surrogate mass matrices. Moreover, $\lambda_{\max}(\tilde{M}^{-1}A)$ is bounded by a term depending only on the number of the mesh elements and the alignment of the shape of the mesh elements with the inverse of the diffusion matrix.

2 Stability condition for explicit time stepping

Let $\set{\mathcal{T}_{h}}$ be a family of simplicial meshes for $\Omega$ and $V^{h}$ the Lagrangian $\mathbb{P}_{m}$ ($m\geq 1$) finite element space associated with $\mathcal{T}_{h}$. Let $K$ be an arbitrary element of $\mathcal{T}_{h}$, $\hat{K}$ the reference element, and $\omega_{i}$ the element patch of the $i^{\text{th}}$ vertex (Fig. 1); element and patch volumes are denoted by ${\left\lvert K\right\rvert}$ and ${\left\lvert\omega_{i}\right\rvert}=\sum_{K\in\omega_{i}}{\left\lvert K\right\rvert}$. For each $K\in\mathcal{T}_{h}$ let $F_{K}\colon\hat{K}\to K$ be an invertible affine mapping and $F_{K}^{\prime}$ its Jacobian matrix which is constant and satisfies $\det(F_{K}^{\prime})={\left\lvert K\right\rvert}$ (for simplicity, we assume that ${\lvert\hat{K}\rvert}=1$). We further assume that the mesh is fixed for all time steps.

With $V^{h}_{D}=V^{h}\cap H^{1}_{D}(\Omega)$, the finite element solution $u^{h}(t)\in V^{h}_{D}$, $t\in\left(0,T\right]$, is defined by

[TABLE]

subject to the initial condition

[TABLE]

Let $N_{\phi}$ be the dimension of the finite element space $V^{h}_{D}$ and denote a nodal basis of $V^{h}_{D}$ by $\set{\phi_{1},\dotsc,\phi_{N_{\phi}}}$, then $u^{h}$ can be expressed as

[TABLE]

Using $\bm{U}={(u^{h}_{1},\dotsc,u^{h}_{N_{\phi}})}^{T}$, 2 and 3 can be written into a matrix form

[TABLE]

where the mass and stiffness matrices $M$ and $A$ are defined by

[TABLE]

for all $i,j=1,\dotsc,N_{\phi}$. We further assume that surrogate mass matrices $\tilde{M}$ considered throughout the paper satisfy

1. (M1)

The reference element matrix $\tilde{M}_{\hat{K}}$ is symmetric positive definite. 2. (M2)

The element matrix $\tilde{M}_{K}$ satisfies $\tilde{M}_{K}=|K|\tilde{M}_{\hat{K}}$.

For example, (M1) and (M2) are satisfied for any mass lumping by means of numerical quadrature with positive weights.

Lemma 1** ([4, Theorem 3.1]).**

For a given explicit RK method with the polynomial stability function $R$ and a symmetric positive definite surrogate matrix $\tilde{M}$ that satisfies222In the following, the less-than-or-equal-to sign for matrices means that the difference between the right-hand side and left-hand side terms is positive semidefinite. $c_{1}\tilde{M}\leq M\leq c_{2}\tilde{M}$ for some positive constants $c_{1}$ and $c_{2}$, the finite element approximation $u^{h}_{n}$ at $t_{n}=n\tau$ satisfies

[TABLE]

if the time step $\tau$ is chosen such that

[TABLE]

This lemma is proven in [4] for the linear finite element discretization. However, from the proof one can see that it is valid for any system in the form of 4 with symmetric positive definite matrices $M$ and $A$. Particularly, it can be used for the system 4 resulting from the $\mathbb{P}_{m}$ finite element discretization. In the following, we establish a series of lemmas for bounds on the stiffness and mass matrices $A$ and $\tilde{M}$ and then develop bounds for $\lambda_{\max}(\tilde{M}^{-1}A)$.

Lemma 2**.**

Let $\eta$ be the maximal number of basis functions per element. Then the stiffness matrix $A$ and its diagonal part $A_{D}$ for $\mathbb{P}_{m}$ finite elements satisfy

[TABLE]

Proof.

Notice that for any positive semi-definite matrix $S$ and any vectors $\bm{u}$ and $\bm{v}$ we have

[TABLE]

From this,

[TABLE]

Lemma 3**.**

Let $\hat{\phi}_{i}$ be the basis functions on the reference element that correspond to $\phi_{i}$ and

[TABLE]

Then the diagonal entries $A_{ii}$ of the stiffness matrix $A$ are bounded by

[TABLE]

Proof.

From the definition of the stiffness matrix we have

[TABLE]

Let $\hat{\nabla}={\partial}/{\partial\bm{\xi}}$ be the gradient operator in $\hat{K}$. The chain rule yields $\nabla={(F_{K}^{\prime})}^{-T}\hat{\nabla}$ and together with $\det(F_{K}^{\prime})={\left\lvert K\right\rvert}$ we obtain

[TABLE]

Lemma 4**.**

Let $\tilde{M}$ be a surrogate $\mathbb{P}_{m}$ finite element mass matrix, $\hat{\Lambda}_{\tilde{M}}$ and $\hat{\lambda}_{\tilde{M}}$ be the largest and smallest eigenvalues of the surrogate mass matrix $\tilde{M}_{\hat{K}}$ on the reference element and

[TABLE]

Then

[TABLE]

Proof.

We have

[TABLE]

The lower bound can be obtained similarly. ∎

Lemma 5**.**

Let $\tilde{M}_{1}$ and $\tilde{M}_{2}$ be two surrogate mass matrices for $\mathbb{P}_{m}$ finite elements. Then

[TABLE]

Proof.

Use Lemma 4 by applying 5 to $\tilde{M}_{1}$ and $\tilde{M}_{2}$. ∎

Corollary 1**.**

Let $\kappa(M_{\hat{K}})$ and $\kappa(\tilde{M}_{\hat{K}})$ be the condition numbers of the full and the surrogate reference element mass matrices. Under the assumptions of Lemma 1 we have

[TABLE]

and

[TABLE]

Proof.

Use $\tilde{M}_{1}=M$ and $\tilde{M}_{2}=\tilde{M}$ in Lemma 5 and apply Lemma 1. ∎

Corollary 2**.**

The surrogate mass matrix $\tilde{M}$ for $\mathbb{P}_{m}$ finite elements and its diagonal part $\tilde{M}_{D}$ satisfy

[TABLE]

Proof.

Using 5 with the canonical basis vector $\bm{e}_{i}$ implies

[TABLE]

which gives

[TABLE]

Since $\tilde{M}_{D}$ and $\mathcal{W}$ are diagonal matrices, this leads to

[TABLE]

The statement now follows from Lemma 5 with $\tilde{M}_{1}=\tilde{M}$ and $\tilde{M}_{2}=\tilde{M}_{D}$. ∎

Having obtained the preliminary bounds on the stiffness and mass matrices $A$ and $\tilde{M}$, we can now give the estimate for the largest eigenvalue of the system matrix $\tilde{M}^{-1}A$ for $\mathbb{P}_{m}$ finite elements.

Theorem 1**.**

The eigenvalues of $\tilde{M}^{-1}A$ are real and positive and the largest eigenvalue is bounded by

[TABLE]

where $\eta$ is the maximal number of basis functions per element. Further,

[TABLE]

Proof.

Since $\tilde{M}$ and $A$ are symmetric positive definite, the eigenvalues of $\tilde{M}^{-1}A$ are real and positive. The lower bound in 7 is obtained by using the canonical basis vectors $\bm{e}_{i}$ and the upper bound follows from Lemmas 2 and 2,

[TABLE]

The geometric bound 8 is a direct consequence of Lemmas 2, 3 and 4,

[TABLE]

∎

Theorem 1 can be used in combination with Lemma 1 or Corollary 1 to derive the stability condition of a given explicit Runge-Kutta scheme, as shown in the next example.

Example 1** (Explicit Euler method).**

The stability region of the explicit Euler method includes the real interval $[-2,0]$. Lemma 1 implies that the method is stable if

[TABLE]

Using Theorem 1, we conclude that the method is stable if the time step $\tau$ satisfies

[TABLE]

or, in terms of mesh geometry,

[TABLE]

Remark 1**.**

Lemmas 2, 3 and 2 are very general and valid for any mesh, any $\mathbb{D}$ and any surrogate mass matrix $\tilde{M}$ satisfying (M1) and (M2). More accurate bounds can be obtained if more information is available about the mesh or the stiffness and mass matrices.

For example, if $A$ is an M-matrix, then the Gershgorin circle theorem yields $\lambda_{\max}(A)\leq 2\max_{i}A_{ii}$ [4, Remark 2.2] and therefore $\eta$ in Theorem 1 can be replaced by $2$.

If $\tilde{M}=M$ (no mass lumping), then, instead of estimating $M_{D}$ through 6, a direct calculation for the standard $\mathbb{P}_{m}$ finite elements yields

[TABLE]

and

[TABLE]

resulting in a slighly more accurate bound in Corollary 2. For simplicity, in Lemma 3 we used $C_{H^{1}}=\max_{i}{\lvert\hat{\phi}_{i}\rvert}_{H^{1}(\hat{K})}^{2}$. A slightly more accurate bound can be derived if we use

[TABLE]

3 Summary and conclusion

Theorem 1 states that the largest eigenvalue of the system matrix and, thus, the largest permissible time step can be bounded by a term depending only on the number of mesh elements and the alignment of the mesh with the diffusion matrix.

The bound in terms of matrix entries is tight within a small factor which depends only on the dimension and the choice of the reference element and basis functions but is independent of the mesh or the coefficients of the IBVP. This is valid for any Lagrangian $\mathbb{P}_{m}$ finite elements with $m\geq 1$.

A similar result is obtained by Zhu and Du [7, Theorem 3.1]. In our notation, it can be written as

[TABLE]

The significant difference between the new bound 8 and the bound 9 is the factor which represents the interplay between the mesh geometry and the diffusion matrix,

[TABLE]

vs.

[TABLE]

For isotropic $\mathbb{D}$ or isotropic meshes both terms are comparable. However, the former is smaller than the latter in general. In particular, if both $\mathbb{D}$ and $K$ are anisotropic, then the difference between 8 and 9 can be very significant (see [4, Sect. 4.4] for a numerical example in case of $\mathbb{P}_{1}$ finite elements). In this sense, Theorem 1 can be seen either as a generalization of [4] to $\mathbb{P}_{m}$ ($m\geq 2$) finite elements or as a more accurate version of [7] for anisotropic meshes and general diffusion coefficients.

Finally, we would like to point out that a similar result can be established for $p$-adaptive finite elements without major modifications.