# Algebraic Conditions for Stability in Runge-Kutta Methods and Their
Certification via Semidefinite Programming

Austin Juhl, David Shirokoff Department of Mathematical Sciences, NJIT,
Newark, NJ 07102,
<EMAIL_ADDRESS>of Mathematical Sciences, NJIT, Newark, NJ 07102,
<EMAIL_ADDRESS>
(May 2024)

###### Abstract

In this work, we present approaches to rigorously certify $A$\- and
$A(\alpha)$-stability in Runge-Kutta methods through the solution of convex
feasibility problems defined by linear matrix inequalities. We adopt two
approaches. The first is based on sum-of-squares programming applied to the
Runge-Kutta $E$-polynomial and is applicable to both $A$\- and
$A(\alpha)$-stability. In the second, we sharpen the algebraic conditions for
$A$-stability of Cooper, Scherer, Türke, and Wendler to incorporate the Runge-
Kutta order conditions. We demonstrate how the theoretical improvement enables
the practical use of these conditions for certification of $A$-stability
within a computational framework. We then use both approaches to obtain
rigorous certificates of stability for several diagonally implicit schemes
devised in the literature.

Keywords: Runge-Kutta, $A$-stability, $A(\alpha)$-stability, algebraic
characterization, semidefinite programming.

AMS Subject Classifications: 65L06, 65L07, 65L20.

##  1 Introduction

In recent years, sum-of-squares optimization and semidefinite programming have
become valuable tools in developing rigorous certificates of stability in
dynamical systems. Such certificates are useful in developing reliable
algorithms and software. Examples range from stochastic linear-quadratic
control [29], switched stability of nonlinear systems [1], reinforcement
learning [20], stability in partial differential equations [13], and
applications for robotic control [22]. Additionally, control systems have
benefited from analogous certification techniques, including verifying the
stability of reinforcement learning policies [12] or variable step-size
convergence bounds for gradient descent [15]. In this work, we adopt a similar
strategy to establish rigorous certificates for $A$\- and $A(\alpha)$\-
stability in Runge-Kutta methods through the computational solution of linear
matrix inequalities via semi-definite programming.

###  1.1 Background

Numerical stability is a critical property for a time-integration scheme. In
the context of Runge-Kutta (RK) methods applied to stiff differential
equations, $A$-stability (or the related $A(\alpha)$-stability), introduced by
Dahlquist [9], is one of the most basic and practically important notions of
stability. As a result of the reduced degrees of freedom in the Butcher
coefficients, Runge-Kutta methods that are $A$-stable with stability functions
that are near-optimal order approximations have been well characterized [6,
28, 16, 23] (cf. [18] for algebraically stable methods). However, many
existing and newly developed Runge-Kutta methods in the literature admit a
diagonally implicit structure and lie outside the application of these
previous results. While one strategy for verifying $A$-stability is via _Sturm
sequences_ [19], our approach here is rooted in convex optimization.

Dating back to the work of Dahlquist (e.g., [10]), it has been known that
$A$-stability is equivalent to the RK stability function satisfying a
particular convex feasibility problem — i.e., up to a transformation, the
stability function lies in the convex cone of _positive functions_ [19,
Chapter IV.5]. Subsequent convex feasibility conditions for $A$-stability also
include: (1) The RK $E$-polynomial lying in the convex cone of non-negative
polynomials; and (2) Satisfying a set of algebraic conditions, which we refer
to as the CSTW conditions, developed by Cooper [7], Scherer and Türke [25],
and Scherer and Wendler [26]. The CSTW conditions form a feasibility problem
over the convex cone of semidefinite matrices. It is worth noting that
semidefinite matrix conditions exist for stronger notions of stability, such
as _algebraic stability_ (or the related concept of $B$-stability) for Runge-
Kutta methods [19, Chapter IV.12], and $G$-stability [11] for linear multistep
methods.

###  1.2 Contributions and Organization of the Paper

The focus of this work lies in the practical aspects of implementing and
generating rigorous certificates of stability through solutions of convex
feasibility problems for $A$\- and $A(\alpha)$-stability. The main
contributions are twofold. First, we provide a theoretical contribution that
sharpens the CSTW conditions. The CSTW conditions do not take into account the
fact that the stability function is a $p$th order approximation to the
exponential — which ultimately limits their practical use in providing
rigorous certificates via computer-assisted means. Our main theoretical result
(Theorem 4.7) modifies the CSTW conditions to account for the RK order
conditions, which then enables the rigorous certification of stability via
computational approaches. In addition to sharpening the CSTW conditions, we
provide details from sum-of-squares optimization for certifying $A$\- and
$A(\alpha)$-stability via testing the non-negativity of the $E$-polynomial.

Our second contribution provides rigorous certificates of stability for
several recently devised schemes in the literature. These schemes include the
7-stage 4th order and 12-stage 5th order schemes from [5] as well as the
9-stage 6th order and 13-stage 8th order schemes from [2].

The paper is organized as follows: §2 and §3 cover background material on
Runge-Kutta methods and linear matrix inequalities, respectively. Linear
matrix inequalities for $A$\- and $A(\alpha)$\- stability of Runge-Kutta
methods are then developed in §4. This section also includes our main
theoretical result along with illustrative examples. We then turn our
attention in §5 to constructing rigorous certificates of $A$\- and
$A(\alpha)$-stability via a computer-assisted solution to the associated
linear matrix inequalities. A summary and outlook are provided in §6. The
Appendix and companion Github repository provide certificates of stability for
the schemes we examine as well as the supporting numerical code.

##  2 Runge-Kutta Background

Numerical time integration of an ordinary differential equation
($\mathbb{V}=\mathbb{R}^{m}$ or $\mathbb{C}^{m}$)

$\boldsymbol{u}^{\prime}(t)=f\big{(}t,\boldsymbol{u}(t)\big{)}\,,\quad\boldsymbol{u}(0)=\boldsymbol{u}_{0}\,;\quad\boldsymbol{u}\in\mathbb{V}\,,\quad
f:\mathbb{R}\times\mathbb{V}\to\mathbb{V}\,,$ (1)

results in a discrete-in-time dynamical system — one of the most common being
Runge-Kutta (RK) methods. Runge-Kutta methods, discretize (1) with $s$ stages
as

$\displaystyle\boldsymbol{h}_{i}$ $\displaystyle=\boldsymbol{u}_{n}+\Delta
t\,\sum_{j=1}^{s}a_{ij}\,f(t_{n}+c_{j}\Delta t,\boldsymbol{h}_{j})\;,\quad
i=1,2,\ldots,s\;$ $\displaystyle\boldsymbol{u}_{n+1}$
$\displaystyle=\boldsymbol{u}_{n}+\Delta
t\,\sum_{j=1}^{s}b_{j}\,f(t_{n}+c_{j}\Delta t,\boldsymbol{h}_{j})\;\,,$

where $\boldsymbol{u}_{n}\approx\boldsymbol{u}(t_{n})$ and $t_{n}$ the $n$th
time step. The RK scheme is defined by coefficients

$\mbox{\boldmath$A$}=[a_{ij}]_{i,j=1}^{s}\,,\qquad\boldsymbol{b}=[b_{1},\ldots,b_{s}]^{T}\,,\quad\textrm{and}\quad\boldsymbol{c}=[c_{1},\ldots,c_{s}]^{T}:=\mbox{\boldmath$A$}\boldsymbol{e}\,;\quad\textrm{where}\quad\boldsymbol{e}=[1,\ldots,1]^{T}\,.$

###  2.1 Linear Stability

For the scalar linear case where $f(\boldsymbol{u})=\lambda\boldsymbol{u}$
with $\lambda\in\mathbb{C}$ and fixed $\Delta t$, the Runge-Kutta dynamics (2)
are expressed as:

$\displaystyle\boldsymbol{u}_{n+1}=W(z)\boldsymbol{u}_{n},\quad\text{where}\quad
z:=\lambda\Delta t\,.$

Here $W(z)$ is the _stability function_ of the scheme given by:

$W(z):=1+z\boldsymbol{b}^{T}(\mbox{\boldmath$I$}-z\mbox{\boldmath$A$})^{-1}\boldsymbol{e}=\frac{\det(\mbox{\boldmath$I$}-z\mbox{\boldmath$A$}+z\boldsymbol{e}\boldsymbol{b}^{T})}{\det(\mbox{\boldmath$I$}-z\mbox{\boldmath$A$})}=\frac{N(z)}{D(z)}.$
(3)

In equation (3), $N(z)$ and $D(z)$ represent polynomials of degree at most
$s$, sharing no common factors, and $I$ is the identity matrix.

A method exhibits a _degenerate_ stability function if $\deg N\leq s-1$ and
$\deg D\leq s-1$; otherwise, the stability function is considered non-
degenerate. Degenerate stability functions occur when
$\det(\mbox{\boldmath$I$}-z\mbox{\boldmath$A$}+z\boldsymbol{e}\boldsymbol{b}^{T})$
and $\det(\mbox{\boldmath$I$}-z\mbox{\boldmath$A$})$ share a common root. It
is worth noting that degenerate stability functions can appear in practice,
for instance, in constructing RK schemes that avoid order reduction (e.g.,
[5]).

The dynamics (2) are _stable_ for a given $z$ if $|W(z)|\leq 1$. We will
primarily be concerned with numerical schemes
$(\mbox{\boldmath$A$},\boldsymbol{b})$ that are $A$-stable, i.e., those where

$A$-stability: $\displaystyle\qquad|W(z)|\leq 1\qquad\textrm{for}\qquad
z\in\mathbb{C}_{-}=\\{z\in\mathbb{C}\;:\;\textrm{Re}\,z\leq 0\\}\,.$

The $A$-stability criteria ensures that the discrete dynamics (2) are stable
whenever the linear ODE (1) is stable. For schemes that are not $A$-stable, it
is useful to characterize the largest $\alpha>0$ for which the RK dynamics are
stable for all $z$ in the sector $S_{\alpha}$ opening into the left-half plane
with angle $\alpha$

$\displaystyle A(\alpha)\textrm{-stability:}$ $\displaystyle|W(z)|\leq
1\qquad\textrm{for}\qquad$ $\displaystyle z\in
S_{\alpha}=\\{z\in\mathbb{C}\,:\,|\arg(-z)|\leq\alpha\\}\,.$

Note that $A$-stability is equivalent to $A(\alpha)$-stability with
$\alpha=\frac{\pi}{2}$. Both $A$\- and $A(\alpha)$\- stability can be recast
in terms of the non-negativity of the $E$-polynomial [19].

In a slight generalization from the standard definition, to allow for
$\alpha\neq\frac{\pi}{2}$, the (generalized) $E$-polynomial is

$\displaystyle
E(y;\alpha)=|D(ye^{-i\alpha})|^{2}-|N(ye^{-i\alpha})|^{2}\,\quad\textrm{where}\quad
0<\alpha\leq\frac{\pi}{2}\,,\;y\in\mathbb{R}\,.$ (4)

A scheme is then $A(\alpha)$-stable if (for $\alpha=\frac{\pi}{2}$ see [19,
Chapter IV.3]) $W(z)$ is analytic in the interior of $S_{\alpha}$ and

$\displaystyle E(y;\alpha)\geq 0\qquad\textrm{for all}\qquad y\geq 0\,.$ (5)

The condition (5) guarantees $A(\alpha)$-stability by ensuring $|W(z)|\leq 1$
for all $z\in\partial S_{\alpha}$. As $W(z)$ is analytic within the interior
of $S_{\alpha}$, the maximum modulus principle implies that $|W(z)|$ is
maximized on $\partial S_{\alpha}$, confirming $A(\alpha)$-stability.

###  2.2 Accuracy and Order Conditions

For an RK scheme to achieve _(classical) order $p$_ on linear, autonomous
problems, the stability function must approximate the exponential function to
order $p$, such that $W(z)=e^{z}+\mathcal{O}(z^{p+1})$ as $z\rightarrow 0$.
This approximation is achieved provided the RK coefficients
$(\mbox{\boldmath$A$},\boldsymbol{b})$ satisfy the _tall-tree_ order
conditions of order $p$:

$\boldsymbol{b}^{T}\mbox{\boldmath$A$}^{j-1}\,\boldsymbol{e}=\frac{1}{j!}\,\quad\text{for}\quad
1\leq j\leq p\,,$ (6)

as outlined in [17]. Additional RK order conditions, i.e., the non-tall-tree
conditions, are further required to achieve accuracy of order $p$ on general
(nonlinear) ODEs.

If $\alpha=\frac{\pi}{2}$, the $p$th order conditions (6) imply that
$E(y;\frac{\pi}{2})$ admits a factor of $y^{2j}$ [19, Chapter IV.3]:

$\displaystyle
E\left(y;\frac{\pi}{2}\right)=\mathcal{O}(y^{2j})\,,\qquad\textrm{as}\quad
y\rightarrow 0\quad\textrm{where},\quad
j\geq\left\lfloor\frac{p}{2}\right\rfloor+1\,.$ (7)

In general for $\alpha<\frac{\pi}{2}$, the polynomial (as implied by the
calculations in [7])

$\displaystyle E\left(y;\alpha\right)=\mathcal{O}(y)\,,\qquad\textrm{as}\quad
y\rightarrow 0\,.$ (8)

##  3 Linear Matrix Inequalities Background

Certifying numerical stability will reduce to a semidefinite feasibility
program involving linear matrix inequalities (LMIs). This section reviews LMIs
and their relationships to convex feasibility as well as non-negative
polynomials.

###  3.1 Feasibility and Convexity of Linear Matrix Inequalities

Given
$\mbox{\boldmath$P$},\mbox{\boldmath$N$}_{1},\,\ldots,\,\mbox{\boldmath$N$}_{n}\in\mathbb{S}^{n}$,
the set of $n\times n$ real symmetric matrices, a _linear matrix inequality_
(LMI) is defined as:

$\displaystyle\mbox{\boldmath$F$}(\boldsymbol{\eta}):=\mbox{\boldmath$P$}+\sum_{j=1}^{d}\eta_{j}\mbox{\boldmath$N$}_{j}\succeq
0\,,$ (9)

where $\mbox{\boldmath$F$}\succeq 0$ indicates that $F$ is positive semi-
definite ($\mbox{\boldmath$F$}\succ 0$ is positive definite). The LMI (9) is
_feasible_ if there exists a vector $\boldsymbol{\eta}$ such that
$\mbox{\boldmath$F$}(\boldsymbol{\eta})\succeq 0$; otherwise the LMI is
_infeasible_. The linearity of $\mbox{\boldmath$F$}(\boldsymbol{\eta})$
ensures that the set:

$\displaystyle\mathcal{C}=\\{\boldsymbol{\eta}\in\mathbb{R}^{d}\;:\;\mbox{\boldmath$F$}(\boldsymbol{\eta})\succeq
0\\}\,,$

containing all $\boldsymbol{\eta}$ that satisfy the LMI (9), is convex. Thus,
assessing the feasibility of $\mbox{\boldmath$F$}(\boldsymbol{\eta})$ —
determining whether $\mathcal{C}$ is non-empty — is a convex feasibility
problem and can be solved via semidefinite programming.

Due to the matrix structure of (9), the feasible set $\mathcal{C}$ may lie in
an affine plane having dimension less than $d$, potentially resulting in
$\mathcal{C}$ having an empty interior. To quantify this feature, the
dimension for the set $\mathcal{C}$ may be characterized in terms of the
affine hull:

$\displaystyle\textrm{aff}(\mathcal{C})=\\{\mu_{1}\boldsymbol{\eta}_{1}+\mu_{2}\boldsymbol{\eta}_{2}\;:\;\mu_{1},\mu_{2}\in\mathbb{R},\;\textrm{and
}\boldsymbol{\eta}_{1},\boldsymbol{\eta}_{2}\in\mathcal{C}\\}=\boldsymbol{\eta}_{0}+\mathbb{V}\,,$

where $\boldsymbol{\eta}_{0}\in\mathcal{C}$ and $\mathbb{V}$ is a vector
space. The dimension of $\mathcal{C}$ is defined as $\dim\mathbb{V}.$

###  3.2 Non-negative Polynomials as Linear Matrix Inequalities

Here, we review characterizing the non-negativity of a polynomial through
linear matrix inequalities.

Let $\mathbb{R}[x]$ denote the set of single-variable polynomials with real
coefficients. Two convex cones within $\mathbb{R}[x]$ include the set of non-
negative polynomials, satisfying $p(x)\geq 0$ for all $x\in\mathbb{R}$, and
the set of sum-of-squares (SOS) polynomials, where each polynomial can be
expressed as $\sum_{j=1}^{\ell}q_{j}^{2}(x)$ for some
$\\{q_{j}\\}_{j=1}^{\ell}\in\mathbb{R}[x]$. In one dimension, these two cones
coincide, meaning that a polynomial $p(x)\in\mathbb{R}[x]$ is non-negative for
all real numbers if and only if it can be decomposed into a sum-of-squares
[21, Theorem 2.5]. For general multivariate polynomials, every SOS polynomial
is non-negative; however, the reverse is not true.

Determining whether a polynomial is SOS and non-negative can be formulated as
an LMI. Consider the subspace of symmetric matrices $\mathcal{N}_{m}$ defined
by:

$\displaystyle\mathcal{N}_{m}:=\big{\\{}\mbox{\boldmath$N$}\in\mathbb{S}^{m}\;:\;\boldsymbol{y}^{T}\mbox{\boldmath$N$}\boldsymbol{y}=0\big{\\}}\,,\quad\textrm{where}\quad\boldsymbol{y}=[1,y,\ldots,y^{m-1}]\,.$

The components $(n_{i,j})_{i,j=1}^{m}$ of a symmetric matrix
$\mbox{\boldmath$N$}\in\mathcal{N}_{m}$ must satisfy:

$\displaystyle\sum_{i+j=r}^{m}n_{i,j}=0\qquad\textrm{for}\quad
r=2,\;3,\;\ldots,\;2m\,.$

Any polynomial $p\in\mathbb{R}[x]$ can then be expressed non-uniquely in
factorized form as:

$\displaystyle
p(x)=p_{0}+p_{1}y+\ldots+p_{2m-2}\,y^{2m-2}=\boldsymbol{y}^{T}(\mbox{\boldmath$P$}+\mbox{\boldmath$N$})\boldsymbol{y}\,,$
(10)

where $\mbox{\boldmath$N$}\in\mathcal{N}_{m}$ and

$\displaystyle\mbox{\boldmath$P$}=\begin{pmatrix}p_{0}&\tfrac{1}{2}p_{1}&&\\\
\tfrac{1}{2}p_{1}&p_{2}&\ddots\vphantom{\ddots}&\\\
&\ddots&\ddots&\tfrac{1}{2}p_{2m-3}\\\
&&\tfrac{1}{2}p_{2m-3}&p_{2m-2}\end{pmatrix}\in\mathbb{R}^{m\times m}\,.$

The polynomial $p$, described in (10), is a sum-of-squares if and only if
there exists an $\mbox{\boldmath$N$}\in\mathcal{N}_{m}$ such that
$\mbox{\boldmath$P$}+\mbox{\boldmath$N$}\succeq 0$. If
$\mbox{\boldmath$Q$}^{T}\mbox{\boldmath$Q$}=\mbox{\boldmath$P$}+\mbox{\boldmath$N$}$
is a Cholesky factorization, then defining
$q_{j}(y)=\boldsymbol{e}_{j}^{T}\mbox{\boldmath$Q$}\boldsymbol{y}$, where
$\boldsymbol{e}_{j}$ is the $j$th unit vector, admits

$\displaystyle
p(y)=\|\mbox{\boldmath$Q$}\boldsymbol{y}\|^{2}=\sum_{j=1}^{m}q_{j}^{2}(y)\,.$

Conversely, if $p$ is an SOS, the coefficients of $q_{j}$ define a positive
definite matrix $(\mbox{\boldmath$P$}+\mbox{\boldmath$N$})$.

To solve subsequent SDPs, it is useful to define a basis for
$\mathcal{N}_{m}$, which has dimension

$d:=\frac{1}{2}(m-1)(m-2).$ (11)

First, consider the index set

$\displaystyle\mathcal{S}_{m}:=\big{\\{}(i,j)\in\mathbb{N}^{2}\;|\;1\leq i\leq
m-2\,,i+2\leq j\leq m\big{\\}}\,\qquad\textrm{for}\qquad m\geq 1\,,$

which has $d$ elements and is empty for $m=1,2$. A basis for $\mathcal{N}_{m}$
is given by $\\{\mbox{\boldmath$N$}_{l}\\}_{l=1}^{d}$:

$\mbox{\boldmath$N$}_{l}=\boldsymbol{e}_{i}\boldsymbol{e}_{j}^{T}+\boldsymbol{e}_{j}\boldsymbol{e}_{i}^{T}-\boldsymbol{e}_{\left\lfloor\frac{i+j}{2}\right\rfloor}\boldsymbol{e}_{\left\lceil\frac{i+j}{2}\right\rceil}^{T}-\boldsymbol{e}_{\left\lceil\frac{i+j}{2}\right\rceil}\boldsymbol{e}_{\left\lfloor\frac{i+j}{2}\right\rfloor}^{T}\,,$

where

$\displaystyle(i,j)\in\mathcal{S}_{m}\qquad\textrm{and}\qquad
l=m(i-1)-\frac{1}{2}i(i+3)+j\,.$

For example, $m=3$ has the basis matrix:

$\mbox{\boldmath$N$}_{1}=\begin{bmatrix}0&0&1\\\ 0&-2&0\\\
1&0&0\end{bmatrix},$

while $m=4$ has three basis matrices:

$\mbox{\boldmath$N$}_{1}=\begin{bmatrix}0&0&1&0\\\ 0&-2&0&0\\\ 1&0&0&0\\\
0&0&0&0\end{bmatrix},\quad\mbox{\boldmath$N$}_{2}=\begin{bmatrix}0&0&0&1\\\
0&0&-1&0\\\ 0&-1&0&0\\\
1&0&0&0\end{bmatrix},\quad\mbox{\boldmath$N$}_{3}=\begin{bmatrix}0&0&0&0\\\
0&0&0&1\\\ 0&0&-2&0\\\ 0&1&0&0\end{bmatrix}.$

With these notations, a polynomial $p(y)$, given by (10), is non-negative if
and only if there exists $\boldsymbol{\eta}\in\mathbb{R}^{d}$ such that the
following LMI is satisfied:

$\displaystyle\mbox{\boldmath$F$}(\boldsymbol{\eta}):=\mbox{\boldmath$P$}+\sum_{l=1}^{d}\eta_{l}\,\mbox{\boldmath$N$}_{l}\succeq
0\,.$ (12)

##  4 $A$\- and $A(\alpha)$-Stability as Linear Matrix Inequalities

To certify the stability of Runge-Kutta methods, we utilize linear matrix
inequalities in two approaches. The first approach leverages the non-
negativity of the (generalized) $E$-polynomial, suitable for both $A$\- and
$A(\alpha)$-stability. The second approach uses the algebraic conditions for
$A$-stability established by Cooper [7], Scherer and Türke [25], and Scherer
and Wendler [26]. Our main theoretical contribution sharpens these algebraic
conditions, which enables their practical use within an SDP framework.

###  4.1 The $E$-polynomial LMI for $A(\alpha)$-stability

An LMI for $A(\alpha)$-stability follows immediately by applying the LMI for
general non-negative polynomials, as described in § 3.2, to $E(y;\alpha)$ in
(5).

In particular, based on the asymptotics (7)–(8), let $F(y)$ be the polynomial
from which the largest even monomial (ensured by the order conditions) has
been factored out:

$\displaystyle\textrm{If}\;\alpha$ $\displaystyle=\frac{\pi}{2}\;$
$\displaystyle:\qquad F(y)$
$\displaystyle:=y^{-2\kappa}\,E(y;\frac{\pi}{2}),\quad\kappa=\left\lfloor\frac{p}{2}\right\rfloor+1\,,$
$\displaystyle\textrm{If}\;\alpha$ $\displaystyle<\frac{\pi}{2}\;$
$\displaystyle:\qquad F(y)$ $\displaystyle:=y^{-2}\,E(y^{2};\alpha)\,.$

In both cases, $F(y)$ is an even polynomial:

$\displaystyle F(y)=p_{0}+p_{2}\,y^{2}+\ldots+p_{2m-2}\,y^{2m-2}\,,$ (13)

where the coefficients of $F$ are polynomial functions of the RK scheme
coefficients and, in the case of $\alpha<\tfrac{\pi}{2}$, polynomial functions
of $\beta:=\cos(\alpha)$.

Combining the polynomial $F(y)$ and the LMI (12), it follows:

######  Lemma 4.1.

A scheme is $A(\alpha)$-stable if:

  1. 1.

$A$ has eigenvalues outside $S_{\alpha}$ (so that $W(z)$ is analytic in
$S_{\alpha}$); and

  2. 2.

The LMI (12) is feasible for the $E$-polynomial (13).

###  4.2 The CSTW Algebraic Conditions and LMI for $A$-stability

Closely related to the $E$-polynomial approach is an LMI based on algebraic
conditions for $A$-stability, laid out in a line of work by Cooper [7],
Scherer and Türke [25] and Scherer and Wendler [26].

Cooper [7] initially established sufficient algebraic conditions for
$A$-stability by factorizing the $E$-polynomial into a quadratic form. Scherer
and Türke [25] later re-derived almost identical conditions by applying the
Kalman-Yakubovich-Popov (KYP) Lemma to $W(z)$; this further showed necessity
(in addition to sufficiency) of the algebraic conditions of Cooper. Subsequent
work by Scherer and Wendler [26] provided even more general algebraic
conditions applicable to degenerate stability functions.

######  Theorem 4.2.

(Scherer–Wendler, Theorem 6.1 in [26]) Let $M$ be any matrix whose column
space is equal to the span of
$[\boldsymbol{e},\mbox{\boldmath$A$}\boldsymbol{e},\mbox{\boldmath$A$}^{2}\boldsymbol{e},\ldots,\mbox{\boldmath$A$}^{s-1}\boldsymbol{e}]$.
The function $W(z)$ in (3) is $A$-stable if and only if there exists a matrix
$\mbox{\boldmath$R$}\in\mathbb{S}^{s}$ such that

$\displaystyle\begin{cases}\mbox{\boldmath$R$}\,\boldsymbol{e}=\boldsymbol{b}\,,\\\
\mbox{\boldmath$X$}=\mbox{\boldmath$R$}\mbox{\boldmath$A$}+\mbox{\boldmath$A$}^{T}\mbox{\boldmath$R$}-\boldsymbol{b}\boldsymbol{b}^{T}\,,\\\
\mbox{\boldmath$M$}^{T}\mbox{\boldmath$R$}\mbox{\boldmath$M$}\succeq 0\,,\\\
\mbox{\boldmath$M$}^{T}\mbox{\boldmath$X$}\mbox{\boldmath$M$}\succeq
0\,.\end{cases}$ (14)

We refer to (14) as the CSTW conditions.

######  Remark 4.3.

When $W(z)$ is nondegenerate, the matrix $M$ can be the identity matrix $I$.
If $W(z)$ is degenerate, then (14) with
$\mbox{\boldmath$M$}=\mbox{\boldmath$I$}$ provides a sufficient condition for
$A$-stability, but it may not be necessary. In the degenerate case, $M$ may be
chosen as
$[\boldsymbol{e},\mbox{\boldmath$A$}\boldsymbol{e},\ldots,\mbox{\boldmath$A$}^{r-1}\boldsymbol{e}]$,
where $r$ is the smallest number for which
$\mbox{\boldmath$A$}^{r}\boldsymbol{e}$ can be expressed in terms of the
vectors
$[\boldsymbol{e},\mbox{\boldmath$A$}\boldsymbol{e},\ldots,\mbox{\boldmath$A$}^{r-1}\boldsymbol{e}]$.

######  Remark 4.4 (Degenerate stability functions).

While degenerate stability functions may seem uncommon, they do appear in the
literature. Recently, Runge-Kutta schemes with degenerate stability functions
have found application in avoiding order reduction via the weak stage order
conditions (also referred to as parabolic order conditions) [5, 4, 3].
$\spadesuit$

The conditions (14) define a convex set for $R$ and are readily converted into
an LMI by parameterizing the equality constraints. Let

$\displaystyle\mbox{\boldmath$B$}:=\textrm{diag}\begin{bmatrix}b_{1},&b_{2},&\ldots,&b_{s}\end{bmatrix}\,,$

and

$\displaystyle\mbox{\boldmath$N$}_{ij}:=\boldsymbol{n}_{ij}\boldsymbol{n}_{ij}^{T}\qquad\textrm{where}\qquad\boldsymbol{n}_{ij}=\boldsymbol{e}_{i}-\boldsymbol{e}_{j}\,,$

for $1\leq i<j\leq s$ with $\boldsymbol{e}_{i}$ being the $i$th unit vector.
Then by construction, $\mbox{\boldmath$N$}_{ij}$ is a basis for the vector
space $\mbox{\boldmath$N$}_{ij}\boldsymbol{e}=0$ (cf. [26]), and $R$ has the
form:

$R$
$\displaystyle=\mbox{\boldmath$B$}+\mbox{\boldmath$N$}(\boldsymbol{\eta})\,\quad\textrm{where}\quad\mbox{\boldmath$N$}(\boldsymbol{\eta})=\sum_{i=1}^{s-1}\sum_{j=i+1}^{s}\eta_{ij}\mbox{\boldmath$N$}_{ij}\,.$
(15)

Condition (14) in LMI form reads:

$\displaystyle\begin{bmatrix}\mbox{\boldmath$B$}&\mbox{\boldmath$0$}\\\
\mbox{\boldmath$0$}&\mbox{\boldmath$B$}\mbox{\boldmath$A$}+\mbox{\boldmath$A$}^{T}\mbox{\boldmath$B$}-\boldsymbol{b}\boldsymbol{b}^{T}\end{bmatrix}+\sum_{i=1}^{s-1}\sum_{j=i+1}^{s}\eta_{ij}\begin{bmatrix}\mbox{\boldmath$N$}_{ij}&\mbox{\boldmath$0$}\\\
\mbox{\boldmath$0$}&\mbox{\boldmath$N$}_{ij}\mbox{\boldmath$A$}+\mbox{\boldmath$A$}^{T}\mbox{\boldmath$N$}_{ij}\end{bmatrix}\succeq
0\,.$ (16)

The CSTW LMI (16) and the $E$-polynomial approach in Lemma 4.1 both ensure
$A$-stability for $W(z)$. A key difference between the two approaches lies in
the incorporation of the order conditions. Satisfying the order conditions
removes lower order terms in the $E$-polynomial. In contrast, the CSTW LMI
(16) does not account for the order conditions, which turns out to be
important in practice.

###  4.3 Main Lemma: Sharpening the CSTW Conditions

The authors in [26] observe (though do not resolve) that zero eigenvalues of
$X$ may limit the practical application of the algebraic conditions (14). The
following Lemma shows that the order conditions result in $X$ always admitting
a family of zero eigenvalues whenever the CSTW LMI is feasible.

######  Lemma 4.5.

Let $(\mbox{\boldmath$A$},\boldsymbol{b})$ satisfy the order conditions (6)
with order $p\geq 2$. If $\mbox{\boldmath$R$}\in\mathbb{S}^{s}$ satisfies
$\mbox{\boldmath$R$}\boldsymbol{e}=\boldsymbol{b}$ and
$\mbox{\boldmath$X$}\succeq 0$ where

$\displaystyle\mbox{\boldmath$X$}:=\mbox{\boldmath$R$}\mbox{\boldmath$A$}+\mbox{\boldmath$A$}^{T}\mbox{\boldmath$R$}-\boldsymbol{b}\boldsymbol{b}^{T}\,,$

then $X$ has the following null vectors:

$\displaystyle\mbox{\boldmath$X$}\mbox{\boldmath$A$}^{j-1}\boldsymbol{e}=0\qquad\textrm{for}\qquad
1\leq j\leq\left\lfloor\frac{p}{2}\right\rfloor\,.$

The proof utilizes the fact that for any matrix $\mbox{\boldmath$X$}\succeq
0$,

$\displaystyle\textrm{if}\qquad\boldsymbol{v}^{\,T}\mbox{\boldmath$X$}\boldsymbol{v}=0\,,\qquad\textrm{then}\qquad\mbox{\boldmath$X$}\boldsymbol{v}=0\,.$

For instance, the Cholesky factorization
$\mbox{\boldmath$X$}=\mbox{\boldmath$Q$}^{T}\mbox{\boldmath$Q$}$ shows that
$\boldsymbol{v}^{T}\mbox{\boldmath$X$}\boldsymbol{v}=\|\mbox{\boldmath$Q$}\boldsymbol{v}\|^{2}=0$.

###### Proof.

To simplify notation, let

$\displaystyle\boldsymbol{v}_{j}=\mbox{\boldmath$A$}^{j}\boldsymbol{e}\qquad
j\geq 0\,.$

We show that for all

$\displaystyle 0\leq n\leq\left\lfloor\frac{p}{2}\right\rfloor-1\,,$

the expression

$\displaystyle\boldsymbol{v}_{n}^{T}\mbox{\boldmath$X$}\boldsymbol{v}_{n}=0\,.$

Since $\mbox{\boldmath$X$}\succeq 0$, it follows that
$\mbox{\boldmath$X$}\boldsymbol{v}_{n}=0$.

For $n=0$, the tall-tree conditions are $\boldsymbol{b}^{T}\boldsymbol{e}=1$
and $\boldsymbol{b}^{T}\boldsymbol{v}_{1}=\tfrac{1}{2}$, hence:

$\displaystyle\boldsymbol{v}_{0}^{T}\mbox{\boldmath$X$}\boldsymbol{v}_{0}$
$\displaystyle=\boldsymbol{b}^{T}\boldsymbol{v}_{1}+\boldsymbol{v}_{1}^{T}\boldsymbol{b}-(\boldsymbol{e}^{T}\boldsymbol{b})^{2}$
$\displaystyle=0\,.$

We now proceed by strong induction: Let $n\leq\lfloor\tfrac{p}{2}\rfloor-1$ be
any positive integer, and assume that
$\mbox{\boldmath$X$}\boldsymbol{v}_{k}=0$ for all $0\leq k\leq n-1$. We show
that $\boldsymbol{v}_{n}^{T}\mbox{\boldmath$X$}\boldsymbol{v}_{n}=0$, which
then completes the proof.

By hypothesis and the $p$th order tall-tree conditions, we have that

$\displaystyle\boldsymbol{b}^{T}\boldsymbol{v}_{j}=\frac{1}{(j+1)!}\qquad\textrm{for}\qquad
j=0,\ldots,2n+1\,,$ (17)

as $2n+1\leq p-1$ by the choice of $n$.

The first step is to obtain an expression for
$\mbox{\boldmath$R$}\boldsymbol{v}_{n}$. Substituting the definition of $X$ in
terms of $R$ into the induction hypothesis
$\mbox{\boldmath$X$}\boldsymbol{v}_{k}=0$, yields the recursion relation

$\mbox{\boldmath$R$}\boldsymbol{v}_{k+1}=\boldsymbol{b}\boldsymbol{b}^{T}\boldsymbol{v}_{k}-\mbox{\boldmath$A$}^{T}\mbox{\boldmath$R$}\boldsymbol{v}_{k}\,\qquad\textrm{for}\qquad
0\leq k\leq n-1\,,$ (18)

so that $\mbox{\boldmath$R$}\boldsymbol{v}_{k+1}$ is written in terms of
$\mbox{\boldmath$R$}\boldsymbol{v}_{k}$. Setting $k=n-1$ in (18), we can use
the conditions (17) on dot products $\boldsymbol{b}^{T}\boldsymbol{v}_{k}$ and
iteratively eliminate $\mbox{\boldmath$R$}\boldsymbol{v}_{j}$ to express
$\mbox{\boldmath$R$}\boldsymbol{v}_{n}$ in the basis
$\\{\boldsymbol{b},\mbox{\boldmath$A$}^{T}\boldsymbol{b},\ldots,(\mbox{\boldmath$A$}^{T})^{n-1}\boldsymbol{b}\\}$:

$\displaystyle\mbox{\boldmath$R$}\boldsymbol{v}_{n}$
$\displaystyle=\sum_{j=0}^{n}\frac{(-1)^{j}}{(n-j)!}(\mbox{\boldmath$A$}^{T})^{j}\,\boldsymbol{b}\,.$
(19)

We then have:

$\displaystyle\boldsymbol{v}_{n}^{T}\mbox{\boldmath$X$}\boldsymbol{v}_{n}$
$\displaystyle=\boldsymbol{v}_{n}^{T}\left(\mbox{\boldmath$R$}\mbox{\boldmath$A$}+\mbox{\boldmath$A$}^{T}\mbox{\boldmath$R$}-\boldsymbol{b}\boldsymbol{b}^{T}\right)\boldsymbol{v}_{n}\,,$
$\displaystyle=2\,\boldsymbol{v}_{n+1}^{T}\mbox{\boldmath$R$}\boldsymbol{v}_{n}-\frac{1}{(n+1)!^{2}}$
$\displaystyle\textrm{(since
}\mbox{\boldmath$A$}\boldsymbol{v}_{n}=\boldsymbol{v}_{n+1}\textrm{)}\,,$
$\displaystyle=2\,\sum_{j=0}^{n}\frac{(-1)^{j}}{(n-j)!(n+j+2)!}-\frac{1}{(n+1)!^{2}}\quad$
$\displaystyle\textrm{(via \eqref{eq:talltreev2} and \eqref{eq:Rv_sum})}\,,$
$\displaystyle=\frac{(-1)^{n}}{(2n+2)!}\,\underbrace{\sum_{j=0}^{2n+2}\begin{pmatrix}2n+2\\\
j\end{pmatrix}1^{2n+2-j}(-1)^{j}}_{=(1-1)^{2n+2}}\,,$ $\displaystyle=0\,.$

∎

######  Remark 4.6.

Lemma 4.5 can be viewed as a generalization that algebraically stable methods
admit a set of null vectors in their algebraic stability matrix
$\mbox{\boldmath$B$}\mbox{\boldmath$A$}+\mbox{\boldmath$A$}^{T}\mbox{\boldmath$B$}-\boldsymbol{b}\boldsymbol{b}^{T}$,
For further reference, see the proof of [19, Lemma 13.14].$\spadesuit$

The following theorem modifies the original CSTW conditions by including the
null vectors of $X$ as additional constraints.

######  Theorem 4.7 (Main Result, CSTW Conditions with order conditions).

Let $W(z)$ be a $p$th order approximation to $e^{z}$, i.e.,
$(\mbox{\boldmath$A$},\boldsymbol{b})$ satisfy (6), and let $M$ be any matrix
whose column space is equal to the span of
$[\boldsymbol{e},\mbox{\boldmath$A$}\boldsymbol{e},\mbox{\boldmath$A$}^{2}\boldsymbol{e},\ldots,\mbox{\boldmath$A$}^{s-1}\boldsymbol{e}]$.
Then $W(z)$ is $A$-stable if and only if there exists a matrix
$\mbox{\boldmath$R$}\in\mathbb{S}^{s}$ such that

$\displaystyle\begin{cases}\mbox{\boldmath$R$}\,\boldsymbol{e}=\boldsymbol{b}\,,&\\\
\mbox{\boldmath$X$}=\mbox{\boldmath$R$}\mbox{\boldmath$A$}+\mbox{\boldmath$A$}^{T}\mbox{\boldmath$R$}-\boldsymbol{b}\boldsymbol{b}^{T}\,,&\\\
\mbox{\boldmath$X$}\mbox{\boldmath$A$}^{j-1}\boldsymbol{e}=0\,&\qquad\textrm{for}\quad
j=1,\,\ldots,\,\left\lfloor\frac{p}{2}\right\rfloor\,,\\\
\mbox{\boldmath$M$}^{T}\mbox{\boldmath$R$}\mbox{\boldmath$M$}\succeq 0\,,&\\\
\mbox{\boldmath$M$}^{T}\mbox{\boldmath$X$}\mbox{\boldmath$M$}\succeq
0\,,&\end{cases}$ (20)

###### Proof.

Note that the proof of Lemma 4.7 holds if $\mbox{\boldmath$X$}\succeq 0$ is
replaced with
$\mbox{\boldmath$M$}^{T}\mbox{\boldmath$X$}\mbox{\boldmath$M$}\succeq 0$ for
any matrix $M$ whose columns span the null vectors
$\mbox{\boldmath$A$}^{j-1}\boldsymbol{e}$ for
$j\leq\left\lfloor\frac{p}{2}\right\rfloor$. By definition, the $M$ appearing
in Theorem 4.7 contains all $\mbox{\boldmath$A$}^{j-1}\boldsymbol{e}$ for
$j\leq\left\lfloor\frac{p}{2}\right\rfloor$. Therefore, any matrix pair
$(\mbox{\boldmath$R$},\mbox{\boldmath$X$})$ satisfying the CSTW conditions
(14) also satisfy (20). ∎

Again, feasibility of (20) is a sufficient condition for $A$-stability when
$\mbox{\boldmath$M$}=\mbox{\boldmath$I$}$.

An LMI for the modified CSTW conditions (20) can be obtained by parameterizing
the $\eta_{ij}$ variables in (15) to satisfy the additional affine constraints
$\mbox{\boldmath$X$}\mbox{\boldmath$A$}^{j-1}\boldsymbol{e}=0$ for
$j=1,\ldots,\left\lfloor\frac{p}{2}\right\rfloor$. The parameterization is
then substituted back into (16).

###  4.4 Examples Highlighting Lemma 4.5

In this section, we provide examples highlighting the significance of Lemma
4.5. The examples demonstrate that without incorporating the constraints
imposed by the approximation property (6), computational approaches are
unlikely to provide rigorous certifications for $A$-stability.

Given a fixed pair $(\mbox{\boldmath$A$},\boldsymbol{b})$, the following
convex set is introduced

$\displaystyle\mathcal{R}(\mbox{\boldmath$A$},\boldsymbol{b}):=\big{\\{}\boldsymbol{\eta}\in\mathbb{R}^{s(s-1)/2}\;:\;\textrm{The
LMI}\;\eqref{LMI:CSTW}\;\textrm{holds}\big{\\}}\,,$

The (original, unmodified) CSTW conditions in Theorem 4.2 are equivalently
rephrased as: A scheme $(\mbox{\boldmath$A$},\boldsymbol{b})$ is $A$-stable if
and only if $\mathcal{R}(\mbox{\boldmath$A$},\boldsymbol{b})$ is non-empty.

It might be expected that the dimension of the convex set $\mathcal{R}$ is
$s(s-1)/2$ (i.e., the dimension of symmetric matrices minus the constraints
imposed by $\mbox{\boldmath$R$}\boldsymbol{e}=\boldsymbol{b}$). However, Lemma
4.5 indicates that the inequality constraint $\mbox{\boldmath$X$}\succeq 0$
(or $\mbox{\boldmath$M$}^{T}\mbox{\boldmath$X$}\mbox{\boldmath$M$}\succeq 0$)
further reduces the dimension of $\mathcal{R}$ to be strictly smaller.
Consequently, computational approaches that seek
$\boldsymbol{\eta}\in\mathcal{R}$ will essentially “never” find feasible
points without correctly characterizing the low dimensional space
$\mathcal{R}$ — which is provided by Theorem 4.7.

The first example demonstrates how the zero eigenvalues of $X$ reduce the
dimension of the convex set $\mathcal{R}$. The second example shows why
$\mbox{\boldmath$X$}\succeq 0$ is a necessary hypothesis in the Lemma.

######  Example 4.8 (SDIRK(3,2)).

This example constructs the set
$\mathcal{R}(\mbox{\boldmath$A$},\boldsymbol{b})$ for the following SDIRK
3-stage $p=2$ method

$\begin{array}[]{c|c}\boldsymbol{c}&\mbox{\boldmath$A$}\\\
\hline\cr&\boldsymbol{b}^{T}\end{array}=\begin{array}[]{c|ccccc}1&1&0&0\\\
\frac{3}{2}&\frac{1}{2}&1&0\\\ 1&1&-1&1\\\ \hline\cr&1&-1&1\end{array}\,.$

The set $\mathcal{R}(\mbox{\boldmath$A$},\boldsymbol{b})$ is defined as
$\boldsymbol{\eta}=\left[\eta_{12},\,\eta_{13},\,\eta_{23}\right]^{T}$ for
which $\mbox{\boldmath$R$}(\boldsymbol{\eta})\succeq 0$ and
$\mbox{\boldmath$X$}(\boldsymbol{\eta})\succeq 0$ where

$\displaystyle\mbox{\boldmath$R$}(\boldsymbol{\eta})$
$\displaystyle=\mbox{\boldmath$B$}+\eta_{12}\,\mbox{\boldmath$N$}_{12}+\eta_{13}\,\mbox{\boldmath$N$}_{13}+\eta_{23}\,\mbox{\boldmath$N$}_{23}\,,$
$\displaystyle\mbox{\boldmath$X$}(\boldsymbol{\eta})$
$\displaystyle=(\mbox{\boldmath$B$}\mbox{\boldmath$A$}+\mbox{\boldmath$A$}^{T}\mbox{\boldmath$B$}-\boldsymbol{b}\boldsymbol{b})+\eta_{12}\,(\mbox{\boldmath$N$}_{12}\mbox{\boldmath$A$}+\mbox{\boldmath$A$}^{T}\mbox{\boldmath$N$}_{12})+\eta_{13}\,(\mbox{\boldmath$N$}_{13}\mbox{\boldmath$A$}+\mbox{\boldmath$A$}^{T}\mbox{\boldmath$N$}_{13})$
$\displaystyle\quad+\eta_{23}\,(\mbox{\boldmath$N$}_{23}\mbox{\boldmath$A$}+\mbox{\boldmath$A$}^{T}\mbox{\boldmath$N$}_{13})\,,$

and $\mbox{\boldmath$B$}=\mathrm{diag}(1,\;-1,\;1)$. Given that
$\boldsymbol{e},\mbox{\boldmath$A$}\boldsymbol{e}$, and
$\mbox{\boldmath$A$}^{2}\boldsymbol{e}$ are linearly independent, the matrix
$M$ in the CSTW conditions is set to the identity matrix.

The CSTW conditions suggest the dimension of
$\mathcal{R}(\mbox{\boldmath$A$},\boldsymbol{b})$ is $3$. However, the
implication of Theorem 4.5 is that the dimension is in fact $1$. Since $p=2$,
Theorem 4.5 allows for the inclusion of the constraint
$\mbox{\boldmath$X$}\boldsymbol{e}=0$ within $\mathcal{R}$, i.e.,

$\displaystyle\mathcal{R}(\mbox{\boldmath$A$},\boldsymbol{b})=\big{\\{}\boldsymbol{\eta}\in\mathbb{R}^{3}\;:\;\mbox{\boldmath$R$}(\boldsymbol{\eta})\succeq
0\,,\;\mbox{\boldmath$X$}(\boldsymbol{\eta})\succeq
0\;,\mbox{\boldmath$X$}\boldsymbol{e}=0\big{\\}}\,.$

The constraint $\mbox{\boldmath$X$}\boldsymbol{e}=0$ imposes two independent
linear equations on $\boldsymbol{\eta}$ whose solution forces $\eta_{12}=3$
and $\eta_{32}=2$. Thus,

$\displaystyle\mathcal{R}(\mbox{\boldmath$A$},\boldsymbol{b})$
$\displaystyle=\big{\\{}\boldsymbol{\eta}\in\mathbb{R}^{3}\;:\;\eta_{12}=3\,,\;\;\eta_{32}=2\,,\;\;\mbox{\boldmath$R$}(\eta_{13})\succeq
0\,,\;\mbox{\boldmath$X$}(\eta_{13})\succeq 0\big{\\}}\,,$

where

$\displaystyle\mbox{\boldmath$R$}(\eta_{13})$
$\displaystyle=\begin{bmatrix}\phantom{-}4&-3&\phantom{-}0\\\
-3&\phantom{-}4&-2\\\
\phantom{-}0&-2&\phantom{-}3\end{bmatrix}+\eta_{13}\begin{bmatrix}\phantom{-}1&\phantom{-}0&-1\\\
\phantom{-}0&\phantom{-}0&\phantom{-}0\\\
-1&\phantom{-}0&\phantom{-}1\end{bmatrix}\,,$
$\displaystyle\mbox{\boldmath$X$}(\eta_{13})$
$\displaystyle=\begin{bmatrix}\phantom{-}4&-5&\phantom{-}1\\\
-5&\phantom{-}11&-6\\\
\phantom{-}1&-6&\phantom{-}5\end{bmatrix}+\eta_{13}\begin{bmatrix}\phantom{-}0&\phantom{-}1&-1\\\
\phantom{-}1&\phantom{-}0&-1\\\ -1&-1&\phantom{-}2\end{bmatrix}\,.$

Figure 1 visualizes non-empty $1$-dimensional set $\mathcal{R}$, which notably
has empty interior.

Figure 1: The blue line is a visualization of
$\mathcal{R}(\mbox{\boldmath$A$},\boldsymbol{b})$ for SDIRK(3,2) in Example
4.8. Note the set $\mathcal{R}$ has dimension $1$, as characterized by Theorem
4.7, which is lower than the expected $3$ dimensional set suggested by Theorem
4.2.

######  Example 4.9 (Hammer & Hollingsworth).

This example demonstrates why $\mbox{\boldmath$X$}\succeq 0$ is required as a
hypothesis in Theorem 4.5 for $X$ to have zero eigenvalues. The Hammer &
Hollingsworth two stage $p=4$ method [17, Table II.7.3] is represented as

$\begin{array}[]{c|c}\boldsymbol{c}&\mbox{\boldmath$A$}\\\
\hline\cr&\boldsymbol{b}^{T}\end{array}=\begin{array}[]{c|cc}\tfrac{1}{2}-\tfrac{\sqrt{3}}{6}&\tfrac{1}{4}&\tfrac{1}{4}-\tfrac{\sqrt{3}}{6}\\\
\tfrac{1}{2}+\tfrac{\sqrt{3}}{6}&\tfrac{1}{4}+\tfrac{\sqrt{3}}{6}&\tfrac{1}{4}\\\
\hline\cr&\tfrac{1}{2}&\tfrac{1}{2}\end{array}\,.$

The scheme yields the following stability function and $E$-polynomial

$\displaystyle
W(z)=\frac{1+z/2+z^{2}/12}{1-z/2+z^{2}/12}\,\qquad\textrm{and}\qquad\quad
E(y)=0\,.$

Since $E(y)\geq 0$, the scheme is $A$-stable and
$\mathcal{R}(\mbox{\boldmath$A$},\boldsymbol{b})$ is not empty.

The set $\mathcal{R}(\mbox{\boldmath$A$},\boldsymbol{b})$ is characterized by
$\eta\in\mathbb{R}$ via

$\displaystyle\mbox{\boldmath$R$}(\eta)$
$\displaystyle=\frac{1}{2}\mbox{\boldmath$I$}+\eta\begin{pmatrix}\phantom{-}1\\\
-1\end{pmatrix}\
\begin{pmatrix}1&-1\end{pmatrix}\,,\qquad\mbox{\boldmath$X$}(\eta)=-\frac{\sqrt{3}}{2}\eta\begin{pmatrix}1&0\\\
0&-1\end{pmatrix}\,.$

The only value of $\eta$ for which $\mbox{\boldmath$X$}(\eta)\succeq
0,\mbox{\boldmath$R$}(\eta)\succeq 0$ is $\eta=0$, i.e.,
$\mathcal{R}(\mbox{\boldmath$A$},\boldsymbol{b})=\\{0\\}$.

Note that $\mbox{\boldmath$X$}\succeq 0$ is required for $X$ to have a non-
trivial null space. For any value of $\eta\neq 0$, the matrix $X$ is
invertible and has no null vectors, even though
$(\mbox{\boldmath$A$},\boldsymbol{b})$ satisfies the order conditions. When
$\mbox{\boldmath$X$}\succeq 0$, it has two linear independent null vectors
predicted by the Theorem 4.5:
$\\{\boldsymbol{e},\mbox{\boldmath$A$}\boldsymbol{e}\\}$ — and thus must be
the zero matrix (i.e., $\eta=0$).

##  5 Examples and Results for Runge-Kutta Schemes

To provide examples of using LMIs to certify the stability of RK methods, we
start with an overview in Subsection 5.1 that details the general numerical
approach. We then apply the approach to an idealized example, the $A$-stable
SDIRK(5,4) scheme, in Subsection 5.2. Additionally, we examine $A$-stability
for several recently devised schemes developed via numerical software [2, 5].
These schemes do not satisfy the tall-tree order conditions exactly but admit
a residual of typical size $\mathcal{O}(10^{-15})$. The section concludes with
examples establishing rigorous bounds on $\alpha$ for $A(\alpha)$-stability.

###  5.1 Computational Details for Rigorous Certification

We describe computational details for rigorously verifying stability via
feasibility of an LMI, proceeding in two steps.

First, we use CVX, a package for specifying and solving convex programs [8,
14], to numerically solve one of the following semi-definite programming
problems:

  * •

To demonstrate the non-negativity of the $E$-polynomial from § 4.1, i.e.,
assess the feasibility of (12), we solve

$\begin{array}[]{ll}\textrm{Minimize: }&1\\\ \textrm{Subject to:
}&\mbox{\boldmath$F$}(\boldsymbol{\eta})=\mbox{\boldmath$P$}+\sum_{j=1}^{d}\eta_{j}\,\mbox{\boldmath$N$}_{j}\succeq
0\,,\end{array}$ (P1)

  * •

To assess the feasibility of the modified CSTW approach (20), we solve

$\begin{array}[]{ll}\textrm{Minimize: }&1\\\ \textrm{Subject to:
}&\mbox{\boldmath$R$}\,\boldsymbol{e}=\boldsymbol{b}\,,\\\
&\mbox{\boldmath$X$}=\mbox{\boldmath$R$}\mbox{\boldmath$A$}+\mbox{\boldmath$A$}^{T}\mbox{\boldmath$R$}-\boldsymbol{b}\boldsymbol{b}^{T}\,,\\\
&\mbox{\boldmath$X$}\mbox{\boldmath$A$}^{j-1}\boldsymbol{e}=0\,,\qquad\qquad\qquad\textrm{for}\quad
j=1,\,\ldots,\,\left\lfloor\frac{p}{2}\right\rfloor\,,\\\
&\mbox{\boldmath$M$}^{T}\mbox{\boldmath$R$}\mbox{\boldmath$M$}\succeq 0\,,\\\
&\mbox{\boldmath$M$}^{T}\mbox{\boldmath$X$}\mbox{\boldmath$M$}\succeq
0\,,\end{array}$ (P2)

where $M$ is any matrix containing the columns where $r$ is the largest
integer for which $\mbox{\boldmath$A$}^{r}\boldsymbol{e}$ can be written as a
linear combination of
$[\boldsymbol{e},\;\mbox{\boldmath$A$}\boldsymbol{e},\;\ldots,\mbox{\boldmath$A$}^{r-1}\boldsymbol{e}]$.
From an implementation standpoint, we convert the constraints in (P2) (as
discussed in § 3.1) into an LMI.

While not a proof, a positive output of CVX provides numerical evidence that
the convex set in question is feasible.

Second, we use the output of CVX to construct a rigorous certificate of
stability in settings where $\mbox{\boldmath$P$}\in\mathbb{Q}^{m\times m}$, or
$\mbox{\boldmath$A$}\in\mathbb{Q}^{s\times
s},\boldsymbol{b}\in\mathbb{Q}^{s}$. We convert the double-precision floating-
point outputs of CVX into symbolic rational entries and perform exact symbolic
LDL factorizations for $F$ in (P1) or
$\mbox{\boldmath$R$},\mbox{\boldmath$X$}$ in (P2) to yield matrices
$\mbox{\boldmath$L$},\mbox{\boldmath$D$}$ with rational entries. A rigorous
certificate is then ensured provided $\mbox{\boldmath$D$}\succeq 0$.

For example, in the case of (P1), CVX yields a set of
$\eta_{j}^{*}\in\mathbb{Q}$ for $(1\leq j\leq d)$; substituting
$\\{\eta_{j}^{*}\\}_{j=1}^{d}$ back into the LMI and symbolically computing an
LDL factorization yields

$\displaystyle\mbox{\boldmath$F$}(\boldsymbol{\eta}^{*})=\mbox{\boldmath$L$}\mbox{\boldmath$D$}\mbox{\boldmath$L$}^{T}\qquad\textrm{with}\qquad\mbox{\boldmath$L$},\mbox{\boldmath$D$}\in\mathbb{Q}^{m\times
m}\,.$

For presentation and simplicity, we often round the output values of CVX
$\eta_{j}^{*}$ to nearby $\bar{\eta}_{j}\in\mathbb{Q}$ with integer numerator
and denominators having fewer digits, yet still yielding positive certificates
of feasibility.

######  Remark 5.1.

Throughout this section, we present coefficient matrices wherever possible.
However, in several of the practical examples, the coefficients of various
matrices, e.g., $P$, $X$, etc., as well as the LDL factorizations of $F$ or
$\mbox{\boldmath$X$},\mbox{\boldmath$R$}$ admit rational values exceeding 100
digits. $\spadesuit$

######  Remark 5.2.

Note that if the coefficients of $P$ or $\mbox{\boldmath$A$},\boldsymbol{b}$
lie in a field extension $\mathbb{F}$ of $\mathbb{Q}$, the LDL factorization
also admits matrices $\mbox{\boldmath$L$},\mbox{\boldmath$D$}$ in
$\mathbb{F}$. Therefore, this approach generalizes in a straightforward way as
long as one can determine the sign of any element $x\in\mathbb{F}$, such as
with interval arithmetic. In addition, it is worth noting that schemes
$\mbox{\boldmath$A$},\boldsymbol{b}$ with irrational entries may have
$E$-polynomials with rational entries and coefficient matrices $P$, allowing
this approach to apply without any further modification (e.g., the scheme (32)
has coefficients in $\mathbb{Q}[\sqrt{2}]$ but $E$-polynomial with
coefficients in $\mathbb{Q}$) $\spadesuit$.

###  5.2 Example of $A$-stability in SDIRK(5,4)

The first example demonstrates the hybrid computational-analytic LMI solution
approach to verify $A$-stability for an idealized example that satisfies the
tall tree order conditions exactly and is known to be $A$-stable. The scheme
is a Diagonally Implicit RK method, represented by the following Butcher
tableau [19, Table 6.5, Chapter IV.6].

$\displaystyle\begin{array}[]{c|c}\boldsymbol{c}&\mbox{\boldmath$A$}\\\
\hline\cr&\boldsymbol{b}^{T}\end{array}\;=\begin{array}[]{c|ccccc}\frac{1}{2}&\frac{1}{4}&&&&\\\
\frac{3}{4}&\frac{1}{2}&\frac{1}{4}&&&\\\
\frac{11}{20}&\frac{17}{50}&-\frac{1}{25}&\frac{1}{4}&&\\\
\frac{1}{2}&\frac{371}{1360}&-\frac{137}{2720}&\frac{15}{544}&\frac{1}{4}&\\\
1&\frac{25}{24}&-\frac{49}{48}&\frac{125}{16}&-\frac{85}{12}&\frac{1}{4}\\\
\hline\cr&\frac{25}{24}&-\frac{49}{48}&\frac{125}{16}&-\frac{85}{12}&\frac{1}{4}\end{array}\,.$
(29)

We provide two proofs: one using the $E$-polynomial LMI and a second based on
modified CSTW LMI conditions.

####  5.2.1 E-ploynomial linear matrix inequality for SDIRK(5,4)

The SDIRK(5,4) scheme in (29) has an $E$-polynomial
$E(y)=y^{6}(9y^{4}-64y^{2}+512)$ ($\alpha=\frac{\pi}{2}$), which after
factoring out the largest monomial factor, yields

$\displaystyle
F(y):=9y^{4}-64y^{2}+512=\boldsymbol{y}^{T}\,\mbox{\boldmath$F$}(\eta)\,\boldsymbol{y}\,,$

where

$\displaystyle\mbox{\boldmath$F$}(\eta)=\begin{bmatrix}512&0&0\\\ 0&-64&0\\\
0&0&9\end{bmatrix}+\eta\begin{bmatrix}0&0&1\\\ 0&-2&0\\\
1&0&0\end{bmatrix}\,.$ (30)

For this example, the space $\mathcal{N}_{3}$ (introduced in § 3.2) has
dimension $d=1$ and is spanned by the second matrix in (30). The numerical
solution (P1) for $\mbox{\boldmath$F$}(\eta)$ yields the following output:

$\displaystyle\eta^{*}=-61.786375823904734\,.$

The fact that CVX obtains a solution is numerical evidence suggesting
$E(y)\geq 0$, indicating that the scheme is $A$-stable.

We can turn this numerical evidence into a rigorous proof by substituting
rational values of $\eta$ close to this computational output and factorizing
exactly
$\mbox{\boldmath$F$}(\eta)=\mbox{\boldmath$L$}\mbox{\boldmath$D$}\mbox{\boldmath$L$}^{T}$
with $\mbox{\boldmath$L$},\mbox{\boldmath$D$}\in\mathbb{Q}^{3\times 3}$. While
there is no requirement to substitute integer values of $\eta$, for the sake
of presentation, we round the CVX output $\eta^{*}$ to $\overline{\eta}=-60$
and obtain:

$\displaystyle\mbox{\boldmath$F$}(-60)=\begin{bmatrix}512&0&-60\\\ 0&56&0\\\
-60&0&9\end{bmatrix}=\mbox{\boldmath$L$}\,\mbox{\boldmath$D$}\,\mbox{\boldmath$L$}^{T}\,,$

where

$\displaystyle\mbox{\boldmath$L$}=\begin{bmatrix}1&0&0\\\\[2.0pt]
0&1&0\\\\[2.0pt]
\text{-}\frac{15}{128}&0&1\end{bmatrix}\,,\qquad\textrm{and}\qquad\mbox{\boldmath$D$}=\begin{bmatrix}512&0&0\\\\[2.0pt]
0&56&0\\\\[2.0pt] 0&0&\frac{63}{32}\end{bmatrix}\succeq 0\,.$

This results in the following SOS representation of $E(y)$, certifying
$A$-stability:

$\displaystyle
E(y)=y^{6}\,\left(512\left(\frac{15}{128}y^{2}-1\right)^{2}+56y^{2}+\frac{63}{32}y^{4}\right)\,.$

######  Remark 5.3.

Since the linear constraints are exactly parameterized, there is robustness in
the choice of $\eta$. For example, the value of $\bar{\eta}=-32$ similarly
yields a certificate of $A$-stability since
$\mbox{\boldmath$F$}(-32)=\mbox{\boldmath$L$}\,\mbox{\boldmath$D$}\,\mbox{\boldmath$L$}^{T}\succeq
0$, where

$L$ $\displaystyle=\left[\begin{array}[]{ccc}1&0&0\\\ 0&1&0\\\
\text{-}\frac{1}{16}&0&1\end{array}\right]\,,\qquad\textrm{and}\qquad\mbox{\boldmath$D$}=\left[\begin{array}[]{ccc}512&0&0\\\
0&0&0\\\ 0&0&7\end{array}\right]\succeq 0\,.$

$\spadesuit$

####  5.2.2 The Modified CSTW-Method for SDIRK(5,4)

Using the SDIRK(5,4) scheme in (29) as an example, this section highlights
several key differences in the modified CSTW approach relative to the
$E$-polynomial.

Since $p=4$, Theorem 4.5 adds two linear constraints on $X$ (corresponding to
null vectors $\boldsymbol{e}$ and $\mbox{\boldmath$A$}\boldsymbol{e}$) into
the associated LMI. Consequently, problem (P2) has three degrees of freedom,
denoted by the vector $\boldsymbol{\eta}$. The rounded numerical solution of
(P2) yields matrices

$\begin{aligned}
\mbox{\boldmath$X$}^{*}&=\left[\begin{array}[]{ccccc}\frac{729823}{97920}&\text{-}\frac{348733}{195840}&\frac{875727}{21760}&\text{-}\frac{334871}{7200}&\frac{237}{400}\\\
\text{-}\frac{348733}{195840}&\frac{170083}{391680}&\text{-}\frac{1259867}{130560}&\frac{160397}{14400}&\text{-}\frac{57}{400}\\\
\frac{875727}{21760}&\text{-}\frac{1259867}{130560}&\frac{5678645}{26112}&\text{-}\frac{241217}{960}&\frac{16}{5}\\\
\text{-}\frac{334871}{7200}&\frac{160397}{14400}&\text{-}\frac{241217}{960}&\frac{1045211}{3600}&\text{-}\frac{1479}{400}\\\
\frac{237}{400}&\text{-}\frac{57}{400}&\frac{16}{5}&\text{-}\frac{1479}{400}&\frac{19}{400}\end{array}\right]\\\
\end{aligned}\quad\textrm{and}\quad\begin{aligned}
\mbox{\boldmath$R$}^{*}&=\left[\begin{array}[]{ccccc}\frac{195061}{16320}&\text{-}\frac{42157}{10880}&\frac{416905}{6528}&-72&\frac{11}{10}\\\
\text{-}\frac{42157}{10880}&\frac{131641}{65280}&\text{-}\frac{324335}{13056}&\frac{1259}{48}&\text{-}\frac{11}{20}\\\
\frac{416905}{6528}&\text{-}\frac{324335}{13056}&\frac{4888637}{13056}&\text{-}\frac{99107}{240}&\frac{73}{10}\\\
-72&\frac{1259}{48}&\text{-}\frac{99107}{240}&\frac{34459}{75}&\text{-}\frac{391}{50}\\\
\frac{11}{10}&\text{-}\frac{11}{20}&\frac{73}{10}&\text{-}\frac{391}{50}&\frac{11}{50}\end{array}\right],\end{aligned}$

which then admits LDL factorizations of the form

$\displaystyle\mbox{\boldmath$X$}^{*}=\mbox{\boldmath$L$}_{X}\mbox{\boldmath$D$}_{X}\mbox{\boldmath$L$}_{X}^{T}\,,\qquad\qquad\mbox{\boldmath$R$}^{*}=\mbox{\boldmath$L$}_{R}\mbox{\boldmath$D$}_{R}\mbox{\boldmath$L$}_{R}^{T}\,,$

where the matrices have coefficients in $\mathbb{Q}$ with

$\displaystyle\mbox{\boldmath$D$}_{X}=\operatorname{diag}\left[\begin{array}[]{c}\frac{729823}{97920}\\\
\frac{2466451}{280252032}\\\ \frac{7352143}{246645100}\\\ 0\\\
0\end{array}\right]\,,\qquad\mbox{\boldmath$D$}_{R}=\operatorname{diag}\left[\begin{array}[]{c}\frac{195061}{16320}\\\
\frac{28479739}{37451712}\\\ \frac{3647946461}{341756868}\\\
\frac{3800443925}{43775357532}\\\ \frac{103805}{2104052}\end{array}\right]\,.$

The fact that $\mbox{\boldmath$D$}_{X}$ admits two zero eigenvalues is by
construction and follows from Theorem 4.7. The certification of $A$-stability
is thus confirmed as the pair
$\mbox{\boldmath$X$}^{*},\mbox{\boldmath$R$}^{*}$ satisfy the CSTW condition
in exact arithmetic.

###  5.3 Schemes Failing to Satisfy The Tall-Tree Order Conditions

Building on the idealized example from the previous subsection, we now focus
on certifying stability for RK schemes developed from numerical solutions of
the order conditions (i.e., using numerical optimization software). The
schemes considered here have coefficients that do not exactly satisfy the
order conditions (6). Instead, the coefficients satisfy these conditions with
a small residual. Failing to satisfy the tall-tree order conditions leads to
two complications: the $E$-polynomial no longer satisfies (7) but instead
includes low-degree monomials with small coefficients; additionally, the
matrix $X$, used in the CSTW conditions, no longer admits exact zero
eigenvalues, but rather small nonzero ones.

We test schemes from two sources:

  * •

Diagonally implicit Runge-Kutta schemes with weak stage order [5] (cf. [4,
3]). These schemes were developed to alleviate the effects of _order
reduction_ on stiff problems, primarily arising from spatial discretizations
of linear partial differential equations. The schemes are denoted as WSO
DIRK(s,p,q), where $s$, $p$, and $q$ denote the number of stages, classical
order, and weak stage order, respectively.

  * •

(Very high order) Diagonally implicit Runge-Kutta schemes with additional
practical properties, developed in [2].

It is worth noting that while both [5, 2] provide strong numerical evidence
for $A$-stability (e.g., a solution with small residual was provided for the
CSTW LMI in [2]) neither work provides a rigorous certificate in the form of
an exact solution to one of the LMI’s.

We adopt two strategies for testing $A$-stability, outlined as follows:

Strategy 1: The Butcher coefficients $(\mbox{\boldmath$A$},\boldsymbol{b})$
are reported as decimal expansions, typically to 16 digits of accuracy (e.g.,
the accuracy of double-precision floating-point arithmetic from which they
were obtained). Treating the coefficients as rational values, we test for the
non-negativity of the associated $E$-polynomial, which also has rational
coefficients. This strategy determines whether the dynamics (2) are $A$-stable
despite the coefficients not exactly satisfying the tall-tree order
conditions.

In contrast to the $E$-polynomial approach, challenges arise when using the
CSTW approach to obtain a rigorous certificate of $A$-stability. The modified
CSTW conditions depend on the exact satisfaction of the tall tree order
conditions. However, in cases where the Butcher coefficients approximate these
conditions with a small residual, the null vectors of $X$ are no longer
applicable constraints in the SDP, and the original CSTW conditions must be
used. The feasible set of the CSTW conditions may then take the form of a
tubular domain of a low-dimensional set, thereby introducing computational
challenges to the numerical solution.

Strategy 2: To overcome the challenges in Strategy 1 due to approximations in
the tall tree order conditions, we assess the $A$-stability of perturbed
schemes $(\mbox{\boldmath$\tilde{A}$},\boldsymbol{\tilde{b}})$ that
simultaneously:

  * •

Have coefficients in $\mathbb{Q}$ and satisfy the tall-tree order conditions
(6) in exact arithmetic;

  * •

Are perturbations of the reported scheme
$(\mbox{\boldmath$A$},\boldsymbol{b})$ in the literature, satisfying an error
bound

$\displaystyle|b_{i}-\tilde{b}_{i}|<\epsilon_{b}\quad\textrm{and}\quad|a_{ij}-{\tilde{a}}_{ij}|<\epsilon_{A}\quad\textrm{for}\quad
i,j=1,\ldots,s\,.$ (31)

Since the perturbed schemes,
$(\mbox{\boldmath$\tilde{A}$},\boldsymbol{\tilde{b}})$, satisfy the tall-tree
order conditions, their $E$-polynomials satisfy (7) and the associated $X$
matrix in the CSTW approach admits null vectors (in exact arithmetic)
characterized by Theorem 4.5. Thus, both the $E$-polynomial and CSTW
approaches provide pathways for certifying $A$-stability for
$(\mbox{\boldmath$\tilde{A}$},\boldsymbol{\tilde{b}})$.

Several schemes that fail to achieve rigorous certificates of $A$-stability
under Strategy 1 are shown to be $(\epsilon_{A},\epsilon_{b})$-close to
schemes that attain rigorous $A$-stability certificates under Strategy 2.

####  5.3.1 Certification of $A$-stability via Strategy 1

Here, we verify $A$-stability of three schemes using Strategy 1. Since the
schemes do not satisfy the order conditions exactly, the $E$-polynomial
($\alpha=\frac{\pi}{2}$) has the form:

$\displaystyle{E}(y)=y^{2}\boldsymbol{y}^{T}\mbox{\boldmath$F$}(\eta)\boldsymbol{y},\qquad\textrm{with}\qquad\mbox{\boldmath$F$}(\eta)=\mbox{\boldmath$P$}+\sum_{j=1}^{d}\eta_{j}\mbox{\boldmath$N$}_{j}\,,$

where $P$ and $\mbox{\boldmath$N$}_{j}$ are as in § 3.2. The details for each
scheme are as follows:

  1. 1.

WSO DIRK(12,5,4) developed in [5]: A symbolic computation of $E(y)$ yields a
diagonal matrix $P$, containing the coefficients of $E(y)$. Since
$\mbox{\boldmath$P$}\succeq 0$, it follows that
$\mbox{\boldmath$F$}(\eta)\succeq 0$ trivially when $\boldsymbol{\eta}=0$.
Consequently, no SDP is required, as the LMI is immediately satisfied.

  2. 2.

WSO DIRK(7,4,4), also developed in [5]: A symbolic computation of
${E}(y)=y^{2}\boldsymbol{y}^{T}\mbox{\boldmath$P$}\boldsymbol{y}$, yields a
diagonal matrix $\mbox{\boldmath$P$}\in\mathbb{Q}^{7\times 7}$ with two
negative coefficients

$\displaystyle p_{6}<0\,,\qquad p_{8}<0\,.$

Since $\mbox{\boldmath$F$}(0)$ does not satisfy the LMI, we seek solutions to
the LMI (which has dimension $d=15$) via an SDP. The $E$-polynomial SDP
identifies a candidate feasible solution $\boldsymbol{\eta^{*}}$, which for
computational simplicity, we round to $\boldsymbol{\bar{\eta}}$, resulting in
the matrix:

$\displaystyle\mbox{\boldmath$F$}(\boldsymbol{\bar{\eta}})=\mbox{\boldmath$P$}-108420\mbox{\boldmath$N$}_{10}+20\mbox{\boldmath$N$}_{12}-3420\mbox{\boldmath$N$}_{13}-30\mbox{\boldmath$N$}_{15}=\mbox{\boldmath$L$}\mbox{\boldmath$D$}\mbox{\boldmath$L$}^{T},$

with $\mbox{\boldmath$L$},\mbox{\boldmath$D$}\in\mathbb{Q}^{7\times 7}$ and
$\mbox{\boldmath$D$}\succ 0$.

  3. 3.

DIRK(13,8)[1]A[(14,6)A] developed in [2]: Similar to the previous case,
symbolic computation of the $E$-polynomial yields
$\mbox{\boldmath$P$}\in\mathbb{Q}^{13\times 13}$ with two negative
coefficients

$\displaystyle p_{14}<0\,,\qquad p_{22}<0\,,$

requiring an SDP to find a potential sum of squares representation for $E(y)$.
The $E$-polynomial LMI has dimension $d=66$, and the corresponding SDP
identifies a feasible solution $\boldsymbol{\eta^{*}}$, which upon rounding
yields $\boldsymbol{\bar{\eta}}$ (reported in Appendix A.1), resulting in the
matrix
$\mbox{\boldmath$F$}(\boldsymbol{\bar{\eta}})=\mbox{\boldmath$L$}\mbox{\boldmath$D$}\mbox{\boldmath$L$}^{T}$,
where the diagonal matrix $\mbox{\boldmath$D$}\succ 0$.

The complete set of coefficients for
$\boldsymbol{\bar{\eta}},\mbox{\boldmath$L$},\mbox{\boldmath$D$}$, as well as
generating Matlab code, can be found in the supplemental material.

The results are formalized with the following Lemma.

######  Lemma 5.4.

The $E$-polynomials for the weak stage order schemes DIRK(12,5,4) and
DIRK(7,4,4) in [5], and DIRK(13,8)[1]A[(14,6)A] in [2] are non-negative, and
the associated schemes are $A$-stable.

####  5.3.2 Certification of $A$-stability via Strategy 2

We now apply Strategy 2 to certify $A$-stability for schemes where Strategy 1
fails due to the $E$-polynomial for $(\mbox{\boldmath$A$},\boldsymbol{b})$
being negative near the origin.

Throughout, tildes denote quantities for the perturbed scheme
$(\mbox{\boldmath$\tilde{A}$},\boldsymbol{\tilde{b}})$, e.g., $\tilde{P}$ is
the matrix containing the coefficients of the $E$-polynomial $\tilde{E}(y)$.

  1. 1.

DIRK(6,6)[1]A[(7,5)A] developed in [2]: The original scheme’s E-polynomial,
${E}(y)=y^{2}\boldsymbol{y}^{T}\mbox{\boldmath$P$}\boldsymbol{y}$, possesses
three negative coefficients in $P$

$\displaystyle p_{0}<0,\qquad p_{2}<0,\qquad p_{4}<0\,.$

Given that $p_{0}<0$, $E(y)$ is negative for values of $y$ near the origin,
indicating that the scheme is not $A$-stable. This inability to achieve
$A$-stability directly follows from the failure to satisfy the tall-tree order
conditions — which, if satisfied, would ensure $p_{0}=p_{2}=p_{4}=0$.

We introduce the perturbed scheme

$\displaystyle\mbox{\boldmath$\tilde{A}$}$
$\displaystyle=\left[\begin{array}[]{cccccc}\frac{33128226}{109158329}&&&&&\\\
\text{-}\frac{254432096}{909477001}&\frac{51289103}{102571593}&&&&\\\
\frac{33289838}{118645151}&\text{-}\frac{825218320}{1881654059}&\frac{130993323}{602959172}&&&\\\
\text{-}\frac{13583292}{200438515}&\frac{156154430}{158643099}&\text{-}\frac{65409371}{245235917}&\frac{81765600}{330141853}&&\\\
\frac{16354062}{130133299}&\text{-}\frac{247816720}{248961507}&\frac{169383005}{222482121}&\text{-}\frac{49241043}{234166886}&\frac{79900588}{92184791}&\\\
\text{-}\frac{34719176}{94331171}&\text{-}\frac{155737141}{155748342}&\frac{42945649}{80312134}&\text{-}\frac{50573402}{289227347}&\frac{678237381}{1102812170}&\frac{31879369}{45767530}\end{array}\right],$

where the coefficient vector $\boldsymbol{\tilde{b}}\in\mathbb{Q}^{6}$ are
chosen to exactly satisfy the tall-tree conditions (6) with $p=6$. The pair
$(\mbox{\boldmath$\tilde{A}$},\boldsymbol{\tilde{b}})$ satisfy the error bound
(31) with $\epsilon_{A}=5\cdot 10^{-17}$, and $\epsilon_{b}=6\cdot 10^{-15}$.

$A$-stability for the perturbed scheme is certified via two approaches:

     * •

The $E$-polynomial
$\tilde{E}(y)=y^{8}\boldsymbol{y}^{T}\mbox{\boldmath$\tilde{P}$}\boldsymbol{y}\geq
0$ follows immediately since
$\mbox{\boldmath$\tilde{P}$}\in\mathbb{Q}^{3\times 3}$ satisfies
$\mbox{\boldmath$\tilde{P}$}\succeq 0$.

     * •

The modified CSTW approach (P2) yields an LMI in 3 variables. The numerical
solution produces a solution pair:

$\displaystyle\tilde{\mbox{\boldmath$X$}}^{*}=\tilde{\mbox{\boldmath$L$}}_{X}\tilde{\mbox{\boldmath$D$}}_{X}\tilde{\mbox{\boldmath$L$}}^{T}_{X}\quad\textrm{and}\quad\tilde{\mbox{\boldmath$R$}}^{*}=\tilde{\mbox{\boldmath$L$}}_{R}\tilde{\mbox{\boldmath$D$}}_{R}\tilde{\mbox{\boldmath$L$}}^{T}_{R}\,,$

with $\tilde{\mbox{\boldmath$D$}}_{X}\succeq 0$
($\tilde{\mbox{\boldmath$D$}}_{X}$ has $4$ zero eigenvalues by Theorem 4.7),
and $\tilde{\mbox{\boldmath$D$}}_{R}\succeq 0$.

  2. 2.

SDIRK(9,6)[1]SAL[(9,5)A] developed in [2]: Similar to the previous example,
this scheme is not $A$-stable as the $E$-polynomial is negative near the
origin. This is reflected in the matrix
$\mbox{\boldmath$P$}\in\mathbb{Q}^{9\times 9}$ having negative coefficients

$\displaystyle p_{0}<0,\qquad p_{2}<0,\qquad p_{4}<0\,,$

a consequence of failing to satisfy the tall-tree order conditions. A
perturbed scheme $(\mbox{\boldmath$\tilde{A}$},\boldsymbol{\tilde{b}})$,
defined in Appendix A.2, satisfies the error bound (31) with
$\epsilon_{A}=\epsilon_{b}=8\cdot 10^{-15}$.

Similar to DIRK(6,6), $A$-stability for the perturbed scheme
$(\mbox{\boldmath$\tilde{A}$},\boldsymbol{\tilde{b}})$ is certified through
two approaches:

     * •

The $E$-polynomial
$\tilde{E}(y)=y^{8}\boldsymbol{y}^{T}\tilde{\mbox{\boldmath$P$}}\boldsymbol{y}\geq
0$, with $\tilde{\mbox{\boldmath$P$}}\in\mathbb{Q}^{6\times 6}$, is
immediately non-negative since $\tilde{\mbox{\boldmath$P$}}\succeq 0$,
avoiding the need for solving an SDP.

     * •

The modified CSTW approach, which utilizes an LMI with $15$-degree-of-freedom,
is solved via SDP. The numerical solver identifies optimal pairs
$(\tilde{\mbox{\boldmath$R$}}^{*},\tilde{\mbox{\boldmath$X$}}^{*})$, yielding
exact LDL factorizations over $\mathbb{Q}$, with all matrices $D$ being
positive semidefinite.

  3. 3.

WSO DIRK(12,5,5) developed in [5]: Similar to the previous two examples, this
scheme also fails to be $A$-stable due to the lowest order terms in the
$E$-polynomial having negative coefficients

$\displaystyle p_{0}<0,\qquad p_{2}<0,$

implying that $E(y)<0$ for small $y$. A perturbed scheme
$(\mbox{\boldmath$\tilde{A}$},\boldsymbol{\tilde{b}})$ is defined in Appendix
A.2 and satisfies (31) with $\epsilon_{A}=\epsilon_{b}=9\cdot 10^{-15}$.

For $(\mbox{\boldmath$\tilde{A}$},\mbox{\boldmath$\tilde{b}$})$, the
$E$-polynomial has the form
$\tilde{E}(y)=y^{6}\boldsymbol{y}^{T}\mbox{\boldmath$\tilde{F}$}(\boldsymbol{\eta})\boldsymbol{y}$,
where $\mbox{\boldmath$\tilde{P}$}\in\mathbb{Q}^{10\times 10}$ is not positive
definite. The associated LMI for the non-negativity of $\tilde{F}$ contains
$\dim\mathcal{N}_{10}=36$ degrees of freedom. The Appendix A.2 presents a
solution $\boldsymbol{\bar{\eta}}$ that yields
$\mbox{\boldmath$\tilde{F}$}(\boldsymbol{\bar{\eta}})=\mbox{\boldmath$\tilde{L}$}\mbox{\boldmath$\tilde{D}$}\mbox{\boldmath$\tilde{L}$}^{T}$,
with $\mbox{\boldmath$\tilde{D}$}\succeq 0$, thereby certifying $A$-stability.

We formalize the result in the following Lemma.

######  Lemma 5.5.

There exist perturbed schemes
$(\mbox{\boldmath$\tilde{A}$},\boldsymbol{\tilde{b}})$ for
DIRK(6,6)[1]A[(7,5)A] and
SDIRK(9,6)[1]SAL[(9,5)A] from [2], and DIRK(12,5,5) from [5], which satisfy
the error bound (31) with $\epsilon_{A},\epsilon_{b}=\mathcal{O}(10^{-15})$.
These schemes simultaneously satisfy the tall-tree conditions and are
$A$-stable.

###  5.4 Examples Establishing $A(\alpha)$ Stability

We now shift attention and establish rigorous bounds on $\alpha$, for
$A(\alpha)$-stability of two schemes that are not $A$-stable.

First note that $E(y^{2};\beta)$, where $\beta:=\cos(\alpha)$, is a bivariate
polynomial in $y,\beta$ with rational coefficients. An initial rational value
of $\beta$ is fixed, and the associated $E$-polynomial SDP (P1) is solved to
determine whether the associated feasible set is nonempty. The value of
$\beta$ is incrementally decreased until (P1) indicates that the feasible set
is empty. The last $\beta$ value before the feasible set is reported non-empty
is $\beta^{*}$. The corresponding $\alpha^{*}=\cos^{-1}(\beta^{*})$ is the
lower bound for the maximal angle $\alpha$ at which the scheme is confirmed to
be $A(\alpha)$-stable. An exact rational certificate in the form of an LDL
factorization is then reported.

The result is formalized as follows:

######  Lemma 5.6.

The schemes (32) and (35) are $A(\alpha)$-stable for some maximum angle
$\alpha$ that is bounded below by the angle $\alpha^{*}$, where
$\alpha^{*}=\cos^{-1}\beta^{*}$ and $\beta^{*}$ is defined in (34) and (36)
respectively.

####  5.4.1 The IRK(4,4) Scheme of Ramos and Vigo

The following scheme, developed by Ramos and Vigo, is a $4$-stage, $4$th order
fully implicit method based on a BDF-type Chebyshev approximation [24]:

$\displaystyle\mathbf{A}=\begin{bmatrix}\frac{22-\sqrt{2}}{96}&\frac{5-8\sqrt{2}}{48}&\frac{22-7\sqrt{2}}{96}&\frac{-1}{16}\\\\[4.0pt]
\frac{4+3\sqrt{2}}{24}&\frac{1}{6}&\frac{4-3\sqrt{2}}{24}&0\\\\[4.0pt]
\frac{22+7\sqrt{2}}{96}&\frac{5+8\sqrt{2}}{48}&\frac{22+\sqrt{2}}{96}&\frac{-1}{16}\\\\[4.0pt]
\frac{1}{3}&\frac{1}{3}&\frac{1}{3}&0\end{bmatrix}\quad\textrm{and}\quad\mathbf{b}=\begin{bmatrix}\frac{1}{3}&\frac{1}{3}&\frac{1}{3}&0\end{bmatrix}^{T}.$
(32)

By defining $\beta:=\cos(\alpha)$, the generalized $E$-polynomial (4) becomes
the following bivariate polynomial in $y$ and $\beta$:

$\displaystyle E(y^{2};\beta)=y^{2}\Big{(}$ $\displaystyle y^{14}+22\beta
y^{12}+(272\beta^{2}-16)y^{10}+(1920\beta^{3}+96\beta)y^{8}$
$\displaystyle+(6144\beta^{4}+3840\beta^{2})y^{6}+(36864\beta^{3}+10752\beta)y^{4}+73728\beta
y^{2}+294912\beta\Big{)}\,.$ (33)

When $\beta=0$ (corresponding to $\alpha=\frac{\pi}{2}$), the $E$-polynomial

$\displaystyle E\left(y;\frac{\pi}{2}\right)=y^{6}\,(y^{2}-16)\,,$

is negative for $|y|<4$, and the scheme is not $A$-stable.

Instead, we determine a lower bound for the maximal value of $\alpha$ for
which $E(y^{2};\beta)\geq 0$. The last $\beta$ value before reaching an empty
feasible set is the bound

$\displaystyle\beta^{*}=\frac{19699132}{4466212691}\qquad(\alpha^{*}\approx
89.74728^{o})\,.$ (34)

Solving Problem (P1) with the value $\beta^{*}$ (details provided in Appendix
B.1) produces the positive certificate

$\displaystyle
E(y^{2};\beta^{*})=y^{2}\boldsymbol{y}^{T}\mbox{\boldmath$F$}(\boldsymbol{\bar{\eta}})\boldsymbol{y}\geq
0$

where
$\mbox{\boldmath$F$}(\boldsymbol{\bar{\eta}})=\mbox{\boldmath$L$}\mbox{\boldmath$D$}\mbox{\boldmath$L$}^{T}$,
and the diagonal matrix $\mbox{\boldmath$D$}\succ 0$, certifying the scheme is
$A(\alpha)$-stable for maximal angle $\alpha\geq\alpha^{*}$.

####  5.4.2 The ESDIRK(8,6) Scheme of Skvortsov

We examine the 8-stage, 6th order ESDIRK scheme of Skvortsov [27]
characterized by its Butcher matrix:

$\mbox{\boldmath$A$}=\left[\begin{array}[]{ccccccccc}0&&&&&&&&\\\
\frac{1}{6}&\frac{1}{6}&&&&&&&\\\
\frac{11}{96}&\text{-}\frac{1}{32}&\frac{1}{6}&&&&&&\\\
\frac{1}{12}&\text{-}\frac{1}{4}&\frac{1}{2}&\frac{1}{6}&&&&&\\\
\text{-}\frac{2015}{15072}&\text{-}\frac{6987}{5024}&\frac{3271}{1884}&\frac{175}{471}&\frac{1}{6}&&&&\\\
\text{-}\frac{326531}{573678}&\text{-}\frac{114988}{31871}&\frac{1208156}{286839}&\frac{132950}{286839}&\frac{68}{203}&\frac{1}{6}&&&\\\
\text{-}\frac{331717945}{2106545616}&\text{-}\frac{480525599}{416107776}&\frac{2240951089}{1404363744}&\frac{394951619}{2808727488}&\text{-}\frac{5160553}{26834976}&\frac{35815}{352512}&\frac{1}{6}&&\\\
\frac{16264655341}{73026914688}&\frac{9786099235}{14425069568}&\text{-}\frac{34306812733}{48684609792}&\text{-}\frac{15985588007}{97369219584}&\frac{37652437}{930279168}&\text{-}\frac{340747}{12220416}&\frac{1}{26}&\frac{1}{6}&\\\
\frac{7}{90}&0&0&0&\frac{16}{45}&\text{-}\frac{4}{45}&\frac{2}{15}&\frac{16}{45}&\frac{1}{6}\end{array}\right]$
(35)

The scheme is stiffly accurate, with $\boldsymbol{b}^{T}$ being the last row
of $A$.

The generalized $E$-polynomial, in $y$ and $\beta$, takes the form:

$\displaystyle E(y^{2},\beta)$
$\displaystyle=y^{2}\sum_{j=0}^{15}q_{2j}(\beta)\,y^{2j}\,\qquad\textrm{where}\qquad\beta=\cos(\alpha)\,.$

Appendix B.2 presents formulas for the polynomials $q_{2j}(\beta)$.

Similar to the scheme of Ramos and Vigo, this scheme fails to be $A$-stable.
When $\beta=0$ ($\alpha=\frac{\pi}{2}$), the $E$-polynomial simplifies to:

$\displaystyle
E\left(y;\frac{\pi}{2}\right)=y^{8}\left(y^{8}+\tfrac{61415271}{616225}y^{6}-\tfrac{91554624}{3925}y^{4}-\tfrac{3175034112}{3925}y^{2}-\tfrac{4152010752}{785}\right)\,.$

Since the lowest order monomial term is negative, e.g., $-4152010752/785<0$,
the polynomial $E(y;\pi/2)$ is negative for values of $y$ near the origin.

After iterating through values of $\beta$, we establish the bound:

$\displaystyle\beta^{*}=\frac{2218472195}{100000000000}\qquad(\alpha^{*}\approx
88.7288^{o})\,.$ (36)

This value is consistent with the $\alpha$ bound reported in [27] and is now
accompanied by a rigorous certificate.

The $E$-polynomial produced by $\beta^{*}$ can be written as

$\displaystyle
E(y^{2};\beta^{*})=y^{2}\boldsymbol{y}^{T}\mbox{\boldmath$F$}(\boldsymbol{\eta})\boldsymbol{y}\,,\quad\textrm{with}\quad\mbox{\boldmath$F$}(\boldsymbol{\eta})=\mbox{\boldmath$P$}+\sum_{j=1}^{105}\eta_{j}\mbox{\boldmath$N$}_{j}\,,$

where $\mbox{\boldmath$P$}\in\mathbb{Q}^{16\times 16}$ is a diagonal matrix
containing three negative coefficients and $d=105$ is the dimension (11).
Numerically solving Problem (P1) yields a solution $\boldsymbol{\eta}^{*}$,
which again we round for simplicity to a nearby rational value
$\bar{\boldsymbol{\eta}}$ presented in Appendix B.2. Substituting
$\bar{\boldsymbol{\eta}}$ into the LMI and factorizing yields

$\displaystyle
E(y^{2};\beta^{*})=y^{2}\boldsymbol{y}^{T}\mbox{\boldmath$L$}\mbox{\boldmath$D$}\mbox{\boldmath$L$}^{T}\boldsymbol{y}\geq
0\,,$ (37)

where $\mbox{\boldmath$D$},\mbox{\boldmath$L$}\in\mathbb{Q}^{16\times 16}$
with $D$ being diagonal with non-negative entries.

######  Remark 5.7.

Even though $\boldsymbol{\eta}^{*}$ has been rounded to a nearby rational
value $\bar{\boldsymbol{\eta}}$ with fewer digits, which we report in the
Appendix, the matrix and $D$ in (37) contains integer denominators with up to
548 digits. $\spadesuit$

##  6 Discussion and Conclusions

This paper adopts a semidefinite optimization approach to rigorously certify
$A$\- and $A(\alpha)$\- stability in Runge-Kutta schemes and improves the
algebraic conditions for $A$-stability by incorporating the order conditions.
A key application of this work is the validation and construction of stable
time-stepping schemes for potential implementation in industrial software.
Looking ahead, these approaches and perspectives could also be useful in
certifying stability in other time-integration schemes. Semidefinite
programming can be used to certify stability in settings where algebraic
conditions for stability are known, such as the algebraic stability of general
linear methods. On the other hand, there are settings, such as $A$-stability
in general linear methods, where, to our knowledge, algebraic conditions for
stability have yet to be formulated, providing opportunities for future work.

## Acknowledgments

The authors would like to thank David Ketcheson for helpful discussions. This
material is based upon work supported by the National Science Foundation under
Grants No. DMS–2012268 and DMS–2309727.

## Data Availability

Certificates of stability and the supporting numerical code are available
online from

https://github.com/ajuhl/Algebraic-Conditions-for-Stability-Certified-via-SDP.

## References

  * [1] Amir Ali Ahmadi and Raphaël M. Jungers  “Switched stability of nonlinear systems via SOS-convex Lyapunov functions and semidefinite programming”  In _52nd IEEE Conference on Decision and Control_ , 2013, pp. 727–732
  * [2] Y. Alamri and D. Ketcheson  “Very high-order A-stable stiffly accurate diagonally implicit Runge-Kutta methods with error estimators” arXiv: 2211.14574, 2024
  * [3] A. Biswas, D. Ketcheson, S. Roberts, B. Seibold and D. Shirokoff  “Explicit Runge Kutta Methods that Alleviate Order Reduction” arXiv:2310.02817, 2023
  * [4] A. Biswas, D. Ketcheson, B. Seibold and D. Shirokoff  “Algebraic Structure of the Weak Stage Order Conditions for Runge–Kutta Methods”  In _SIAM Journal on Numerical Analysis_ 62.1, 2024, pp. 48–72  DOI: 10.1137/22M1483943
  * [5] A. Biswas, D. Ketcheson, B. Seibold and D. Shirokoff  “Design of DIRK schemes with high weak stage order”  In _Commun. Appl. Math. Comput._ 18, 2023, pp. 1–28
  * [6] J.C. Butcher  “On A-stable implicit Runge-Kutta methods”  In _BIT Numerical Mathematics_ 17, 1977, pp. 375–378
  * [7] G. Cooper  “An algebraic condition for A-stable Runge-Kutta methods”  In _Numerical Analysis_ Longman Scientific & Technical, 1986, pp. 32–46
  * [8] Inc. CVX Research  “CVX: Matlab Software for Disciplined Convex Programming, version 2.0”, http://cvxr.com/cvx, 2012
  * [9] G. Dahlquist  “A special stability problem for linear multistep methods”  In _BIT Numerical Mathematics_ 3, 1963, pp. 27–43
  * [10] G. Dahlquist  “Positive Functions and some Applications to Stability Questions for Numerical Methods”  In _Recent Advances in Numerical Analysis_ Academic Press, 1978, pp. 1–29
  * [11] Germund Dahlquist  “On stability and error analysis for stiff non-linear problems part i”, 1975
  * [12] Y. Drori and M. Teboulle  “Performance of first-order methods for smooth convex minimization: a novel approach”  In _Mathematical Programming_ 145, 2014, pp. 1436–4646  DOI: 10.1007/s10107-013-0653-0
  * [13] A. Gahlawat and G. Valmorbida  “A semi-definite programming approach to stability analysis of linear partial differential equations”  In _2017 IEEE 56th Annual Conference on Decision and Control (CDC)_ , 2017, pp. 1882–1887
  * [14] M. Grant and S. Boyd  “Graph implementations for nonsmooth convex programs”  In _Recent Advances in Learning and Control_ , Lecture Notes in Control and Information Sciences  Springer-Verlag Limited, 2008, pp. 95–110
  * [15] B. Grimmer  “Provably Faster Gradient Descent via Long Steps” arxiv.org/abs/2307.06324, 2024
  * [16] E. Hairer  “Constructive characterization of A-stable approximations to exp z and its connection with algebraically stable Runge-Kutta methods”  In _Numer. Math._ 39, 1982, pp. 247–258
  * [17] E. Hairer, S.. Nørsett and G. Wanner  “Solving Ordinary Differential Equations I (2nd Revised. Ed.): Nonstiff Problems”  Springer-Verlag, New York, 1993
  * [18] E. Hairer and G. Wanner  “Algebraically stable and implementable Runge-Kutta methods of high order”  In _SIAM J. Numer. Anal._ 18, 1981, pp. 1098–1108
  * [19] E. Hairer and G. Wanner  “Solving Ordinary Differential Equations II(2nd Revised. Ed.): Stiff and Differential-Algebraic Problems”  Springer Berlin Heidelberg, 1996
  * [20] M. Jin and J. Lavaei  “Stability-Certified Reinforcement Learning: A Control-Theoretic Perspective”  In _IEEE Access_ 8, 2020, pp. 229086
  * [21] J.. Lasserre  “Moments, Positiv Polynomials and Their Applications”  Imperial College Press, 2010
  * [22] Anirudha Majumdar, Amir Ali Ahmadi and Russ Tedrake  “Control design along trajectories with sums of squares programming”  In _2013 IEEE International Conference on Robotics and Automation_ , 2013, pp. 4054–4061
  * [23] S.P. Nørsett  “C-polynomials for rational approximations to the exponential function”  In _Numer. Math._ 25, 1975, pp. 39–56
  * [24] H. Ramos and J. Vigo-Aguiar  “A fourth-order Runge-Kutta method based on BDF-type Chebyshev approximations”  In _Journal of Computational and Applied Mathematics_ 204 Elsevier, 2007, pp. 124–136
  * [25] R. Scherer and H. Turke  “Algebraic Characterization of A-stable Runge-Kutta Methods”  In _Applied Numerical Mathematics_ 5, 1989, pp. 133–144
  * [26] R. Scherer and W. Wendler  “Complete Algebraic Characterization of A-stable Runge-Kutta Methods”  In _SIAM Journal on Numerical Analysis_ 31.2 Society for IndustrialApplied Mathematics, 1994, pp. 540–551
  * [27] L.M. Skortsov  “Diagonally implicit Runge-Kutta methods for stiff problems”  In _Comput. Math. and Math. Phys._ 46, 2006, pp. 2110–2123
  * [28] G. Wanner  “Characterization of allA-stable methods of order 2m-4”  In _BIT Numerical Mathematics_ 20, 1980, pp. 367–374
  * [29] David D. Yao, Shuzhong Zhang and Xun Yu Zhou  “Stochastic Linear-Quadratic Control via Semidefinite Programming”  In _SIAM Journal on Control and Optimization_ 40.3, 2001, pp. 801–823

##  Appendix Appendix A: Supplemental Details for $A$-stable Schemes

Here, we present numerical coefficients verifying $A$-stability.

###  A.1 DIRK (13,8)[1]A[(14,6)]A

The scheme has $E$-polynomial
$E(y)=y^{2}\boldsymbol{y}^{T}\mbox{\boldmath$F$}(\boldsymbol{\bar{\eta}})\boldsymbol{y}$
where

$\displaystyle\mbox{\boldmath$F$}(\boldsymbol{\bar{\eta}})=\mbox{\boldmath$P$}+\sum_{j=1}^{66}\overline{\eta}_{j}\mbox{\boldmath$N$}_{j}\succeq
0\,,$ (38)

with coefficients of $\boldsymbol{\bar{\eta}}$ given by:

$\displaystyle\overline{\eta}_{1}$ $\displaystyle=-8470700$
$\displaystyle\qquad\overline{\eta}_{23}$ $\displaystyle=-2465700$
$\displaystyle\qquad\overline{\eta}_{45}$ $\displaystyle=\phantom{-}16400$
$\displaystyle\overline{\eta}_{2}$ $\displaystyle=\phantom{-}0$
$\displaystyle\qquad\overline{\eta}_{24}$
$\displaystyle=\phantom{-}665333207306300$
$\displaystyle\qquad\overline{\eta}_{46}$
$\displaystyle=\phantom{-}4362777200$ $\displaystyle\overline{\eta}_{3}$
$\displaystyle=-3700$ $\displaystyle\qquad\overline{\eta}_{25}$
$\displaystyle=-639200$ $\displaystyle\qquad\overline{\eta}_{47}$
$\displaystyle=\phantom{-}0$ $\displaystyle\overline{\eta}_{4}$
$\displaystyle=-17219137300$ $\displaystyle\qquad\overline{\eta}_{26}$
$\displaystyle=-27185350128356800$ $\displaystyle\qquad\overline{\eta}_{48}$
$\displaystyle=\phantom{-}2843900$ $\displaystyle\overline{\eta}_{5}$
$\displaystyle=-2500$ $\displaystyle\qquad\overline{\eta}_{27}$
$\displaystyle=-58914800$ $\displaystyle\qquad\overline{\eta}_{49}$
$\displaystyle=-1564541946300$ $\displaystyle\overline{\eta}_{6}$
$\displaystyle=\phantom{-}1552662793500$
$\displaystyle\qquad\overline{\eta}_{28}$
$\displaystyle=\phantom{-}540480365078400$
$\displaystyle\qquad\overline{\eta}_{50}$ $\displaystyle=-4500$
$\displaystyle\overline{\eta}_{7}$ $\displaystyle=\phantom{-}1100$
$\displaystyle\qquad\overline{\eta}_{29}$ $\displaystyle=-3060809665186300$
$\displaystyle\qquad\overline{\eta}_{51}$ $\displaystyle=\phantom{-}338036200$
$\displaystyle\overline{\eta}_{8}$ $\displaystyle=-291771102600$
$\displaystyle\qquad\overline{\eta}_{30}$ $\displaystyle=\phantom{-}933800$
$\displaystyle\qquad\overline{\eta}_{52}$
$\displaystyle=\phantom{-}26403742783600$ $\displaystyle\overline{\eta}_{9}$
$\displaystyle=\phantom{-}3000$ $\displaystyle\qquad\overline{\eta}_{31}$
$\displaystyle=-101765843690800$ $\displaystyle\qquad\overline{\eta}_{53}$
$\displaystyle=-61100$ $\displaystyle\overline{\eta}_{10}$
$\displaystyle=\phantom{-}5801136200$
$\displaystyle\qquad\overline{\eta}_{32}$ $\displaystyle=\phantom{-}0$
$\displaystyle\qquad\overline{\eta}_{54}$ $\displaystyle=-20078540400$
$\displaystyle\overline{\eta}_{11}$ $\displaystyle=-10101540029400$
$\displaystyle\qquad\overline{\eta}_{33}$ $\displaystyle=-16794400$
$\displaystyle\qquad\overline{\eta}_{55}$ $\displaystyle=\phantom{-}0$
$\displaystyle\overline{\eta}_{12}$ $\displaystyle=-949900$
$\displaystyle\qquad\overline{\eta}_{34}$ $\displaystyle=\phantom{-}0$
$\displaystyle\qquad\overline{\eta}_{56}$ $\displaystyle=-256500$
$\displaystyle\overline{\eta}_{13}$ $\displaystyle=\phantom{-}45726880755500$
$\displaystyle\qquad\overline{\eta}_{35}$ $\displaystyle=\phantom{-}3600$
$\displaystyle\qquad\overline{\eta}_{57}$
$\displaystyle=\phantom{-}294637386500$ $\displaystyle\overline{\eta}_{14}$
$\displaystyle=\phantom{-}244100$ $\displaystyle\qquad\overline{\eta}_{36}$
$\displaystyle=\phantom{-}19168235300$
$\displaystyle\qquad\overline{\eta}_{58}$ $\displaystyle=\phantom{-}600$
$\displaystyle\overline{\eta}_{15}$ $\displaystyle=-2222445564500$
$\displaystyle\qquad\overline{\eta}_{37}$ $\displaystyle=\phantom{-}0$
$\displaystyle\qquad\overline{\eta}_{59}$ $\displaystyle=-63663500$
$\displaystyle\overline{\eta}_{16}$ $\displaystyle=-3000$
$\displaystyle\qquad\overline{\eta}_{38}$ $\displaystyle=-14572100$
$\displaystyle\qquad\overline{\eta}_{60}$ $\displaystyle=-1286411886200$
$\displaystyle\overline{\eta}_{17}$ $\displaystyle=-1556643190544700$
$\displaystyle\qquad\overline{\eta}_{39}$ $\displaystyle=\phantom{-}0$
$\displaystyle\qquad\overline{\eta}_{61}$ $\displaystyle=\phantom{-}200$
$\displaystyle\overline{\eta}_{18}$ $\displaystyle=\phantom{-}1744900$
$\displaystyle\qquad\overline{\eta}_{40}$ $\displaystyle=\phantom{-}18000$
$\displaystyle\qquad\overline{\eta}_{62}$ $\displaystyle=\phantom{-}978576800$
$\displaystyle\overline{\eta}_{19}$ $\displaystyle=\phantom{-}292506767260200$
$\displaystyle\qquad\overline{\eta}_{41}$
$\displaystyle=\phantom{-}16836266900$
$\displaystyle\qquad\overline{\eta}_{63}$ $\displaystyle=\phantom{-}0$
$\displaystyle\overline{\eta}_{20}$ $\displaystyle=-1976300$
$\displaystyle\qquad\overline{\eta}_{42}$ $\displaystyle=\phantom{-}0$
$\displaystyle\qquad\overline{\eta}_{64}$ $\displaystyle=-5871439100$
$\displaystyle\overline{\eta}_{21}$ $\displaystyle=-5815692149200$
$\displaystyle\qquad\overline{\eta}_{43}$ $\displaystyle=-3637800$
$\displaystyle\qquad\overline{\eta}_{65}$ $\displaystyle=\phantom{-}0$
$\displaystyle\overline{\eta}_{22}$ $\displaystyle=-13689040894555700$
$\displaystyle\qquad\overline{\eta}_{44}$ $\displaystyle=-5738476103900$
$\displaystyle\qquad\overline{\eta}_{66}$ $\displaystyle=\phantom{-}1269200$

###  A.2 The Perturbed Schemes from § 5.3.2

In this appendix, we report the coefficients of the perturbed schemes
$(\mbox{\boldmath$\tilde{A}$},\boldsymbol{\tilde{b}})$ for which the lemmas in
§ 5.3.2 were established.

SDIRK(9,6)[1]SAL[(9,5)A]: The perturbed scheme coefficients are

$\displaystyle\begin{aligned}
\mbox{\boldmath$\tilde{A}$}&=\left[\begin{array}[]{ccccccccc}\frac{87518253}{401224696}&&&&&&&&\\\
\text{-}\frac{109147862}{1208036163}&\frac{87518253}{401224696}&&&&&&&\\\
\frac{70447391}{407323275}&\text{-}\frac{60128027}{170018875}&\frac{87518253}{401224696}&&&&&&\\\
\frac{258011928}{503929669}&\frac{32776647}{1131632696}&\text{-}\frac{13636148}{946751265}&\frac{87518253}{401224696}&&&&&\\\
\frac{2139251}{459753907}&\text{-}\frac{18361775}{242765601}&\frac{13889605}{63926963}&\text{-}\frac{17789601}{861400843}&\frac{87518253}{401224696}&&&&\\\
\frac{48479320}{54097599}&\frac{222681723}{1598951647}&\text{-}\frac{6683180}{35754039}&\frac{14834219}{220428796}&\text{-}\frac{67840169}{193336343}&\frac{87518253}{401224696}&&&\\\
\frac{67080581}{121311880}&\text{-}\frac{401007739}{912707597}&\frac{169517869}{507988720}&\text{-}\frac{10552416}{310889555}&\text{-}\frac{54396621}{357996284}&\frac{12212839}{571158716}&\frac{87518253}{401224696}&&\\\
\frac{176730716}{279920507}&\frac{185311713}{255696311}&\text{-}\frac{40314204}{93283073}&\frac{92464054}{154464243}&\text{-}\frac{281536253}{397040384}&\text{-}\frac{15560941}{32151589}&\frac{170501635}{450595763}&\frac{87518253}{401224696}&\\\
0&\tilde{b}_{2}&\tilde{b}_{3}&\tilde{b}_{4}&\tilde{b}_{5}&\tilde{b}_{6}&\tilde{b}_{7}&\tilde{b}_{8}&\frac{87518253}{401224696}\end{array}\right]\end{aligned}$

where $\tilde{b}_{j}$ ($3\leq j\leq 8$) is provided by solving for the tall
tree conditions (6) (for $p=6$), leaving a free variable $\tilde{b}_{2}$ which
is selected to minimize the $\ell_{2}$ difference
$\boldsymbol{b}-\boldsymbol{\tilde{b}}$.

WSO DIRK(12,5,5) The perturbed scheme coefficients are

$\mbox{\boldmath$\tilde{A}$}={\left[\begin{array}[]{cccccccccccc}\frac{10747729}{261281103}&&&&&&&&&&&\\\
\frac{39255733}{244819013}&\frac{19264472}{289086473}&&&&&&&&&&\\\
\text{-}\frac{63980223}{186836899}&\frac{70601555}{81544818}&\frac{48699329}{492234648}&&&&&&&&&\\\
\frac{757656751}{80284215}&\text{-}\frac{5452955845}{500830197}&\frac{261637874}{98954371}&\frac{86994158}{471217857}&&&&&&&&\\\
\text{-}\frac{290059410}{846787661}&\frac{142027951}{274596637}&\frac{63720572}{69536691}&\frac{8719292}{166871841}&\frac{230505997}{1977768146}&&&&&&&\\\
\text{-}\frac{96530823}{46089059}&\frac{278501442}{108044467}&\frac{186959114}{327745145}&\frac{31757051}{261668409}&\text{-}\frac{65608216}{138056009}&\frac{50822223}{96152122}&&&&&&\\\
\frac{110645970}{326232259}&\text{-}\frac{52198210}{186593643}&\frac{48069176}{46243347}&\frac{36028733}{602611029}&\text{-}\frac{42230947}{197997752}&\frac{11234921}{134641567}&\frac{339062341}{1406835502}&&&&&\\\
\frac{432515173}{73254485}&\frac{1037694284}{327224921}&\text{-}\frac{1241126819}{100347987}&\text{-}\frac{76295713}{152911958}&\frac{155195477}{71832145}&\frac{268359096}{140054533}&\frac{309243168}{155550259}&\frac{266682747}{1194765721}&&&&\\\
\frac{131193951}{284188360}&\text{-}\frac{74904386}{387416395}&\text{-}\frac{136922649}{1129220324}&\frac{46345993}{695639065}&\frac{61016892}{143403385}&\frac{74333617}{94618599}&\frac{157610940}{188314681}&\frac{31846751}{198449271}&\frac{305893355}{845914548}&&&\\\
\text{-}\frac{82024283}{115728139}&\frac{91013047}{140744863}&\frac{63214225}{132835881}&\text{-}\frac{178521351}{694495505}&\frac{96297873}{85736426}&\frac{97069417}{174996915}&\frac{25465272}{79767817}&\frac{11384038}{31516593}&\frac{232759857}{396742103}&\frac{83426711}{354434193}&&\\\
\frac{100055236}{234642175}&\frac{213090564}{161088509}&\frac{96986289}{228435568}&\text{-}\frac{174464145}{68947184}&\text{-}\frac{67238506}{859605737}&\frac{2068579853}{1961737581}&\frac{170367931}{366730407}&\frac{198092237}{172991608}&\frac{312091183}{725567705}&\frac{159786147}{106558690}&\frac{9565123}{660600961}&\\\
\frac{5486027}{454369097}&\frac{9757227}{18810635}&\frac{62278071}{555407431}&\text{-}\frac{1694527}{341651848}&\text{-}\frac{133046372}{98916929}&\frac{323864293}{952870301}&\frac{261379716}{320347663}&\tilde{b}_{8}&\tilde{b}_{9}&\tilde{b}_{10}&\tilde{b}_{11}&\tilde{b}_{12}\end{array}\right]}$

The coefficients $\tilde{b}_{j}\in\mathbb{Q}$ ($8\leq j\leq 12$) are chosen to
satisfy the order conditions exactly (6) for $p=5$.

The $E$-polynomial
$E(y)=y^{2}\boldsymbol{y}^{T}\mbox{\boldmath$F$}(\boldsymbol{\overline{\eta}})\boldsymbol{y}=\mbox{\boldmath$L$}\mbox{\boldmath$D$}\mbox{\boldmath$L$}^{T}\geq
0$ where $\boldsymbol{\overline{\eta}}\cdot 10^{-3}$ has coefficients:

$\displaystyle\overline{\eta}_{1}$ $\displaystyle=-6922561820555$
$\displaystyle\qquad\overline{\eta}_{7}$ $\displaystyle=\phantom{-}3934894$
$\displaystyle\qquad\overline{\eta}_{13}$ $\displaystyle=-550$
$\displaystyle\qquad\overline{\eta}_{19}$ $\displaystyle=-40414145977$
$\displaystyle\qquad\overline{\eta}_{25}$ $\displaystyle=\phantom{-}692417$
$\displaystyle\qquad\overline{\eta}_{31}$ $\displaystyle=-10$
$\displaystyle\overline{\eta}_{2}$ $\displaystyle=-6159041$
$\displaystyle\overline{\eta}_{8}$ $\displaystyle=-6775128853059$
$\displaystyle\overline{\eta}_{14}$ $\displaystyle=-1017731680875$
$\displaystyle\overline{\eta}_{20}$ $\displaystyle=-57578$
$\displaystyle\overline{\eta}_{26}$ $\displaystyle=-2928211$
$\displaystyle\overline{\eta}_{32}$ $\displaystyle=-40848$
$\displaystyle\overline{\eta}_{3}$ $\displaystyle=\phantom{-}401988060958$
$\displaystyle\overline{\eta}_{9}$ $\displaystyle=-4571676$
$\displaystyle\overline{\eta}_{15}$ $\displaystyle=-1687736$
$\displaystyle\overline{\eta}_{21}$ $\displaystyle=\phantom{-}139245638$
$\displaystyle\overline{\eta}_{27}$ $\displaystyle=-2$
$\displaystyle\overline{\eta}_{33}$ $\displaystyle=\phantom{-}0$
$\displaystyle\overline{\eta}_{4}$ $\displaystyle=\phantom{-}209390$
$\displaystyle\overline{\eta}_{10}$ $\displaystyle=\phantom{-}142233029108$
$\displaystyle\overline{\eta}_{16}$ $\displaystyle=\phantom{-}8353405937$
$\displaystyle\overline{\eta}_{22}$ $\displaystyle=\phantom{-}147$
$\displaystyle\overline{\eta}_{28}$ $\displaystyle=-5711$
$\displaystyle\overline{\eta}_{34}$ $\displaystyle=\phantom{-}859$
$\displaystyle\overline{\eta}_{5}$ $\displaystyle=-3199255341$
$\displaystyle\overline{\eta}_{11}$ $\displaystyle=\phantom{-}198557$
$\displaystyle\overline{\eta}_{17}$ $\displaystyle=\phantom{-}14765$
$\displaystyle\overline{\eta}_{23}$ $\displaystyle=-553788073$
$\displaystyle\overline{\eta}_{29}$ $\displaystyle=\phantom{-}3$
$\displaystyle\overline{\eta}_{35}$ $\displaystyle=\phantom{-}0$
$\displaystyle\overline{\eta}_{6}$ $\displaystyle=-4095$
$\displaystyle\overline{\eta}_{12}$ $\displaystyle=-490076976$
$\displaystyle\overline{\eta}_{18}$ $\displaystyle=-10433650$
$\displaystyle\overline{\eta}_{24}$ $\displaystyle=-636$
$\displaystyle\overline{\eta}_{30}$ $\displaystyle=\phantom{-}143788$
$\displaystyle\overline{\eta}_{36}$ $\displaystyle=-3$

##  Appendix Appendix B: Supplemental Details for $A(\alpha)$-stable Schemes

###  B.1 The IRK(4,4) Scheme of Ramos and Vigo

The scheme has $E$-polynomial
$E(y^{2},\beta^{*})=y^{2}\boldsymbol{y}^{T}\mbox{\boldmath$F$}(\boldsymbol{\bar{\eta}})\boldsymbol{y}$
where

$\displaystyle\mbox{\boldmath$F$}(\boldsymbol{\bar{\eta}})=\mbox{\boldmath$P$}+\sum_{j=1}^{21}\overline{\eta}_{j}\mbox{\boldmath$N$}_{j}\succeq
0\,,$ (39)

with coefficients of $\boldsymbol{\bar{\eta}}$ given by
$\overline{\eta}_{2}=\overline{\eta}_{4}=\overline{\eta}_{6}=\overline{\eta}_{8}=\overline{\eta}_{10}=\overline{\eta}_{13}=\overline{\eta}_{15}=\overline{\eta}_{17}=\overline{\eta}_{20}=0$
and:

$\displaystyle\overline{\eta}_{1}$ $\displaystyle=-\tfrac{343818785}{387257}$
$\displaystyle\qquad\overline{\eta}_{7}$
$\displaystyle=-\tfrac{1002782638}{963823}$
$\displaystyle\qquad\overline{\eta}_{12}$
$\displaystyle=-\tfrac{140623753}{190944}$
$\displaystyle\qquad\overline{\eta}_{18}$
$\displaystyle=\phantom{-}\tfrac{1407711}{121108}$
$\displaystyle\overline{\eta}_{3}$
$\displaystyle=\phantom{-}\tfrac{44352332}{270307}$
$\displaystyle\overline{\eta}_{9}$
$\displaystyle=\phantom{-}\tfrac{26195675}{165379}$
$\displaystyle\overline{\eta}_{14}$
$\displaystyle=\phantom{-}\tfrac{19205029}{233487}$
$\displaystyle\overline{\eta}_{19}$ $\displaystyle=-\tfrac{23169437}{338293}$
$\displaystyle\overline{\eta}_{5}$ $\displaystyle=-\tfrac{7044484}{1620291}$
$\displaystyle\overline{\eta}_{11}$ $\displaystyle=-\tfrac{526928}{268115}$
$\displaystyle\overline{\eta}_{16}$ $\displaystyle=-\tfrac{55546025}{187031}$
$\displaystyle\overline{\eta}_{21}$ $\displaystyle=-\tfrac{2727674}{410745}$

###  B.2 ESDIRK(8,6) Skvortsov Scheme in §5.4.2

We present coefficients for the scheme in subsection 5.4.2. The generalized
$E$-polynomial is

$\displaystyle E(y,\beta)$
$\displaystyle=y^{2}\sum_{j=0}^{15}q_{2j}(\beta)\,y^{2j}\,\qquad\textrm{where}\qquad\beta=\cos(\alpha)\,,$

with polynomial coefficients of $\beta$:

$\displaystyle q_{30}(\beta)$ $\displaystyle=1$ $\displaystyle q_{28}(\beta)$
$\displaystyle=96\beta\vphantom{\beta^{1}}$ $\displaystyle q_{26}(\beta)$
$\displaystyle=4032\beta^{2}+\tfrac{61415271}{616225}$ $\displaystyle
q_{24}(\beta)$ $\displaystyle=96768\beta^{3}+\tfrac{73235232}{3925}\beta$
$\displaystyle q_{22}(\beta)$
$\displaystyle=1451520\beta^{4}+\tfrac{3567255552}{3925}\beta^{2}-\tfrac{91554624}{3925}$
$\displaystyle q_{20}(\beta)$
$\displaystyle=13934592\beta^{5}+\tfrac{14120096256}{785}\beta^{3}+\tfrac{9673437312}{3925}\beta$
$\displaystyle q_{18}(\beta)$
$\displaystyle=83607552\beta^{6}+\tfrac{158718486528}{785}\beta^{4}+\tfrac{70172863488}{785}\beta^{2}-\tfrac{3175034112}{3925}$
$\displaystyle q_{16}(\beta)$
$\displaystyle=286654464\beta^{7}+\tfrac{1149206704128}{785}\beta^{5}+\tfrac{985309774848}{785}\beta^{3}+\tfrac{136916137728}{785}\beta$
$\displaystyle q_{14}(\beta)$
$\displaystyle=429981696\beta^{8}+\tfrac{5111514906624}{785}\beta^{6}+\tfrac{8045011279872}{785}\beta^{4}+\tfrac{694870576128}{157}\beta^{2}-\tfrac{4152010752}{785}$
$\displaystyle q_{12}(\beta)$
$\displaystyle=\tfrac{10416951558144}{785}\beta^{7}+\tfrac{41311620145152}{785}\beta^{5}+\tfrac{34328744275968}{785}\beta^{3}+\tfrac{5012232804864}{785}\beta$
$\displaystyle q_{10}(\beta)$
$\displaystyle=\tfrac{650132324352}{5}\beta^{6}+232190115840\beta^{4}+121899810816\beta^{2}$
$\displaystyle q_{8}(\beta)$
$\displaystyle=\tfrac{3064909529088}{5}\beta^{5}+835884417024\beta^{3}+121899810816\beta$
$\displaystyle q_{6}(\beta)$
$\displaystyle=2298682146816\beta^{4}+1671768834048\beta^{2}$ $\displaystyle
q_{4}(\beta)$ $\displaystyle=6269133127680\beta^{3}+1253826625536\beta$
$\displaystyle q_{2}(\beta)$ $\displaystyle=9403699691520\beta^{2}$
$\displaystyle q_{0}(\beta)$ $\displaystyle=5642219814912\beta\,.$

Solution coefficients $\boldsymbol{\bar{\eta}}$ for non-negativity of
$E(y,\beta^{*})$.

$\displaystyle\overline{\eta}_{1}$ $\displaystyle=-\tfrac{544417542815}{1496}$
$\displaystyle\qquad\overline{\eta}_{28}$
$\displaystyle=-\tfrac{12225979229715}{323}$
$\displaystyle\qquad\overline{\eta}_{55}$
$\displaystyle=-\tfrac{3023415560208}{431}$
$\displaystyle\qquad\overline{\eta}_{82}$
$\displaystyle=\phantom{-}\tfrac{2813400132854}{117}$
$\displaystyle\overline{\eta}_{3}$
$\displaystyle=\phantom{-}\tfrac{29842486335}{2687}$
$\displaystyle\overline{\eta}_{30}$
$\displaystyle=\phantom{-}\tfrac{433008073700}{13}$
$\displaystyle\overline{\eta}_{57}$
$\displaystyle=\phantom{-}\tfrac{667361119609}{1153}$
$\displaystyle\overline{\eta}_{84}$ $\displaystyle=-\tfrac{958434399311}{491}$
$\displaystyle\overline{\eta}_{5}$ $\displaystyle=-\tfrac{49740795}{931}$
$\displaystyle\overline{\eta}_{32}$ $\displaystyle=-\tfrac{8424611000191}{98}$
$\displaystyle\overline{\eta}_{59}$ $\displaystyle=-\tfrac{15606394145}{932}$
$\displaystyle\overline{\eta}_{85}$
$\displaystyle=\phantom{-}\tfrac{11434535023}{204}$
$\displaystyle\overline{\eta}_{7}$ $\displaystyle=-\tfrac{28918182864}{1151}$
$\displaystyle\overline{\eta}_{34}$
$\displaystyle=-\tfrac{25273367886229}{395}$
$\displaystyle\overline{\eta}_{61}$
$\displaystyle=\phantom{-}\tfrac{70764854}{907}$
$\displaystyle\overline{\eta}_{87}$ $\displaystyle=-\tfrac{220302170}{863}$
$\displaystyle\overline{\eta}_{9}$
$\displaystyle=\phantom{-}\tfrac{699763151}{1893}$
$\displaystyle\overline{\eta}_{36}$
$\displaystyle=\phantom{-}\tfrac{9262966388511}{292}$
$\displaystyle\overline{\eta}_{63}$
$\displaystyle=-\tfrac{6003091166555}{1293}$
$\displaystyle\overline{\eta}_{89}$
$\displaystyle=\phantom{-}\tfrac{9715693027109}{1173}$
$\displaystyle\overline{\eta}_{11}$ $\displaystyle=-\tfrac{611431235}{606}$
$\displaystyle\overline{\eta}_{38}$
$\displaystyle=-\tfrac{11180085898009}{131}$
$\displaystyle\overline{\eta}_{65}$
$\displaystyle=\phantom{-}\tfrac{145331883768}{577}$
$\displaystyle\overline{\eta}_{91}$ $\displaystyle=-\tfrac{261194717165}{607}$
$\displaystyle\overline{\eta}_{13}$
$\displaystyle=\phantom{-}\tfrac{4436247}{898}$
$\displaystyle\overline{\eta}_{40}$
$\displaystyle=-\tfrac{10519983461815}{799}$
$\displaystyle\overline{\eta}_{67}$ $\displaystyle=-\tfrac{2769061905}{757}$
$\displaystyle\overline{\eta}_{93}$
$\displaystyle=\phantom{-}\tfrac{5649757805}{927}$
$\displaystyle\overline{\eta}_{15}$ $\displaystyle=-\tfrac{20553023}{989}$
$\displaystyle\overline{\eta}_{42}$
$\displaystyle=-\tfrac{13505692581791}{106}$
$\displaystyle\overline{\eta}_{69}$
$\displaystyle=\phantom{-}\tfrac{2674611184492}{349}$
$\displaystyle\overline{\eta}_{95}$
$\displaystyle=-\tfrac{5175464369937}{284}$
$\displaystyle\overline{\eta}_{17}$ $\displaystyle=-\tfrac{184449}{1142}$
$\displaystyle\overline{\eta}_{44}$
$\displaystyle=\phantom{-}\tfrac{13635317782915}{474}$
$\displaystyle\overline{\eta}_{70}$ $\displaystyle=-\tfrac{435375081998}{497}$
$\displaystyle\overline{\eta}_{96}$
$\displaystyle=\phantom{-}\tfrac{1361933288513}{864}$
$\displaystyle\overline{\eta}_{19}$ $\displaystyle=-\tfrac{5040480296101}{64}$
$\displaystyle\overline{\eta}_{46}$
$\displaystyle=-\tfrac{22015444860842}{179}$
$\displaystyle\overline{\eta}_{72}$
$\displaystyle=\phantom{-}\tfrac{28638357917}{941}$
$\displaystyle\overline{\eta}_{98}$
$\displaystyle=-\tfrac{161915093932}{3421}$
$\displaystyle\overline{\eta}_{21}$
$\displaystyle=-\tfrac{20253725215145}{857}$
$\displaystyle\overline{\eta}_{48}$
$\displaystyle=-\tfrac{23644956443647}{379}$
$\displaystyle\overline{\eta}_{74}$ $\displaystyle=-\tfrac{18964563}{119}$
$\displaystyle\overline{\eta}_{100}$
$\displaystyle=\phantom{-}\tfrac{405541643}{1810}$
$\displaystyle\overline{\eta}_{23}$
$\displaystyle=\phantom{-}\tfrac{240883590679}{566}$
$\displaystyle\overline{\eta}_{50}$
$\displaystyle=-\tfrac{2159073125275}{1597}$
$\displaystyle\overline{\eta}_{76}$ $\displaystyle=-\tfrac{883878192398}{307}$
$\displaystyle\overline{\eta}_{102}$
$\displaystyle=-\tfrac{4892528411838}{1501}$
$\displaystyle\overline{\eta}_{25}$
$\displaystyle=\phantom{-}\tfrac{4696001199805}{754}$
$\displaystyle\overline{\eta}_{51}$
$\displaystyle=\phantom{-}\tfrac{58606387358}{645}$
$\displaystyle\overline{\eta}_{78}$
$\displaystyle=\phantom{-}\tfrac{38657713081}{326}$
$\displaystyle\overline{\eta}_{103}$
$\displaystyle=\phantom{-}\tfrac{127826367589}{733}$
$\displaystyle\overline{\eta}_{27}$
$\displaystyle=-\tfrac{26257051251763}{262}$
$\displaystyle\overline{\eta}_{53}$ $\displaystyle=-\tfrac{2209952042}{1439}$
$\displaystyle\overline{\eta}_{80}$ $\displaystyle=-\tfrac{1721374819}{1269}$
$\displaystyle\overline{\eta}_{105}$ $\displaystyle=-\tfrac{1622486899}{642}$

All $\overline{\eta}_{j}=0$ for any $j=1,\ldots,105$ not defined above.