An analysis of discontinuous Galerkin methods for the compressible Euler equations: entropy and $L_2$ stability

Abstract

The objective of this article is to characterize the entropy and $L_2$ stability of several representative discontinuous Galerkin (DG) methods for solving the compressible Euler equations. Towards this end, three DG methods are constructed: one DG method with entropy variables as its unknowns, and two DG methods with conservative variables as their unknowns. These methods are employed in order to discretize the compressible Euler equations in space. Thereafter, the resulting semi-discrete formulations are analyzed, and the entropy and $L_2$ stability characteristics are evaluated. It is shown that the semi-discrete formulation of the DG method with entropy variables is entropy and $L_2$ stable. Furthermore, it is shown that the semi-discrete formulations of the DG methods with conservative variables are only guaranteed to be entropy and $L_2$ stable under the following assumptions: the entropy projection errors vanish, or the terms containing the entropy projection errors are non-positive. Thereafter, the semi-discrete formulation with entropy variables, and one of the semi-discrete formulations with conservative variables, are discretized in time with an ‘algebraically stable’ Runge–Kutta (RK) scheme. The resulting formulations are fully-discrete and can be immediately applied to practical problems. In this article, they are employed to simulate a vortex propagating for long distances. It is shown that temporal stability is maintained by the DG method with entropy variables, but the DG method with conservative variables exhibits instability.

Appendices

Appendix A: Notational conventions

The principal results in this paper are expressed via a mixture of index notation and vector notation. This combination of notation is best explained with an example. Consider the term $\left( \varvec{w}^{h}_{,\varvec{x}_i} \right) ^T \varvec{f}^{i} \left( \varvec{u}^h \right)$. This can be expanded as follows when $d =2$

$$\begin{aligned} \left( \varvec{w}^{h}_{,\varvec{x}_i} \right) ^T \varvec{f}^{i} \left( \varvec{u}^h \right) = \left( \frac{\partial \varvec{w}^{h}}{\partial \varvec{x}_1} \right) ^T \varvec{f}^{1} \left( \varvec{u}^h \right) + \left( \frac{\partial \varvec{w}^{h}}{\partial \varvec{x}_2} \right) ^T \varvec{f}^{2} \left( \varvec{u}^h \right) , \end{aligned}$$

where the repeated index i facilitates the standard Einstein summation over d dimensions, and the transpose facilitates the standard dot product between a pair of m-vectors [for instance $\varvec{w}^{h}_{,\varvec{x}_1}$ and $\varvec{f}^{1} \left( \varvec{u}^h \right)$].

Appendix B: Supporting lemmas

Lemma B.1

Suppose that the pressure p and density $\rho$ are bounded in the following fashion

$$\begin{aligned} 0< p \le M, \quad 0 < \rho \le N, \end{aligned}$$

(B.1)

where M and N are positive numbers. Then, the symmetric Jacobian matrices ${\widetilde{A}}_0$ and ${\widetilde{A}}_{0}^{-1}\ ($defined in Eq. (3.7)) are guaranteed to be positive definite (PD).

Proof

In what follows, we present the proof for the 2D case. The proof for the 3D case and for higher dimensions, is very similar. In 2D, one may begin by examining the matrix ${\widetilde{A}}_0$, which can be expressed as follows

$$\begin{aligned} {\widetilde{A}}_0 = \begin{bmatrix} \rho&\quad \rho u&\quad \rho v&\quad \rho H-p \\ \rho u&\quad \rho u^2 + p&\quad \rho u v&\quad \rho u H \\ \rho v&\quad \rho u v&\quad \rho v^2 + p&\quad \rho v H \\ \rho H - p&\quad \rho u H&\quad \rho v H&\quad \rho H^2 - \gamma e p \end{bmatrix}, \end{aligned}$$

where

$$\begin{aligned} H = \gamma e + \frac{1}{2} \left( u^2 + v^2 \right) . \end{aligned}$$

The characteristic polynomial of ${\widetilde{A}}_0$ takes the following form

$$\begin{aligned} \frac{1}{4} \left( \lambda - p \right) \left( a \lambda ^3 + b \lambda ^2 + c \lambda +d \right) = 0, \end{aligned}$$

(B.2)

where

$$\begin{aligned} a&= 4, \nonumber \\ b&= - \rho \left( \left( 2 + u^2 + v^2 \right) ^2 + 4 e \left( \gamma \left( 1 + u^2 + v^2 + e \right) - 1 \right) \right) , \nonumber \\ c&= \rho p \left( \left( 2 + u^2 + v^2 \right) ^2 + 4e \left( \gamma e + 1 \right) \right) , \nonumber \\ d&= -4 \rho e p^2. \end{aligned}$$

(B.3)

It immediately follows from Eq. (B.2), that the pressure p is an eigenvalue. Evidently, this eigenvalue is positive since Eq. (B.1) holds.

One may determine the signs of the remaining eigenvalues of ${\widetilde{A}}_0$ by first observing that the matrix is symmetric and will have all real entries if Eq. (B.1) is satisfied. As a result, the eigenvalues of ${\widetilde{A}}_0$ will be real, i.e. the roots of its characteristic polynomial will be real. In accordance with Decartes’ rule of signs, a cubic polynomial with real roots, and with coefficients $a > 0$, $b < 0$, $c > 0$, and $d < 0$ is guaranteed to have three positive roots. By inspection of Eq. (B.3), it immediately follows that the coefficients a, b, c, and d will have the required signs if Eq. (B.1) holds, if $\gamma > 1$ (which is true for virtually all gases), and if one notes that $e = p /(\rho \left( \gamma -1 \right) ) > 0$. This completes the proof that ${\widetilde{A}}_0$ has positive eigenvalues.

The eigenvalues of ${\widetilde{A}}_{0}^{-1}$ are obtained by taking the reciprocals of the eigenvalues of ${\widetilde{A}}_{0}$. The eigenvalues are positive as long as Eq. (B.1) holds, as under these circumstances, the eigenvalues of ${\widetilde{A}}_0$ remain positive, and the reciprocals of positive numbers are positive. This completes the proof that ${\widetilde{A}}_{0}^{-1}$ has positive eigenvalues. $\square$

Lemma B.2

The jump in the entropy functional $F^i \left( \varvec{v}\right)$ across an interface is governed by the following

$$\begin{aligned}&\left[ \!\left[ F^i \left( \varvec{u}^h \right) \right] \! \right] ^{-}_{+} + \left( \left\{ \!\left\{ \varvec{v}\left( \varvec{u}^h \right) \right\} \! \right\} ^{-}_{+}\right) ^T \left[ \!\left[ \varvec{f}^i \left( \varvec{u}^h \right) \right] \! \right] ^{+}_{-} \nonumber \\&\quad = \frac{1}{2} \int _{0}^{1} \left( 1 - \theta \right) \left( \left[ \!\left[ \varvec{v}\left( \varvec{u}^h \right) \right] \! \right] ^{+}_{-}\right) ^T \left( {\widetilde{A}}_i \left( \overline{\overline{\varvec{v}}} \left( \theta \right) \right) - {\widetilde{A}}_i \left( \overline{\varvec{v}} \left( \theta \right) \right) \right) \left[ \!\left[ \varvec{v}\left( \varvec{u}^h \right) \right] \! \right] ^{+}_{-}d \theta .\nonumber \\ \end{aligned}$$

(B.4)

Proof

The proof is virtually identical to the proofs of similar statements in [69] and [5]. It is repeated here for the sake of completeness.

Recall that ${\mathcal {F}}^{i} = {\mathcal {F}}^i \left( \varvec{v}\left( \varvec{u}\right) \right)$. On utilizing this fact in conjunction with Taylor’s theorem, one can obtain the following

$$\begin{aligned}&{\mathcal {F}}^i \left( \varvec{v}\left( \varvec{u}^h_{+} \right) \right) - {\mathcal {F}}^i \left( \varvec{v}\left( \varvec{u}_{-}^h \right) \right) - {\mathcal {F}}^{i}_{, \varvec{v}} \left( \varvec{v}\left( \varvec{u}^h_{+} \right) \right) \left( \varvec{v}\left( \varvec{u}^h_{+} \right) - \varvec{v}\left( \varvec{u}_{-}^h \right) \right) \nonumber \\&\quad + \int _{0}^{1} \left( 1- \theta \right) \left( \varvec{v}\left( \varvec{u}^h_{+} \right) - \varvec{v}\left( \varvec{u}_{-}^h \right) \right) ^T {\mathcal {F}}^{i}_{,\varvec{v}, \varvec{v}} \left( \overline{\varvec{v}} \left( \theta \right) \right) \left( \varvec{v}\left( \varvec{u}^h_{+} \right) - \varvec{v}\left( \varvec{u}_{-}^h \right) \right) d \theta = 0,\nonumber \\ \end{aligned}$$

(B.5)

$$\begin{aligned}&{\mathcal {F}}^i \left( \varvec{v}\left( \varvec{u}^h_{+} \right) \right) - {\mathcal {F}}^i \left( \varvec{v}\left( \varvec{u}^h_{-} \right) \right) - {\mathcal {F}}^{i}_{, \varvec{v}} \left( \varvec{v}\left( \varvec{u}^h_{-} \right) \right) \left( \varvec{v}\left( \varvec{u}^h_{+} \right) -\varvec{v}\left( \varvec{u}^h_{-} \right) \right) \nonumber \\&\quad - \int _{0}^{1} \left( 1- \theta \right) \left( \varvec{v}\left( \varvec{u}^h_{+} \right) - \varvec{v}\left( \varvec{u}^h_{-} \right) \right) ^T {\mathcal {F}}^{i}_{,\varvec{v}, \varvec{v}} \left( \overline{\overline{\varvec{v}}} \left( \theta \right) \right) \left( \varvec{v}\left( \varvec{u}^h_{+} \right) - \varvec{v}\left( \varvec{u}^h_{-} \right) \right) d \theta = 0. \nonumber \\ \end{aligned}$$

(B.6)

Upon multiplying Eqs. (B.5) and (B.6) by (1/2) and summing the results, one obtains

$$\begin{aligned}&{\mathcal {F}}^i \left( \varvec{v}\left( \varvec{u}^h_{+} \right) \right) - {\mathcal {F}}^i \left( \varvec{v}\left( \varvec{u}^h_{-} \right) \right) - \frac{1}{2} \left( \varvec{f}^i \left( \varvec{u}^h_{+} \right) + \varvec{f}^i \left( \varvec{u}^h_{-} \right) \right) ^T \left( \varvec{v}\left( \varvec{u}^h_{+} \right) - \varvec{v}\left( \varvec{u}^h_{-} \right) \right) \nonumber \\&\quad = \frac{1}{2} \int _{0}^{1} \left( 1- \theta \right) \left( \varvec{v}\left( \varvec{u}^h_{+} \right) - \varvec{v}\left( \varvec{u}^h_{-} \right) \right) ^T \left( {\widetilde{A}}_i \left( \overline{\overline{\varvec{v}}} \left( \theta \right) \right) - {\widetilde{A}}_i \left( \overline{\varvec{v}} \left( \theta \right) \right) \right) \left( \varvec{v}\left( \varvec{u}^h_{+} \right) - \varvec{v}\left( \varvec{u}^h_{-} \right) \right) d \theta ,\nonumber \\ \end{aligned}$$

(B.7)

where the fact that $\left( \varvec{f}^i \right) ^T = {\mathcal {F}}^{i}_{, \varvec{v}}$ and ${\widetilde{A}}_{i} = {\mathcal {F}}^{i}_{,\varvec{v}, \varvec{v}}$ has been used (cf. Eqs. (3.3) and (3.7)). Setting Eq. (B.7) aside for the moment, consider the following jump identity that derives from Eq. (3.5)

$$\begin{aligned}&F^i \left( \varvec{u}^h_{+} \right) - F^i \left( \varvec{u}^h_{-} \right) + {\mathcal {F}}^i \left( \varvec{v}\left( \varvec{u}^h_{+} \right) \right) - {\mathcal {F}}^i \left( \varvec{v}\left( \varvec{u}^h_{-} \right) \right) \nonumber \\&\quad = \frac{1}{2} \left( \varvec{v}\left( \varvec{u}^h_{+} \right) + \varvec{v}\left( \varvec{u}^h_{-} \right) \right) ^T \left( \varvec{f}^i \left( \varvec{u}^h_{+} \right) - \varvec{f}^i \left( \varvec{u}^h_{-} \right) \right) \nonumber \\&\qquad + \frac{1}{2} \left( \varvec{f}^i \left( \varvec{u}^h_{+} \right) + \varvec{f}^i \left( \varvec{u}^h_{-} \right) \right) ^T \left( \varvec{v}\left( \varvec{u}^h_{+} \right) - \varvec{v}\left( \varvec{u}^h_{-} \right) \right) . \end{aligned}$$

(B.8)

Substituting Eq. (B.7) into Eq. (B.8) yields

$$\begin{aligned}&F^i \left( \varvec{u}^h_{-} \right) - F^i \left( \varvec{u}^h_{+} \right) + \frac{1}{2} \left( \varvec{v}\left( \varvec{u}^h_{-} \right) + \varvec{v}\left( \varvec{u}^h_{+} \right) \right) ^T \left( \varvec{f}^i \left( \varvec{u}^h_{+} \right) - \varvec{f}^i \left( \varvec{u}^h_{-} \right) \right) \nonumber \\&\quad = \frac{1}{2} \int _{0}^{1} \left( 1- \theta \right) \left( \varvec{v}\left( \varvec{u}^h_{+} \right) - \varvec{v}\left( \varvec{u}^h_{-} \right) \right) ^T \left( {\widetilde{A}}_i \left( \overline{\overline{\varvec{v}}} \left( \theta \right) \right) - {\widetilde{A}}_i \left( \overline{\varvec{v}} \left( \theta \right) \right) \right) \left( \varvec{v}\left( \varvec{u}^h_{+} \right) - \varvec{v}\left( \varvec{u}^h_{-} \right) \right) d \theta .\nonumber \\ \end{aligned}$$

(B.9)

Upon substituting the spatial jump [Eq. (5.1)] and average [Eq. (5.2)] operators into Eq. (B.9) one obtains Eq. (B.4). $\square$

Lemma B.3

Suppose that the initial condition $\varvec{u}\left( t_0 \right) \in \varvec{L}_2 \left( \varOmega \right)$, and the pressure p and density $\rho$ are positive and bounded for all convex combinations of states $\varvec{u}\left( t_0 \right)$ and $\varvec{u}^{*}$ defined by $\widehat{\widehat{\varvec{u}}} \left( \theta \right) = \varvec{u}^{*} + \theta \left( \varvec{u}\left( t_0 \right) - \varvec{u}^{*} \right)$, where $0 \le \theta \le 1$. Under these circumstances, the following statements hold: (i) the $L_2$ norm of $\varvec{u}\left( t_0 \right) - \varvec{u}^{*}$ is well-defined; (ii) the eigenvalues of ${\widetilde{A}}_0^{-1} \left( \widehat{\widehat{\varvec{u}}} \left( \theta \right) \right)$ are real and positive; and (iii) the broken integral of $H \left( \varvec{u}\left( t_0 \right) , \varvec{u}^{*} \right)$ over the domain ${\mathcal {T}}_h$ is bounded in the following fashion

$$\begin{aligned} \sum _{T_k \in {\mathcal {T}}_h} \int _{T_k} H \left( \varvec{u}\left( t_0 \right) , \varvec{u}^{*} \right) dx \le \frac{\lambda _{1,_{{\mathcal {T}}_h}}}{2} \Vert \varvec{u}\left( t_0 \right) - \varvec{u}^{*} \Vert _{L_2, {\mathcal {T}}_h}^2, \end{aligned}$$

(B.10)

where $\lambda _{1,_{{\mathcal {T}}_h}}$ is the maximum eigenvalue of ${\widetilde{A}}_0^{-1} \left( \widehat{\widehat{\varvec{u}}} \left( \theta \right) \right)$ over all elements in the domain.

Note: throughout this lemma we implicitly assume that $\varvec{u}\left( t_0 \right) = \varvec{u}\left( \varvec{v}^h \left( t_0 \right) \right)$.

Proof

The proof of part (i) follows from inspection, and the proof of part (ii) is given in Lemma B.1. Therefore, it remains to prove part (iii). One may begin the proof by utilizing Taylor’s Theorem in order to express $U \left( \varvec{u}\left( t_0 \right) \right)$ in terms of $U \left( \varvec{u}^{*} \right)$ as follows

$$\begin{aligned} U \left( \varvec{u}\left( t_0 \right) \right)&= U \left( \varvec{u}^{*} \right) + \left( U_{,\varvec{u}} \left( \varvec{u}^{*} \right) \right) ^T \left( \varvec{u}\left( t_0 \right) - \varvec{u}^{*} \right) \nonumber \\&\quad + \int _{0}^{1} \left( 1 - \theta \right) \left( \varvec{u}\left( t_0 \right) - \varvec{u}^{*} \right) ^T U_{,\varvec{u},\varvec{u}} \left( \widehat{\widehat{\varvec{u}}} \left( \theta \right) \right) \left( \varvec{u}\left( t_0 \right) - \varvec{u}^{*} \right) d \theta . \end{aligned}$$

(B.11)

Upon substituting the expression for $H \left( \varvec{u}\left( t_0 \right) , \varvec{u}^{*} \right)$ (as given by definition 7.3) into Eq. (B.11), one obtains

$$\begin{aligned} H \left( \varvec{u}\left( t_0 \right) , \varvec{u}^{*} \right)&= \int _{0}^{1} \left( 1 - \theta \right) \left( \varvec{u}\left( t_0 \right) - \varvec{u}^{*} \right) ^T U_{,\varvec{u},\varvec{u}} \left( \widehat{\widehat{\varvec{u}}} \left( \theta \right) \right) \left( \varvec{u}\left( t_0 \right) - \varvec{u}^{*} \right) d \theta \nonumber \\&= \int _{0}^{1} \left( 1 - \theta \right) \left( \varvec{u}\left( t_0 \right) - \varvec{u}^{*} \right) ^T {\widetilde{A}}_0^{-1} \left( \widehat{\widehat{\varvec{u}}} \left( \theta \right) \right) \left( \varvec{u}\left( t_0 \right) - \varvec{u}^{*} \right) d \theta , \end{aligned}$$

(B.12)

where the definition of the inverse Jacobian matrix ${\widetilde{A}}_0^{-1} = U_{,\varvec{u}, \varvec{u}}$ has been used in the last line (see Eqs. (3.4) and (3.7)). The matrix ${\widetilde{A}}_0^{-1} \left( \widehat{\widehat{\varvec{u}}} \left( \theta \right) \right)$ has a maximum eigenvalue denoted by $\lambda _1 \left( \theta \right) = \lambda _1 \left( {\widetilde{A}}_0^{-1} \left( \widehat{\widehat{\varvec{u}}} \left( \theta \right) \right) \right)$. Therefore, the last three terms that appear under the integral sign in Eq. (B.12) can be bounded from above as follows

$$\begin{aligned}&\left( \varvec{u}\left( t_0 \right) - \varvec{u}^{*} \right) ^T {\widetilde{A}}_0^{-1} \left( \widehat{\widehat{\varvec{u}}} \left( \theta \right) \right) \left( \varvec{u}\left( t_0 \right) - \varvec{u}^{*} \right) \nonumber \\&\quad \le \lambda _1 \left( \theta \right) \left( \varvec{u}\left( t_0 \right) - \varvec{u}^{*} \right) ^T \left( \varvec{u}\left( t_0 \right) - \varvec{u}^{*} \right) , \nonumber \\&\quad \le \lambda _1 \left( \theta \right) \left\| \varvec{u}\left( t_0 \right) - \varvec{u}^{*} \right\| _{L_2}^2, \end{aligned}$$

(B.13)

where several standard results from Linear Algebra (cf. for instance [61]) have been used. On substituting Eq. (B.13) into Eq. (B.12) and integrating both sides over the entire domain, one obtains

$$\begin{aligned} \sum _{T_k \in {\mathcal {T}}_h} \int _{T_k} H \left( \varvec{u}\left( t_0 \right) , \varvec{u}^{*} \right) dx&\le \sum _{T_k \in {\mathcal {T}}_h} \int _{T_k} \int _{0}^{1} \left( \left( 1 - \theta \right) \lambda _1 \left( \theta \right) \left\| \varvec{u}\left( t_0 \right) - \varvec{u}^{*} \right\| _{L_2}^2 \right) d \theta \, dx, \nonumber \\&\le \sum _{T_k \in {\mathcal {T}}_h} \left( \int _{T_k} \lambda _1 \left\| \varvec{u}\left( t_0 \right) - \varvec{u}^{*} \right\| _{L_2}^2 dx \right) \int _{0}^{1} \left( 1 - \theta \right) d \theta , \nonumber \\&\le \frac{1}{2} \lambda _{1,_{{\mathcal {T}}_h}} \left\| \varvec{u}\left( t_0 \right) - \varvec{u}^{*} \right\| _{L_2, {\mathcal {T}}_h}^2, \end{aligned}$$

(B.14)

where we have defined $\lambda _1$ to be the maximum of $\lambda _1 \left( \theta \right)$ over all $0 \le \theta \le 1$, and $\lambda _{1,_{{\mathcal {T}}_h}}$ to be the maximum of $\lambda _1$ over all elements in the domain. Finally, the proof is completed by comparing Eq. (B.14) with Eq. (B.10). $\square$

Lemma B.4

Suppose that the initial condition $\varvec{v}^h \left( t_0 \right) \in \varvec{L}_2 \left( \varOmega \right)$, the initial condition $\varvec{u}\left( t_0 \right) \in \varvec{L}_1 \left( \varOmega \right)$, each component of $\varvec{v}\left( \varvec{u}^{*} \right)$ is bounded, and the pressure p and density $\rho$ are positive and bounded for all convex combinations of states $\varvec{v}^{h} \left( t_0 \right)$ and $\varvec{v}\left( \varvec{u}^{*} \right)$ defined by $\widehat{\widehat{\varvec{v}}} \left( \theta \right) = \varvec{v}^{h} \left( t_0 \right) + \theta \left( \varvec{v}\left( \varvec{u}^{*} \right) - \varvec{v}^{h} \left( t_0 \right) \right)$, where $0 \le \theta \le 1$. Under these circumstances, the following statements hold: (i) the $L_2$ norm of $\varvec{v}^h \left( t_0 \right) - \varvec{v}\left( \varvec{u}^{*} \right)$ is well-defined; (ii) the eigenvalues of ${\widetilde{A}}_0 \left( \widehat{\widehat{\varvec{v}}} \left( \theta \right) \right)$ are real and positive; and (iii) the broken integral of ${\mathcal {H}} \left( \varvec{v}\left( \varvec{u}^* \right) , \varvec{v}^h \left( t_0 \right) \right)$ over the domain ${\mathcal {T}}_h$ is bounded in the following fashion

$$\begin{aligned} \sum _{T_k \in {\mathcal {T}}_h} \int _{T_k} {\mathcal {H}} \left( \varvec{v}\left( \varvec{u}^* \right) , \varvec{v}^h \left( t_0 \right) \right) dx \le \frac{\lambda _{1,_{{\mathcal {T}}_h}}}{2} \Vert \varvec{v}^h \left( t_0 \right) - \varvec{v}\left( \varvec{u}^{*} \right) \Vert _{L_2, {\mathcal {T}}_h}^2, \end{aligned}$$

(B.15)

where $\lambda _{1,_{{\mathcal {T}}_h}}$ is the maximum eigenvalue of ${\widetilde{A}}_0 \left( \widehat{\widehat{\varvec{v}}} \left( \theta \right) \right)$, over all elements in the domain.

Note: throughout this lemma we implicitly assume that $\varvec{u}\left( t_0 \right) = \varvec{u}\left( \varvec{v}^h \left( t_0 \right) \right)$.

Proof

The proof of this lemma is very similar to the proof of Lemma B.3. In particular, only part (iii) requires proof. One may begin the proof by using Taylor’s Theorem in order to express ${\mathcal {U}} \left( \varvec{v}\left( \varvec{u}^* \right) \right)$ in terms of ${\mathcal {U}} \left( \varvec{v}^h \left( t_0 \right) \right)$ as follows

$$\begin{aligned} {\mathcal {U}} \left( \varvec{v}\left( \varvec{u}^* \right) \right)&= {\mathcal {U}} \left( \varvec{v}^h \left( t_0 \right) \right) + \left( {\mathcal {U}}_{,\varvec{v}} \left( \varvec{v}^h \left( t_0 \right) \right) \right) ^T \left( \varvec{v}\left( \varvec{u}^* \right) - \varvec{v}^h \left( t_0 \right) \right) \nonumber \\&\quad + \int _{0}^{1} \left( 1 - \theta \right) \left( \varvec{v}\left( \varvec{u}^* \right) - \varvec{v}^h \left( t_0 \right) \right) ^T {\mathcal {U}}_{,\varvec{v},\varvec{v}} \left( \widehat{\widehat{\varvec{v}}} \left( \theta \right) \right) \left( \varvec{v}\left( \varvec{u}^* \right) - \varvec{v}^h \left( t_0 \right) \right) d \theta . \end{aligned}$$

(B.16)

On substituting the expression for ${\mathcal {H}} \left( \varvec{v}\left( \varvec{u}^* \right) , \varvec{v}^h \left( t_0 \right) \right)$ (as given by Definition 7.6) into Eq. (B.16), one obtains

$$\begin{aligned}&{\mathcal {H}} \left( \varvec{v}\left( \varvec{u}^* \right) , \varvec{v}^h \left( t_0 \right) \right) \nonumber \\&\quad = \int _{0}^{1} \left( 1 - \theta \right) \left( \varvec{v}\left( \varvec{u}^* \right) - \varvec{v}^h \left( t_0 \right) \right) ^T {\mathcal {U}}_{,\varvec{v},\varvec{v}} \left( \widehat{\widehat{\varvec{v}}} \left( \theta \right) \right) \left( \varvec{v}\left( \varvec{u}^* \right) - \varvec{v}^h \left( t_0 \right) \right) d \theta \nonumber \\&\quad = \int _{0}^{1} \left( 1 - \theta \right) \left( \varvec{v}\left( \varvec{u}^* \right) - \varvec{v}^h \left( t_0 \right) \right) ^T {\widetilde{A}}_0 \left( \widehat{\widehat{\varvec{v}}} \left( \theta \right) \right) \left( \varvec{v}\left( \varvec{u}^* \right) - \varvec{v}^h \left( t_0 \right) \right) d \theta , \end{aligned}$$

(B.17)

where the definition of the Jacobian matrix ${\widetilde{A}}_0 = {\mathcal {U}}_{,\varvec{v}, \varvec{v}}$ has been used in the last line [see Eqs. (3.2) and (3.7)]. The remainder of the proof requires a simple eigenvalue analysis in order to find the maximum eigenvalue $\lambda _1 \left( \theta \right) = \lambda _1 \left( {\widetilde{A}}_0 \left( \widehat{\widehat{\varvec{v}}} \left( \theta \right) \right) \right)$, and thereafter, the construction of an upper bound with the maximum eigenvalue and the $L_2$ norm

$$\begin{aligned} \Vert \varvec{v}^h \left( t_0 \right) - \varvec{v}\left( \varvec{u}^{*} \right) \Vert _{L_2}, \end{aligned}$$

in accordance with the proof of Lemma B.3. $\square$

Lemma B.5

The functional $\; U \left( \varvec{u}\right)$ is strongly convex on the set ${{\mathscr {C}} \subset \varvec{\mathrm {dom}} \; U \left( \varvec{u}\right) }$, under the assumptions that ${\mathscr {C}}$ is a convex set, and that the eigenvalues of the matrix ${\widetilde{A}}_{0}^{-1} \left( \varvec{u}\right)$ are bounded away from zero for all $\varvec{u}\in {\mathscr {C}}$, where ${\widetilde{A}}_{0}^{-1}$ is defined in Eq. (3.7).

Proof

In order to prove the strong convexity of U, one must obtain a particular inequality governing the Hessian, $U_{, \varvec{u}, \varvec{u}}$, for all $\varvec{u}\in {\mathscr {C}}$. Towards this end, one may first note that the following identity holds in accordance with Eqs. (3.4) and (3.7)

$$\begin{aligned} U_{, \varvec{u}, \varvec{u}} = \varvec{v}_{,\varvec{u}} = {\widetilde{A}}_{0}^{-1}. \end{aligned}$$

Here, the inverse Jacobian (${\widetilde{A}}_{0}^{-1}$) is SPD because the Jacobian itself (${\widetilde{A}}_{0}$) is SPD under the assumptions of Lemma B.1. As a result, the eigenvalues of $U_{, \varvec{u}, \varvec{u}} = {\widetilde{A}}_{0}^{-1}$ are greater than zero, and are guaranteed to satisfy the following inequality

$$\begin{aligned} \lambda _{1} \ge \cdots \ge \lambda _{m} > 0. \end{aligned}$$

Furthermore, it is common practice (cf. [6, 62]) to assume that the minimum eigenvalue of ${\widetilde{A}}_{0}^{-1}$, for all $\varvec{u}\in {\mathscr {C}}$, is bounded away from zero

$$\begin{aligned} \lambda _{m} \ge 2 C, \end{aligned}$$

where $C > 0$ is a constant independent of $\varvec{u}$. Setting this result aside for the moment, one may utilize the Hessian $U_{, \varvec{u}, \varvec{u}}$ in order to construct a quadratic form

$$\begin{aligned} \varvec{u}_{b}^T \left( U_{, \varvec{u}, \varvec{u}} \left( \varvec{u}_{a} \right) \right) \varvec{u}_{b}, \end{aligned}$$

(B.18)

where $\varvec{u}_{a}$ and $\varvec{u}_{b}$ are arbitrary values of $\varvec{u}\in {\mathscr {C}}$. In accordance with several standard results from Linear Algebra (cf. for instance [61]), the quadratic form in Eq. (B.18) satisfies the following inequalities

$$\begin{aligned} \varvec{u}_{b}^T \left( U_{, \varvec{u}, \varvec{u}} \left( \varvec{u}_{a} \right) \right) \varvec{u}_{b}&\ge \lambda _m \varvec{u}_{b}^T \varvec{u}_{b}, \nonumber \\&= \lambda _m \left\| \varvec{u}_{b} \right\| _{L_2}^2, \nonumber \\&\ge 2 C \left\| \varvec{u}_{b} \right\| _{L_2}^2. \end{aligned}$$

(B.19)

From Eq. (B.19) and the strong convexity criterion in [7], p. 459, it immediately follows that $U \left( \varvec{u}\right)$ is strongly convex. $\square$

Lemma B.6

The functional $\; {\mathcal {U}} \left( \varvec{v}\right)$ is strongly convex on the set ${{\mathscr {C}} \subset \varvec{\mathrm {dom}} \; {\mathcal {U}} \left( \varvec{v}\right) }$, under the assumptions that ${\mathscr {C}}$ is a convex set, and that the eigenvalues of the matrix ${\widetilde{A}}_{0} \left( \varvec{v}\right)$ are bounded away from zero for all $\varvec{v}\in {\mathscr {C}}$, where ${\widetilde{A}}_0$ is defined in Eq. (3.7).

Proof

The proof of this lemma is very similar to the proof of Lemma B.5. The proof begins with the observation that ${\widetilde{A}}_{0}$ is SPD under the assumptions of Lemma B.1, and that

$$\begin{aligned} {\mathcal {U}}_{, \varvec{v}, \varvec{v}} = \varvec{u}_{,\varvec{v}} = {\widetilde{A}}_{0}, \end{aligned}$$

in accordance with Eqs. (3.2) and (3.7). The rest of the proof requires a simple eigenvalue analysis, and the construction of inequalities involving the quadratic form,

$$\begin{aligned} \varvec{v}_{b}^T \left( {\mathcal {U}}_{, \varvec{v}, \varvec{v}} \left( \varvec{v}_{a} \right) \right) \varvec{v}_{b}, \end{aligned}$$

in accordance with the proof of Lemma B.5. $\square$

Appendix C: Algebraically stable SDIRK methods

In this section, the Butcher tables for some well-known algebraically stable SDIRK methods are presented.

The 1-stage 1st-order ‘Backward Euler’ method has the following Butcher table

$$\begin{aligned} \begin{array}{c|c} 1 &{} 1 \\ \hline &{} 1 \end{array} \end{aligned}$$

The 2-stage 2nd-order method due to [1] and [8] has the following Butcher table

$$\begin{aligned} \begin{array}{c|ccc} \zeta &{} \zeta &{} 0\\ c_2 &{} a_{21} &{} \zeta \\ \hline &{}b_1 &{} b_2 \end{array} \end{aligned}$$

where

$$\begin{aligned} \zeta&= \frac{\left( 2 \pm \sqrt{2} \right) }{2}, \\ a_{21}&= 1 - 2 \zeta , \qquad b_1 = \frac{1}{2}, \qquad b_2 = \frac{1}{2}, \qquad c_2 = 1 - \zeta . \end{aligned}$$

The 3-stage 3rd-order method due to [1] and [8] has the following Butcher table

$$\begin{aligned} \begin{array}{c|ccc} \zeta &{} \zeta &{} 0 &{} 0\\ c_2 &{} a_{21} &{} \zeta &{} 0 \\ c_3 &{} a_{31} &{} a_{32} &{} \zeta \\ \hline &{}b_1 &{} b_2 &{} b_3 \end{array} \end{aligned}$$

where

$$\begin{aligned} \zeta&= 0.43586652150845899941601945, \\ | d_1 |&\ge 1.774294247072785073096332, \qquad d_2 = \left( \frac{1}{2} - \zeta - d_1 \right) \frac{b_2}{b_3}, \\ a_{21}&= \frac{1}{2} - \zeta +d_1, \qquad a_{31} = 1 - 2 \zeta - d_2, \qquad a_{32} = d_2, \\ b_1&= \frac{d_1 \left( 1 -2 \zeta \right) + \frac{1}{6}}{\left( 2 \zeta - 1\right) \left( 2 \zeta - 1 -2 d_1 \right) }, \qquad b_2 = \frac{\zeta ^2 - \zeta + \frac{1}{6}}{\left( \zeta - \frac{1}{2} \right) ^2 - d_1^2}, \\&\qquad b_3 = \frac{d_1 \left( 2 \zeta - 1 \right) + \frac{1}{6}}{\left( 2 \zeta -1 \right) \left( 2 \zeta -1 + 2 d_1 \right) }, \\ c_2&= \frac{1}{2} + d_1, \qquad c_3 = 1 - \zeta . \end{aligned}$$

Note that the expression for $d_{2}$ given here is different than the corresponding expression in [8]. The definition given here is the correct definition, as the original definition contains a typographical error that prevents the method from being algebraically stable.

Finally, the 4-stage 4th-order method due to [1, 8] has the following Butcher table

$$\begin{aligned} \begin{array}{c|cccc} \zeta &{} \zeta &{} 0 &{} 0 &{} 0\\ c_2 &{} a_{21} &{} \zeta &{} 0 &{} 0\\ c_3 &{} a_{31} &{} a_{32} &{} \zeta &{} 0 \\ c_4 &{} a_{41} &{} a_{42} &{} a_{43} &{} \zeta \\ \hline &{}b_1 &{} b_2 &{} b_3 &{} b_4 \end{array} \end{aligned}$$

where

$$\begin{aligned} \zeta&= 0.5728160624821348554080014, \\ d_1&= \zeta ^2 - \frac{\zeta }{2} + \frac{1}{12}, \qquad d_2 = -\frac{\zeta \left( \zeta -\frac{1}{3} \right) }{\left( 2 \zeta ^2 - \zeta + \frac{1}{6} \right) }, \qquad d_3 =\zeta ^3 -\frac{3 \zeta ^2}{2} + \frac{\zeta }{2} - \frac{1}{24}, \\ a_{21}&= \frac{\zeta \left( \frac{1}{3} - \zeta \right) }{2 d_1}, \qquad a_{31} = \frac{1}{2} - \zeta - \frac{d_3}{d_1} + \frac{8 d_3 \left( \zeta - \frac{1}{4} \right) }{\zeta \left( \zeta - \frac{1}{3} \right) } , \qquad a_{32} = -\frac{8 d_3 \left( \zeta - \frac{1}{4} \right) }{\zeta \left( \zeta - \frac{1}{3} \right) }, \\ a_{41}&= 1 - 2 \zeta - a_{42} - \frac{d_1}{\left( 4 \zeta -1 \right) \left( b_2 - \frac{1}{2} \right) }, \qquad a_{42} = \frac{ \left( \zeta ^2 - \zeta + \frac{1}{6} \right) /2 + d_3/d_2}{d_2 \left( \frac{1}{2} - b_2 \right) }, \\ a_{43}&= \frac{d_1}{\left( 4 \zeta -1 \right) \left( b_2 - \frac{1}{2} \right) }, \qquad b_1 = \frac{1}{2} - b_2, \qquad b_2 = \frac{d_1^2}{2 \zeta \left( \zeta - \frac{1}{3} \right) \left( \zeta - \frac{1}{4} \right) }, \qquad b_3 = b_2, \\ b_4&= \frac{1}{2} - b_2, \qquad c_2 = \frac{ \zeta \left( \zeta - \frac{1}{2} \right) ^2}{d_1}, \qquad c_3 = \frac{1}{2} - \frac{d_3}{d_1}, \qquad c_4 = 1- \zeta . \end{aligned}$$

Note that the expression for $a_{31}$ given here is different than the corresponding expression in [8]. The definition given here is the correct definition, as the original definition contains a typographical error that prevents the method from being algebraically stable.