# Algebraic Structure of the Weak Stage Order Conditions for Runge-Kutta
Methods

Abhijit Biswas  Computer, Electrical, and Mathematical Sciences & Engineering
Division
King Abdullah University of Science and Technology
Thuwal 23955
Saudi Arabia<EMAIL_ADDRESS>https://math.temple.edu/~tug14809 ,
David Ketcheson  Computer, Electrical, and Mathematical Sciences & Engineering
Division
King Abdullah University of Science and Technology
Thuwal 23955
Saudi Arabia<EMAIL_ADDRESS>https://www.davidketcheson.info ,
Benjamin Seibold  Department of Mathematics
Temple University
1805 North Broad Street
Philadelphia, PA 19122<EMAIL_ADDRESS>http://www.math.temple.edu/~seibold
and  David Shirokoff  Corresponding author, Department of Mathematical
Sciences
New Jersey Institute of Technology
University Heights
Newark, NJ 07102<EMAIL_ADDRESS>https://web.njit.edu/~shirokof

###### Abstract.

Runge-Kutta (RK) methods may exhibit order reduction when applied to stiff
problems. For linear problems with time-independent operators, order reduction
can be avoided if the method satisfies certain weak stage order (WSO)
conditions, which are less restrictive than traditional stage order
conditions. This paper outlines the first algebraic theory of WSO, and
establishes general order barriers that relate the WSO of a RK scheme to its
order and number of stages for both fully-implicit and DIRK schemes. It is
shown in several scenarios that the constructed bounds are sharp. The theory
characterizes WSO in terms of orthogonal invariant subspaces and associated
minimal polynomials. The resulting necessary conditions on the structure of RK
methods with WSO are then shown to be of practical use for the construction of
such schemes.

###### Key words and phrases:

Weak stage order, Runge-Kutta, order-reduction, B-convergence, DIRK methods

###### 2000 Mathematics Subject Classification:

65L04; 65L20; 65M12.

##  1\. Introduction

Runge-Kutta (RK) methods are widely used to numerically advance initial value
problems in time. It is well known that the rate of convergence of the
numerical solution may be lower than the design order of the Runge-Kutta
method, particularly when integrating stiff ODEs, stiff PDEs, or initial-
boundary value problems (IBVPs) with time-dependent boundary conditions [9, 4,
14, 40, 42, 28, 25, 30]. In general, this _order reduction_ phenomenon can be
avoided by the use of fully-implicit RK methods with high stage order [33,
18]. Order reduction stemming from time-dependent PDE boundary conditions can
also be avoided by various problem- and method-dependent modifications [13, 1,
2, 11, 3, 12].

In this work we study conditions on the Runge-Kutta coefficients that are
weaker than those of high stage order and guarantee high-order convergence for
a large class of problems where order reduction would typically appear. We are
particularly motivated by the existence of diagonally implicit RK (DIRK)
methods that satisfy these _weak stage order_ (WSO) conditions. The WSO
approach is appealing since it involves only the design of RK schemes, which
can then be applied in standard fashion without any problem-specific
modifications. DIRK methods are particularly appealing due to their
computational efficiency and simplicity. Herein we study the relationship
between the WSO and other salient properties—such as order, number of stages,
and stiff accuracy—of RK methods in general and DIRK methods in particular.
These relationships constrain and guide the design of such methods in a way
that has greatly facilitated their construction (see [22, 6]).

###  1.1. Fundamentals and Notation

We consider the initial value problem

$u^{\prime}(t)=f(t,u(t)),\ u(0)=u_{0};\ u\in\mathbb{R}^{m},\
f:\mathbb{R}\times\mathbb{R}^{m}\to\mathbb{R}^{m}$ (1)

and its numerical approximation $u_{n}\approx u(t_{n})$, where $t_{n}=n\Delta
t$, via an $s$-stage Runge-Kutta (RK) scheme, where the one-step update rule

$\displaystyle g_{i}$ $\displaystyle=u_{n}+\Delta
t\sum_{j=1}^{s}a_{ij}f(t_{n}+c_{j}\Delta t,g_{j})\;,\quad i=1,2,\ldots,s\;$
(2a) $\displaystyle u_{n+1}$ $\displaystyle=u_{n}+\Delta
t\sum_{j=1}^{s}b_{j}f(t_{n}+c_{j}\Delta t,g_{j})\;,$ (2b)

is concisely represented via the Butcher tableau

$\begin{array}[]{c|c}\vec{c}&A\\\
\hline\cr&\vec{b}^{\,T}\end{array}\;,\;\textrm{where}\;A=(a_{ij})_{i,j=1}^{s},\quad\vec{b}=(b_{1},\ldots,b_{s})^{T},\quad\vec{c}=(c_{1},\ldots,c_{s})^{T}.$

Herein we always assume

$\vec{c}=A\vec{e}\;,$ (3)

where $\vec{e}=(1,\ldots,1)^{T}\in\mathbb{R}^{s}$, and denote
$C=\mathrm{diag}(c_{1},c_{2},\ldots,c_{s})$ so that $\vec{c}=C\vec{e}$. We let
$n_{c}$ denote the number of distinct abscissas $c_{j}$ (e.g., if
$c_{1}=\dots=c_{s}$ then $n_{c}=1$; in turn, if all $c_{j}$ are distinct (the
scheme is non-confluent) then $n_{c}=s$).

Schemes for which $A$ is lower-triangular ($a_{ij}=0$ for $j>i$) are called
_diagonally implicit Runge-Kutta_ (DIRK) methods. If furthermore all diagonal
entries are equal, i.e., $a_{ii}=\gamma$ for $i=1,\dots,s$, then the scheme is
called _singly diagonally implicit_ (SDIRK). We refer to DIRK methods with
$a_{11}=0$ as EDIRKs. Implicit schemes that are not diagonally implicit are
referred to as _fully implicit_.

When a RK scheme is applied to the scalar ODE $u^{\prime}(t)=\lambda u$, it
results in the iteration $u_{n+1}=R(z)u_{n}$, where $z:=\lambda\Delta t$ and

$R(z):=1+z\vec{b}^{\,T}(I-zA)^{-1}\vec{e}=\frac{\det(I-zA+z\vec{e}\vec{b}^{\,T})}{\det(I-zA)}$
(4)

is the _stability function_ of the scheme.

A RK scheme has _(classical) order $p$_ if for sufficiently smooth solutions,
the error incurred in a single step (2) is ${\mathcal{O}}(\Delta t^{p+1})$ as
$\Delta t\rightarrow 0$. This imposes a set of _order conditions_ [17] on the
Butcher coefficients. For $p\leq 3$, the order conditions are

$\vec{b}^{\,T}A^{j}\,\vec{c}^{\,k}=\frac{1}{(j+k+1)\dots(k+1)}\quad\text{for}\quad
0\leq j+k\leq p-1\;\mbox{and}\;j\/,\,k\geq 0\;,$ (5)

where $\vec{c}^{\,k}:=\left(c_{1}^{k},c_{2}^{k},\ldots,c_{s}^{k}\right)^{T}$
denotes component-wise exponentiation. For $p>3$, additional order conditions
beyond (5) are required, e.g., $\vec{b}^{\,T}CA\vec{c}=\frac{1}{8}$ for $p=4$.
Related to the order conditions are the widely used _simplifying assumptions_
:

$\displaystyle B(\xi):$ $\displaystyle\vec{b}^{\,T}\vec{c}^{\,k-1}$
$\displaystyle=\tfrac{1}{k}$ for $\displaystyle k=1,2,\ldots,\xi\;,$ (6)
$\displaystyle C(\xi):$ $\displaystyle\vec{\tau}^{(k)}$ $\displaystyle=0$ for
$\displaystyle k=1,2,\ldots,\xi\;.$ (7)

Here we use (slightly deviating from standard notation) the $k$-th _stage
order residual_

$\displaystyle\vec{\tau}^{(k)}$
$\displaystyle:=A\vec{c}^{\,k-1}-\tfrac{1}{k}\vec{c}^{\,k}=\left(AC^{k-1}-\tfrac{1}{k}C^{k}\right)\vec{e},\quad\textrm{for
}k\geq 1\;.$ (8)

Note that $B(\xi)$ and $C(\xi)$ determine the order of accuracy of the
quadrature and subquadrature rules on which the RK method is based [10].

The _stage order_ of a RK scheme is $\tilde{q}=\min\\{q_{1},q_{2}\\}$, where
$q_{1}$, $q_{2}$ are the largest integers such that $B(q_{1})$ and $C(q_{2})$
hold. Stage order $\tilde{q}$ guarantees that each RK stage is accurate to
order $\tilde{q}$, implying that the method itself has order at least
$\tilde{q}$. Notice that (3) implies $\tau^{(1)}=0$, hence $\tilde{q}\geq 1$.
It is worth noting that even when $\vec{\tau}^{(k)}\neq 0$, the $p$-th order
conditions (5) imply the orthogonality relations

$\vec{b}^{\,T}A^{j}\vec{\tau}^{(k)}=0,\quad\textrm{for}\;1\leq j+k\leq
p-1,\;\textrm{with}\;j\geq 0\;\textrm{and}\;k\geq 1\;,$ (9)

and in particular, $\vec{b}^{\,T}\vec{\tau}^{(k)}=0$ for $1\leq k\leq p-1$.

Lastly, a RK scheme is called _stiffly accurate_ [33] if the last row of $A$
equals $\vec{b}^{\,T}$, i.e., $a_{sj}=b_{j}$ for $j=1,\dots,s$. When $A$ is
invertible, stiff accuracy implies $\vec{b}^{\,T}A^{-1}\vec{e}=1$, and hence
$\lim_{z\rightarrow\infty}R(z)=0$.

###  1.2. Weak Stage Order, Stiff Limits, and Error Formulas

Analogous to (6) and (7), one can introduce the conditions on the RK scheme
$(A,\vec{b})$:

$\displaystyle S(\xi):$
$\displaystyle\vec{b}^{\,T}A^{j}\vec{\tau}^{(k)}=0,\quad\textrm{for
all}\;0\leq j\leq s-1,\;1\leq k\leq\xi\;.$ (10)

Based on these conditions (10), we define

######  Definition 1.1.

(Weak stage order, WSO) The weak stage order $q$ of an $s$-stage RK scheme
$(A,\vec{b})$ is the largest integer for which condition $S(q)$ in (10) holds.
If $S(\xi)$ holds for every $\xi\geq 1$, then $q=\infty$.

As the name suggests, the conditions for WSO $q$ relax those for stage order
$\tilde{q}$, i.e., $\tilde{q}\leq q$ always. The relevance of WSO in the
context of order reduction can be understood by examining the error for the
(stiff) Prothero-Robinson problem [33, 18],

$u^{\prime}=\lambda\left(u-\phi(t)\right)+\phi^{\prime}(t);\quad
u(0)=\phi(0),\ \text{with}\ \text{Re}(\lambda)\leq 0\;,$ (11)

where the given smooth function $\phi(t)$ is a global attractor of the
solution $u(t)=\phi(t)$ when $\textrm{Re}(\lambda)<0$. Problem (11) can be
made arbitrarily stiff by choosing $|\lambda|$ large. As in [18, (15.9) Chap.
IV.15], one can derive exact formulas for the RK scheme’s _global error_ (with
equal time-steps):

$u_{n+1}-\phi(t_{n+1})=(R(z))^{n+1}(u_{0}-\phi(0))+\sum_{j=0}^{n}(R(z))^{n-j}\vec{\delta}_{\Delta
t}(t_{j})\;,$

where $z=\lambda\Delta t$ and $R(z)$ are defined as in (4), and
$\vec{\delta}_{\Delta t}(t_{n})$ is the _local error_. For a $p$-th order
method, $\vec{\delta}_{\Delta t}(t_{n})$ has the form

$\displaystyle\vec{\delta}_{\Delta t}(t_{n})$
$\displaystyle=-z\vec{b}^{\,T}(I-zA)^{-1}\vec{\mathcal{E}}(t_{n})+{\mathcal{O}}(\Delta
t^{p+1})\;,$ (12) $\displaystyle\vec{\mathcal{E}}(t_{n})$
$\displaystyle=\sum_{k\geq 1}\frac{(\Delta
t)^{k}}{(k-1)!}\vec{\tau}^{(k)}\phi^{(k)}(t_{n})\;,$ (13)

where $\phi^{(k)}(t_{n})$ is the $k$-th derivative of $\phi$ at $t_{n}$.

In the classical RK convergence theory one examines the non-stiff limit where
both the time scale $1/|\lambda|$ and the solution $\phi(t)$ are fully
resolved. Considering $\Delta t\to 0$ with $z=\mathcal{O}(\Delta t)$ allows to
(i) Taylor-expand (12) in $z$ which leads to terms like
$\vec{b}^{\,T}A^{\ell}\vec{\tau}^{(k)}$, and (ii) use the relations (9) to
show that the local error $\vec{\delta}_{\Delta t}(t_{n})={\mathcal{O}}(\Delta
t^{p+1})$.

In a stiff problem setting, one fully resolves the slow time scale defined by
the solution $\phi(t)$ (i.e., $\Delta t|\phi^{\prime}(t)|\ll 1$) but under-
resolves the fast time scale $1/|\lambda|$ (i.e., $|z|=\Delta
t|\lambda|>{\mathcal{O}}(1)$). One can then investigate the local error
$\vec{\delta}_{\Delta t}(t_{n})$ under the _semi-stiff limit_ given by $\Delta
t\rightarrow 0$ while simultaneously holding $z=\textrm{const.}$ with
$\textrm{Re}(z)<0$; or even the stiff limit in the sense of Prothero-Robinson,
given by $\Delta t\rightarrow 0,\lambda\rightarrow-\infty$ such that
$z\rightarrow-\infty$ (assuming $I-zA$ is invertible).

In the semi-stiff limit (for each fixed $z$), the local error
$\vec{\delta}_{\Delta t}(t_{n})$ is at least order $\tilde{q}+1$ because
$\vec{\tau}^{(k)}=0$ for $k\leq\tilde{q}$, but in general it will fail to be
of order $p+1$. The idea behind WSO is to ensure that the functions

$\displaystyle\widetilde{W}_{k}(z):=z\vec{b}^{\,T}(I-zA)^{-1}\vec{\tau}^{(k)}\equiv
0\textrm{~{}~{}for~{}~{}}k=1,\ldots,q$

vanish, yielding a local error $\vec{\delta}_{\Delta
t}(t_{n})={\mathcal{O}}(\Delta t^{q+1})$ in the semi-stiff limit. WSO is the
least restrictive condition to guarantee $\widetilde{W}_{k}(z)=0$ for $k\leq
q$ [22, Thm. 2].

###  1.3. Prior Work on Avoiding Order Reduction

The prior analysis highlights that in the semi-stiff limit, (i) the order of
the local error is reduced; (ii) methods with high stage order $\tilde{q}$
avoid this order reduction; and (iii) methods with high WSO $q$ avoid order
reduction, at least for a class of problems including (11).

Ostermann & Roche [28] showed that, for linear problems, methods satisfying

$W_{k}(z):=\frac{k\;\vec{b}^{\,T}(I-zA)^{-1}\vec{\tau}^{(k)}}{R(z)-1}\equiv
0\quad\textrm{for}\;k=1,2,\ldots,q\;,$ (14)

along with other technical assumptions, avoid order reduction111i.e., attain a
convergence rate of $\min\\{p,q+2+\bar{\nu}\\}$ where $\bar{\nu}$ depends on
the PDE, see [28, Thm. 2]. when applied to linear parabolic PDEs; see [28] for
RK methods and [29] for the related Rosenbrock W- methods (ROW). Hundsdorfer
showed that a generalization of (14) implies a convergence rate of $O(\Delta
t^{q+1})$ for linear-implicit RK schemes applied to linear ODEs [19, Thm.
3.3]222See also [8, Thm. 3.3] which used $W_{k}(z)$ to achieve improved
convergence estimates.. Montijano showed that a necessary condition for
convergence at order $q$ on the problem (11) is $b^{T}A^{-1}c^{j}=1$ for
$0\leq j\leq q$ [26, 15].

Because (14) contains the variable $z$, it is somewhat impractical as a means
to design RK schemes. Scholz [41] converted (14) into direct “order
conditions” on the scheme coefficients. By expanding terms like
$z\vec{b}^{\,T}(I-zA)\vec{\tau}^{(k)}$ in the variable $H=z/(1-z\gamma)$
(which is possible only for SDIRK methods), Scholz developed B-convergence
order barriers (up to order $4$) for ROW methods and demonstrated that the
weaker conditions are compatible with ROW (SDIRK) structures by constructing
schemes with up to $p=q=4$ in (14). The complexity of Scholz’s conditions for
ROW methods were subsequently simplified by Ostermann & Roche [29, Eq.
(3.12’)] to (effectively) the conditions (10), where the $A$ and
$\vec{\tau}^{(k)}$ in (10) are ROW analogues of the corresponding RK scheme
quantities. For further discussion of order reduction in the context of stiff
problems we refer to [7, Chapter 4], [15, Chapter 7] and [18, Chapter IV.15]
and references therein.

Although originally derived only for SDIRK methods, the “order conditions”
(10) are in fact necessary conditions for any RK scheme to avoid order
reduction in (11); see §1.2 and [34, Eqs. (20)–(21)], [35, Thm. 3], [38, 22].
A geometric interpretation of (10) is provided in [38]. Recently, [36, Chap.
6] and [37] extended (10) to generalized-structure additively partitioned
Runge-Kutta (GARK) methods.

###  1.4. Contribution of this Work

While much is known about RK order conditions, less is known about the
solvability of the WSO conditions, and more generally, the simultaneous
solution of both sets of order conditions. The DIRK structure further limits
the solvability of the WSO conditions relative to fully implicit methods.

Early attempts at constructing DIRK schemes with WSO [38, 22] relied upon the
simplifying assumption that the vectors $\tau^{(k)}$ in (10) be eigenvectors
of $A$. For instance, DIRK schemes up to order $4$ and WSO $3$ have been
provided in [22]. This simplifying “eigenvector criterion” limits the
construction of WSO schemes for DIRKs (with invertible $A$) to $q\leq 3$
(shown below). A key motivating factor for this work is the need to construct
DIRK schemes with high WSO $q>3$; in order to enable this, in the present work
the vectors $\vec{\tau}^{(k)}$ are not assumed to be eigenvectors of $A$.
Moreover, many fundamental results in this work are not restricted to DIRK
schemes.

This paper provides the first known bounds relating a RK scheme’s WSO $q$ to
the order $p$ and number of stages $s$. The starting point (§2) is to recast
WSO in terms of orthogonal invariant subspaces related to the matrix $A$. From
there we establish general bounds on $s$, $q$, $p$ for (arbitrary) RK schemes
(§3). In §4 we obtain a new formula for the stability function in terms of a
family of polynomials (orthogonal with respect to a linear functional). We
then prove stricter bounds on $s$, $p$, $q$ for DIRK schemes (§5). The study
also reveals necessary conditions for WSO $q>3$ (§5.2), which provides the
theoretical foundation for devising new schemes in the companion paper [6]. We
conclude with examples (§6) demonstrating the sharpness of the bounds.

##  2\. Weak Stage Order, Invariant Subspaces, and Minimal Polynomials

In this section we recast weak stage order (WSO) in terms of orthogonal left
and right $A$-invariant subspaces and associated minimal polynomials. This
provides a framework to obtain the main results (e.g., bounds on the order $p$
and WSO $q$ in terms of the number of stages $s$ for both DIRKs and fully
implicit schemes).

The algebraic definition (Def. 1.1) of WSO has a corresponding geometric
interpretation. The combination of the Krylov subspaces generated by the first
$m$ stage order residuals
$\\{\vec{\tau}^{(1)},\vec{\tau}^{(2)},\ldots,\vec{\tau}^{(m)}\\}$ defines the
space

$K_{m}:=\operatorname{span}\left\\{\vec{\tau}^{(1)},A\vec{\tau}^{(1)},\ldots,A^{s-1}\vec{\tau}^{(1)},\vec{\tau}^{(2)},A\vec{\tau}^{(2)},\ldots,A^{s-1}\vec{\tau}^{({m})}\right\\}.$
(15)

If $m=\infty$, we define $K_{\infty}$ as the smallest $A$-invariant space
containing $\vec{\tau}^{(k)}$ for all $k\geq 1$. Note that due to the Cayley-
Hamilton theorem, $K_{m}$ is an $A$-invariant subspace (i.e., $A\vec{v}\in
K_{m}$ for any $\vec{v}\in K_{m}$). A RK scheme then satisfies condition
$S(q)$, given by (10), if $\vec{b}\perp K_{q}$. This leads to an alternative
definition of WSO as the largest integer $q$ for which there exists an
$A$-invariant vector space $V$ with $\vec{\tau}^{(k)}\in V,\ \text{for}\ 1\leq
k\leq q$ and $\vec{b}^{\,T}y=0\ \text{for all}\ y\in V$.

In addition to the right Krylov subspace $K_{m}$ defined in (15), we introduce
the left Krylov subspace generated by $A$ and $\vec{b}$:

$Y:=\operatorname{span}\left\\{\vec{b},A^{T}\vec{b},\ldots,(A^{T})^{s-1}\vec{b}\right\\}.$
(16)

The Cayley-Hamilton theorem implies that $Y$ is left $A$-invariant, i.e.,
$A^{T}\vec{v}\in Y$ for all $\vec{v}\in Y$.

######  Lemma 2.1.

Given a Runge-Kutta scheme defined by $A,\vec{b}$, let $q\geq 1$ be an integer
or $\infty$, and $K_{q}$ and $Y$ be the subspaces (15) and (16) of
$\mathbb{R}^{s}$. Then the scheme has weak stage order of at least $q$ if and
only if $Y$ and $K_{q}$ are orthogonal, in which case

$\operatorname{dim}(Y)+\operatorname{dim}(K_{q})\leq s\;.$ (17)

This lemma is proved by observing that the WSO conditions (10) imply that
$\vec{v}^{T}\vec{w}=0$ for every $\vec{v}\in Y$ and $\vec{w}\in K_{q}$.

RK methods with (classical) stage order $q$ have $\dim(K_{q})=0$ and
$\dim(K_{q+1})>0$. Except for schemes that were specifically designed with
high WSO, almost all irreducible schemes in the literature have $\dim(Y)=s$,
so Lemma 2.1 implies that their WSO is no greater than their stage order.
However, some existing schemes have $\dim(Y)<s$; for example, explicit methods
with one or more entries on the first subdiagonal equal to zero. These include
the 8th-order method of Prince & Dormand [32], as well as some of the fifth-
order SSP methods of Ruuth & Spiteri [39].

Our results in this section will be based on the linear algebra of orthogonal
left and right invariant subspaces applied to $Y$ and $K_{m}$. We therefore
first review and establish key results (e.g., eigenspace, eigenvalue and
minimal polynomials) pertaining to the action of a matrix on its orthogonal
invariant subspaces. We employ the standard terminology: a _monic polynomial_
$p(x)$ is defined as having a leading (i.e., highest degree term) coefficient
of one; the _characteristic polynomial_ of a matrix is defined as the (monic)
polynomial $\textrm{char}_{A}(x):=\det(xI-A)$; the _minimal polynomial_ $p(x)$
of a matrix $A$ is the monic polynomial of smallest degree for which $p(A)=0$.
While it is often the case that (for instance when the eigenvalues of $A$ are
distinct) the minimal polynomial is the characteristic polynomial, i.e.,
$p(x)=\textrm{char}_{A}(x)$, in general $p(x)$ is of lower degree than
$\textrm{char}_{A}(x)$ when $A$ has repeated eigenvalues.

For a real matrix $A$, [5, Chap. 8 & 9A] the minimal polynomial $p(x)$ of $A$:
(i) is unique; (ii) has real coefficients333More precisely, when $A$ is real
with $p_{r}(x)$ the minimal polynomial with coefficients in $\mathbb{R}$, and
$p_{c}(x)$ is the minimal polynomial with coefficients over $\mathbb{C}$, then
$p_{r}(x)=p_{c}(x)$.; (iii) divides the characteristic polynomial of A (thus
every root of $p(x)$ is an eigenvalue of $A$ and the
$\textrm{deg}\;p\leq\textrm{deg}\;\textrm{char}_{A}(x)$); moreover, (iv) every
eigenvalue of $A$ (not including multiplicity) is a root of $p(x)$; (v) every
polynomial $\tilde{p}(x)$ satisfying $\tilde{p}(A)=0$ is divisible by $p(x)$.
The following two theorems concern minimal polynomials for matrices restricted
to invariant subspaces and will be employed in our study of WSO.

######  Theorem 2.2.

($A$-invariant subspaces) Let $A\in\mathbb{R}^{s\times s}$, and
$U\subseteq\mathbb{R}^{s}$ an $A$-invariant subspace with dimension
$d:=\operatorname{dim}(U)$ ($d=0$ is possible). Then there exists a unique
monic polynomial $p(x)$ of minimal degree and coefficients in $\mathbb{R}$,
such that:

$\displaystyle p(A)\vec{u}=0,\quad\forall\vec{u}\in U\;.$ (18)

This polynomial has the following properties:

  1. (a)

$p(x)$ divides the characteristic polynomial $\textrm{char}_{A}(x)$;

  2. (b)

every root of $p(x)$ is an eigenvalue of $A$; and

  3. (c)

$\deg p\leq\operatorname{dim}(U)$.

If, in addition, $U$ has the form
$U=\textrm{span}\\{\vec{v},A\vec{v},\ldots,A^{d-1}\vec{v}\\}$ (or $U=\\{0\\}$
when $\operatorname{dim}(U)=0$) where $A^{j}\vec{v}$ are linearly independent
for $j=0,\ldots,d-1$, then

  1. (d)

$\deg p=\operatorname{dim}(U)$.

  2. (e)

Condition (18) is equivalent to $p(A)\vec{v}=0$.

We include a proof of Theorem 2.2 in Appendix A, using straightforward
generalizations of textbook arguments: $p(x)$ is the minimal polynomial for
$A$ restricted to the subspace $U$ (from which properties (a–e) follow).
Theorem 2.2 can then be extended to the case when $A$ has both a left and
right orthogonal invariant subspace.

######  Theorem 2.3.

(Left and right orthogonal $A$-invariant subspaces) Let
$A\in\mathbb{R}^{s\times s}$ be a real matrix and $U\subseteq\mathbb{R}^{s}$
and $V\subseteq\mathbb{R}^{s}$ orthogonal subspaces, where $U$ is
$A$-invariant and $V$ is $A^{T}$-invariant. Denote the minimal polynomials in
(18) from Theorem 2.2 as $p(x)$ for $U$ and $q(x)$ for $V$, i.e.,
$p(A)\vec{u}=0\;\forall\vec{u}\in U$ and $q(A^{T})\vec{v}=0\;\forall\vec{v}\in
V$. Then the product $p(x)q(x)$ divides the characteristic polynomial
$\textrm{char}_{A}(x)$.

Again, a proof of Theorem 2.3 is included in Appendix A for completeness.
These results have a direct consequence for RK schemes with high WSO.

######  Example 2.4.

To highlight Theorem 2.2 and Theorem 2.3, consider

$A=\begin{pmatrix}\phantom{-}1&\phantom{-}0&0&0\\\
\phantom{-}1&\phantom{-}1&0&0\\\ \phantom{-}2&-1&1&0\\\
-2&\phantom{-}0&1&2\end{pmatrix},\quad U=\begin{pmatrix}1&0&\phantom{-}1\\\
0&1&-1\\\ 0&1&\phantom{-}1\\\
1&0&-1\end{pmatrix},\quad\vec{v}=\begin{pmatrix}\phantom{-}1\\\
\phantom{-}1\\\ -1\\\ -1\end{pmatrix}.$

The matrix $A$ has an $A$-invariant subspace spanned by the columns of $U$,
and an orthogonal $A^{T}$-invariant space spanned by $\vec{v}$. The associated
minimal polynomials are $p(x)=(x-1)^{3}$, s.t. $p(A)U=0$, and $q(x)=(x-2)$,
s.t. $\vec{v}^{T}q(A)=0$. The product
$p(x)q(x)=(x-1)^{3}(x-2)=\textrm{char}_{A}(x)$ equals, and thus divides, the
characteristic polynomial. The example here will also be used to illustrate
the lemmas in §5.2.

######  Corollary 2.5.

Let $(A,\vec{b})$ be a Runge-Kutta method. Then there exists a unique (non-
zero, real coefficient) monic polynomial $Q(x)$ of minimal degree such that

$\vec{b}^{\,T}Q(A)=0\;.$ (19)

The polynomial $Q(x)$ has $\deg Q=\dim Y$ and divides $\textrm{char}_{A}(x)$.

If, in addition, the RK method has WSO $q\geq 2$, then there exists a unique
(non-zero, real coefficient) monic polynomial $P(x)$ of minimal degree, such
that

$P(A)\vec{\tau}^{(k)}=0,\quad\mathrm{for}\;k=2,\ldots,q\;.$ (20)

Moreover, $\deg P\leq\dim K_{q}$, and the product $P(x)Q(x)$ divides
$\textrm{char}_{A}(x)$.

In (19–20), $K_{q}$ and $Y$ are the subspaces defined by (15) and (16); and
$\vec{\tau}^{(k)}$ are the stage order residuals. We refer to $Q(x)$ and
$P(x)$ in Corollary 2.5 as the _minimal polynomials_ for $Y$ and $K_{q}$,
respectively.

###### Proof.

For the existence of $Q$ and its properties, set $U=Y$ and apply Theorem
2.2(d–e). For the existence of $P$ and its properties, set $U=K_{q}$ and $V=Y$
and apply Theorem 2.3, where we denote the polynomials $p(x)$, $q(x)$ in
Theorem 2.3 as $P(x)$, $Q(x)$ respectively. The only point to prove is that
the condition

$P(A)\vec{v}=0,\;\forall\vec{v}\in K_{q}\;,$ (21)

from Theorem 2.3 is equivalent to (20) where $\vec{v}$ is restricted to be the
stage order residuals $\vec{\tau}^{(k)}$ for $k=1,\ldots,q$. Clearly, (21)
implies (20) since $K_{q}$ trivially includes $\vec{\tau}^{(k)}$ for
$k=1,\ldots,q$. Conversely, (20) implies that every basis vector in $K_{q}$,
i.e., $A^{j}\vec{\tau}^{(k)}$, satisfies $P(A)A^{j}\vec{\tau}^{(k)}=0$. Hence,
(20) may be used in lieu of (21) to define the minimal polynomial $P(x)$. ∎

######  Remark 2.6.

(Eigenvalues of $A$) A consequence of Corollary 2.5 is that the eigenvalues of
$A$ (including multiplicity) for a RK scheme can be partitioned into three
sets: the roots of $P(x)$, the roots of $Q(x)$, and the roots of $N(x)$, where
$\textrm{char}_{A}(x)=P(x)Q(x)N(x)$. $\spadesuit$

Broadly speaking, the polynomials $Q(x)$ and $P(x)$ are algebraic objects
associated with the corresponding geometric (orthogonal) spaces $K_{q}$ and
$Y$. In the next sections, we examine the solvability of algebraic equations
related to $Q(x)$ and $P(x)$ as well as the spaces $K_{m}$ and $Y$. This will
lead to necessary conditions for high WSO schemes.

##  3\. Weak Stage Order Barriers for General Runge-Kutta Schemes

In this section we systematically study the subspaces $Y$ and $K_{m}$ for an
arbitrary Runge-Kutta scheme, i.e., the results apply to fully implicit RK
schemes, DIRKs, schemes with singular $A$ (e.g., EDIRKs), etc. Later, in §5,
we then restrict to DIRK schemes, which will lead to a refinement of the
results obtained here.

###  3.1. The Subspace $Y$

Here we show that the degrees of the numerator and the denominator of the
Runge-Kutta stability function are determined by the dimension of the subspace
$Y$. When combined with classical results on rational approximations to
$e^{z}$, we obtain a bound on $p$ in terms of $\dim Y$.

######  Theorem 3.1.

(Rational approximations of $e^{z}$, [18, Thm. IV.3.11,Thm. IV.4.18], and [27,
Thm. 2.1]) Let $R(z)=N_{r}(z)/D_{r}(z)$ be a rational approximation to $e^{z}$
of order $p$. Then $p\leq\deg N_{r}+\deg D_{r}$. If, further, $R(z)$ has real
poles then $p\leq\deg N_{r}+1$.

For an $s$-stage RK method, the stability function is a rational approximation
of $e^{z}$ with numerator and denominator of degree at most $s$. In fact, we
can replace $s$ in this statement with $\dim Y$.

######  Theorem 3.2.

Let a RK method $(A,b)$ be given, with $Y$ defined as in (16). Then the
stability function of the method is $R(z)=N_{r}(z)/D_{r}(z)$, where $\deg
N_{r}\leq\dim Y$ and $\deg D_{r}\leq\dim Y$.

###### Proof.

Define an orthogonal matrix $L=(U\;|\;V)\in\mathbb{R}^{s\times s}$, where the
columns of $U\in\mathbb{R}^{s\times d}$, with $d:=\dim Y$, are a basis of $Y$
and the columns of $V$ are a basis of its orthogonal complement. Then
$\vec{b}^{\,T}V=\vec{0}$ and $\vec{b}^{\,T}(I-zA)^{-1}V=\vec{0}$. So the
stability function (4) is

$R(z)=1+z\vec{b}^{\,T}LL^{T}(I-zA)^{-1}LL^{T}\vec{e}=1+z\vec{b}^{\,T}_{U}MU^{T}\vec{e}\;,$
(22)

where $\vec{b}^{\,T}_{U}=\vec{b}^{\,T}U\in\mathbb{R}^{d}$ and
$M=U^{T}(I-zA)^{-1}U\in\mathbb{R}^{d\times d}$. Using this and Cramer’s rule
it can be shown that the numerator and denominator of $R(z)$ are given by
determinants of $d\times d$ matrices whose entries are linear functions of
$z$. Thus $R(z)$ is a rational function with numerator and denominator of
degree at most $d$. ∎

For stiffly accurate schemes with $A$ invertible, $R(z)$ is further
constrained to satisfy

$\displaystyle\lim_{z\rightarrow\infty}R(z)=\lim_{z\rightarrow\infty}N_{r}(z)/D_{r}(z)=0\;,$

so that $\deg N_{r}\leq\deg D_{r}-1\leq\dim Y-1$. For DIRK schemes, the
stability function $R(z)$ has real poles since the eigenvalues of $A$ are
confined to the diagonal entries—which are real. Combining these observations
with Theorem 3.1 and Theorem 3.2, we have proved:

######  Theorem 3.3.

Let a RK scheme $(A,b)$ be given with $Y$ defined as in (16) and stability
function $R(z)=e^{z}+{\mathcal{O}}(z^{p+1})$. Then

$\left\lfloor\frac{p+1+\sigma}{2}\right\rfloor\leq\dim Y\;.$ (23)

If the scheme is diagonally implicit, (23) we have the stronger bound

$p\leq\dim(Y)+1-\sigma\;.$ (24)

In both (23) and (24), $\sigma=1$ if the scheme is stiffly accurate and $A$ is
invertible, and $\sigma=0$ otherwise.

Here $\lfloor x\rfloor$ is the standard floor function for the real number
$x$. Combined with (17) and Corollary 2.5, Theorem 3.3 implies the following:

######  Corollary 3.4.

For an $s$-stage Runge-Kutta scheme with order $p\geq 1$ and weak stage order
$q$ (with $K_{q}$ and $P(x)$ defined in (15) and (20)), one has

$\displaystyle\deg(P)\leq\dim{(K_{q})}\leq\left\\{\begin{array}[]{cl}s-\left\lfloor\tfrac{p+1+\sigma}{2}\right\rfloor,&\textrm{for
fully implicit schemes},\\\ s-p-1+\sigma,&\textrm{for diagonally implicit
schemes},\end{array}\right.$

where $\sigma=1$ if $A$ is invertible and the method is stiffly accurate, and
$\sigma=0$ otherwise.

###  3.2. The Subspace $K_{m}$ for General RK Schemes

Here we provide a lower bound on $\dim(K_{m})$ in terms of $m$ and the number
of distinct abscissas $n_{c}$. The key idea is that $K_{m}$ contains linear
combinations of the vectors $\vec{\tau}^{(k)}$, and for large enough $m$ the
addition of each new vector increases the dimension of the space $K_{m}$.

Let $w(x)=w_{m-1}x^{m-1}+w_{m-2}x^{m-2}+\ldots+w_{0}$ denote an arbitrary
(real) polynomial and let $W(x):=\int_{0}^{x}w(s)\;{\rm ds}$. The space
$K_{m}$ then contains the linear combination

$\displaystyle\sum_{j=1}^{m}w_{j-1}\vec{\tau}^{(j)}$
$\displaystyle=\left(\sum_{j=1}^{m}Aw_{j-1}C^{j-1}-\frac{1}{j}w_{j-1}C^{j}\right)\vec{e},\quad\textrm{using
\eqref{Eq:SOV_MatrixForm}}\;,$
$\displaystyle=\big{(}Aw(C)-W(C)\big{)}\vec{e}\in K_{m}\;.$ (25)

Setting $w(x)=\alpha(x)P_{\rm C}(x)$ in (25), where $P_{\rm C}(x)$ is the
minimal polynomial for $C$ (so that $w(C)=0$), shows that $W(C)\vec{e}\in
K_{m}$. Since $C$ is diagonal, $P_{\rm C}(x)$ has simple roots (one for each
distinct entry of $\vec{c}$) and so $\deg P_{\rm C}=n_{c}$. Thus, for $m\geq
n_{c}+1$, $K_{m}$ contains the subspace $\mathbf{W}_{m}$, where

$\displaystyle\mathbf{W}_{m}$
$\displaystyle:=\left\\{W(C)\vec{e}\;\big{|}\;\textrm{for
all}\;W(x)=\int_{0}^{x}\alpha(x)P_{\rm C}(x)\;{\rm
dx}\;\;\textrm{and}\;\;\alpha(x)\in\mathbf{V}_{m-1-n_{c}}\right\\},$
$\displaystyle\mathbf{V}_{d}$
$\displaystyle=\left\\{\alpha_{0}+\ldots+\alpha_{d}x^{d}\;|\;\alpha_{j}\in\mathbb{R},\textrm{
for }j=0,\ldots,d\right\\}.$

The next theorem follows directly from properties of the space
$\mathbf{W}_{m}$.

######  Theorem 3.5.

(Properties of $K_{m}$) Consider a Runge-Kutta scheme with $n_{c}$ distinct
abscissas and the space $K_{m}$ defined in (15). Then:

  1. (a)

If $m\leq 2n_{c}-1$, then $\dim K_{m}\geq\max\\{m-n_{c},0\\}$.

  2. (b)

If $m=2n_{c}$ and all $c_{j}\neq 0$ (i.e., $P_{\rm C}(0)\neq 0$), then
$\mathbf{W}_{m}=\mathrm{span}\\{\vec{e},\vec{c},\ldots,\vec{c}^{\,n_{c}-1}\\}$
with $\mathbf{W}_{m}\subseteq K_{m}$; in particular $\vec{e}\in K_{m}$ and
$\dim K_{m}\geq n_{c}$.

  3. (c)

If $m=2n_{c}-1$ and $c_{j}=0$ for at least one value of $j$ (i.e., $P_{\rm
C}(0)=0$), then
$\mathbf{W}_{m}=\mathrm{span}\\{\vec{c},\ldots,\vec{c}^{\,n_{c}-1}\\}\subseteq
K_{m}$, in particular $\vec{c}\in K_{m}$ and $\dim K_{m}\geq n_{c}-1$;

  4. (d)

When case (b) or (c) holds, then $K_{m}$ is equal to the smallest
$A$-invariant space containing $\mathbf{W}_{m}$ and $K_{m}=K_{m^{\prime}}$ for
any $m^{\prime}\geq m$ (i.e., increasing $m$ further does not enlarge the
space $K_{m}$).

###### Proof.

For (a), we compute $\dim\mathbf{W}_{m}$ which bounds $\dim K_{m}$ from below.
First note that vectors of the form $w(C)\vec{e}$ (for any polynomial $w(x)$)
are in one-to-one correspondence with their associated polynomial
remainders444Here $w(x)\mod P_{\rm C}(x)$ is the remainder of $w(x)$ when
divided by $P_{\rm C}(x)$; the set of all such remainders can be identified
with the vector space $\mathbf{V}_{n_{c}-1}$. $w(x)\mod P_{\rm C}(x)$ (i.e.,
as vector spaces they are isomorphic555For any polynomial $w(x)\mod P_{\rm
C}(x)$, the identification of $w(x)\mapsto w(C)\vec{e}$ is surjective and
preserves the operations of the vector space. It is also injective: since
$w(C)$ is a diagonal matrix, the vector $w(C)\vec{e}=0$ if and only if
$w(C)=0$, and $w(C)=0$ if and only if $w(x)$ is divisible by $P_{\rm C}(x)$.
Hence, the vectors $w_{1}(C)\vec{e}=w_{2}(C)\vec{e}$ if and only if
$w_{1}(x)-w_{2}(x)=\alpha(x)P_{\rm C}(x)$ for some $\alpha(x)$, i.e.,
$w_{1}(x)\equiv w_{2}(x)\mod P_{\rm C}(x)$.). Thus, $\mathbf{W}_{m}$ has the
same dimension as the range of
$I:\mathbf{V}_{m-1-n_{c}}\rightarrow\mathbf{V}_{n_{c}-1}$ defined by

$I[\alpha]:=\left(\int_{0}^{x}\alpha(x)P_{\rm C}(x)\;{\rm dx}\right)\mod
P_{\rm C}(x)\;,$

where $I[\cdot]$ is the linear (integral) operator that maps
$\mathbf{V}_{m-1-n_{c}}$ (which has dimension $m-n_{c}$) into polynomials
modulo $P_{\rm C}(x)$, (which has dimension $n_{c}$). We now claim that
$I[\cdot]$ has a trivial null space whenever $m\leq 2n_{c}-1$. Suppose
not—then there are polynomials $\alpha(x)\neq 0$ ($\deg\alpha\leq n_{c}-2$)
and $\beta(x)$ ($\deg\beta\leq n_{c}-1$) such that

$\int_{0}^{x}\alpha(x)P_{\rm C}(x){\rm dx}=\beta(x)P_{\rm C}(x)\;.$

Differentiating both sides yields:

$\left(\alpha(x)-\beta^{\prime}(x)\right)P_{\rm C}(x)=\beta(x)P_{\rm
C}^{\prime}(x)\;.$ (26)

The minimal polynomial $P_{\rm C}(x)$ contains $n_{c}\geq 1$ distinct roots,
and hence $P_{\rm C}^{\prime}(x)$ has no common root with $P_{\rm C}(x)$
(i.e., $P_{\rm C}(x)$ and $P_{\rm C}^{\prime}(x)$ are relatively prime in
$\mathbb{R}[x]$). Thus, (26) requires that $P_{\rm C}^{\prime}(x)$ divide
$\left(\alpha(x)-\beta^{\prime}(x)\right)$ and $P_{\rm C}(x)$ divide
$\beta(x)$, which is impossible in light of the degrees of $\alpha(x)$,
$\beta(x)$. Thus no $\alpha(x)$ exists, and whenever $n_{c}+1\leq m\leq
2n_{c}-1$ one has the following dimension:

$\dim\mathbf{W}_{m}=\dim\mathbf{V}_{m-1-n_{c}}=(m-n_{c})\leq\dim K_{m}\;.$
(27)

Combining (27) with the trivial $\dim K_{m}\geq 0$ when $m<n_{c}+1$ proves
(a).

For (b), when $m=2n_{c}$, (26) has the (unique) family of solutions
$\beta(x)=\gamma P_{\rm C}(x)$ and $\alpha(x)=2\gamma P_{\rm C}^{\prime}(x)$
($\gamma\in\mathbb{R}$). However, $I[\alpha]=\gamma P_{\rm C}(x)^{2}-\gamma
P_{\rm C}(0)^{2}=-\gamma P_{\rm C}(0)^{2}\mod P_{\rm C}(x)$, showing that
$\gamma=0$ is still the only solution to $I[\alpha]=0$. Hence, $I[\cdot]$ has
a trivial null space, $\mathbf{W}_{m}$ contains the vector $\vec{e}$ (e.g.,
take $\gamma=-P_{\rm C}(0)^{-2}$), and $\dim\mathbf{W}_{m}=n_{c}$.

For (c), if $P_{\rm C}(0)=0$, then every $W(x)$ appearing in $\mathbf{W}_{m}$
is divisible by $x^{2}$ so that $W(x)$ and $P_{\rm C}(x)$ share a common
factor of $x$. Thus $x$ divides the remainder $W(x)\mod P_{\rm C}(x)$ as well
as every polynomial in the range of $I[\cdot]$, i.e., for every
$\alpha\in\mathbf{V}_{m-1-n_{c}}$, $I[\alpha]=xr(x)$ for some $r(x)$ with
$\deg r\leq n_{c}-2$. When $m=2n_{c}-1$ the range of $I[\cdot]$ has dimension
$n_{c}-1$ and includes all monomials $x^{j}$ for $1\leq j\leq n_{c}-1$.

The vector $\vec{c}$ has an explicitly construction: since $P_{\rm C}(x)$ and
$P_{\rm C}^{\prime}(x)$ are relatively prime, polynomials $\beta(x)$
($\deg\beta\leq n_{c}-1$) and $\gamma(x)$ ($\deg\gamma(x)\leq n_{c}-2$) exist
such that

$\gamma(x)P_{\rm C}(x)=\beta(x)P_{\rm C}^{\prime}(x)+1\;.$

Set $\alpha(x):=\gamma(x)-\beta^{\prime}(x)$ ($\deg\alpha\leq n_{c}-2$), so
that $\alpha(x)P_{\rm C}(x)=\frac{d}{dx}\left(\beta(x)P_{\rm C}(x)\right)+1$.
Since $P_{\rm C}(0)=0$, integrating both sides shows that $I[\alpha]=x$ and
hence: $\vec{c}\in\mathbf{W}_{m}$.

For (d), let $V$ be the smallest $A$-invariant space containing
$\mathbf{W}_{m}$; note that $K_{m}$ always contains $V$. When (b) or (c)
holds, every $\vec{\tau}^{(k)}$ for arbitrary $k$ is a linear combination of
$A\vec{c}^{\,k}\in V$ and $\vec{c}^{\,k+1}\in V$, hence $K_{m^{\prime}}=V$ for
all $m^{\prime}\geq m$. ∎

When combined with Corollary 3.4, and the fact that the order conditions
$\vec{b}^{\,T}\vec{e}=1$ (for $p=1$) and $\vec{b}^{\,T}\vec{c}=1/2$ (for
$p=2$) are incompatible with $\vec{b}\perp K_{m}$ whenever $\vec{e}\in K_{m}$
or $\vec{c}\in K_{m}$, Theorem 3.5 directly implies the following limitation
of WSO:

######  Theorem 3.6.

(Main Result) A Runge-Kutta scheme with $s$ stages, $n_{c}$ distinct
abscissas, order $p$, and weak stage order $q$ satisfies:

  1. (1)

If the abscissa values are all non-zero ($c_{j}\neq 0$), then $q\leq
2n_{c}-1$.

  2. (2)

If some of the abscissa values are zero ($c_{j}=0$), then either (i) $q\leq
2n_{c}-2$, or (ii) $q=\infty$ in which case $p=1$.

  3. (3)

If $q\leq 2n_{c}-1$, then

$q+\left\lfloor\frac{p+1+\sigma}{2}\right\rfloor\leq s+n_{c}\;,$

where $\sigma=1$ for stiff accuracy with $A$ invertible, and $\sigma=0$
otherwise.

######  Remark 3.7.

For simplicity we have written “order p” in Theorem 3.6 (and also in the later
Theorem 5.4); however, both results hold under the weaker condition that
$R(z)=e^{z}+{\mathcal{O}}(z^{p+1})$ is a $p$-th order approximation as
$z\rightarrow 0$. $\spadesuit$

##  4\. Stability Function in a Basis for $Y$

In general, the Runge-Kutta stability function has a restricted form when $p$
is close to (the optimal) $2s$, and also when $p\geq s$ (e.g., see [18, Chap.
IV.3, p. 47]). In this section, we extend this constrained structure for
$R(z)$ to the more general setting $p\geq\dim Y$ and to when $p$ is close to
$2\dim Y$. The key idea is to write $Q(x)$, and $Y$, in a basis of orthogonal
polynomials related to the (_tall-tree_) $p$-th order conditions

$\vec{b}^{\,T}A^{j}\vec{e}=\frac{1}{(j+1)!},\quad 0\leq j\leq p-1\;.$ (28)

Specifically, define the linear functional $\mathcal{L}(\cdot)$ on the space
of monomials $x^{n}$ with _moments_ $(\mu_{n})_{n\geq 0}$ as

$\mathcal{L}(x^{n}):=\mu_{n},\;\textrm{where}\;\mu_{n}:=\frac{1}{(n+1)!},\;\textrm{for}\;n\geq
0\;.$ (29)

Associated to $\mathcal{L}$ are the Hankel moment matrices
$\mathbf{H}_{n}=(\mu_{i+j-2})_{i,j=1}^{n}$. Appendix B shows that
$\det\mathbf{H}_{n}\neq 0$ for all $n\geq 0$, which is a sufficient condition
to construct a basis of polynomials $(Q_{j})_{j\geq 0}$ satisfying the
orthogonality relation

$\mathcal{L}\left(Q_{i}Q_{j}\right)=\zeta_{i}\delta_{ij},\;\textrm{where}\;\zeta_{0}=1,\;\textrm{and}\;\zeta_{i}=\frac{\det(\boldsymbol{H}_{i+1})}{\det(\boldsymbol{H}_{i})},\;\textrm{if}\;i\geq
1\;.$

Here $\delta_{ij}$ is the Kronecker delta. For the moments in (29), the first
two orthogonal polynomials are $Q_{0}(x)=1$ and $Q_{1}(x)=x-1/2$; subsequent
polynomials $Q_{j}(x)$ satisfy the three term recursion relation given by (see
Appendix B for details):

$Q_{n+1}(x)=xQ_{n}(x)+\xi_{n}^{2}Q_{n-1}(x),\quad\xi_{n}^{2}=\frac{1}{4(4n^{2}-1)},\;\textrm{for}\;n\geq
1\;.$ (30)

For instance, the next two polynomials are

$Q_{2}(x)=x^{2}-\tfrac{1}{2}x+\tfrac{1}{12}\quad\textrm{and}\quad
Q_{3}(x)=x^{3}-\tfrac{1}{2}x^{2}+\tfrac{1}{10}x-\tfrac{1}{120}\;.$

Any degree $d$ polynomial can then be written in terms of the basis
$\\{Q_{1},\dots,Q_{d}\\}$. We use $\alpha_{j}$ to denote the coefficients of
the minimal polynomial for $Y$ in this basis:

$Q(x)=Q_{d}(x)+\alpha_{d-1}Q_{d-1}(x)+\ldots+\alpha_{0}Q_{0}(x),\quad\textrm{where}\;d:=\dim
Y\;.$ (31)

The order conditions then constrain (31) as follows.

######  Lemma 4.1.

($Q$ is “orthogonal” to $Q_{j}$) Consider a RK method with coefficients
$(A,\vec{b})$, and let $Q$ be the minimal polynomial for the subspace $Y$, as
defined in (16), written in the form (31). If the method’s stability function
is $R(z)=e^{z}+{\mathcal{O}}(z^{p+1})$ as $z\rightarrow 0$, then
$\alpha_{j}=0$ for $j\leq p-d-1$.

###### Proof.

Set $N:=p-d-1\geq 0$ (if $N<0$ there is nothing to prove). The $p$-th order
conditions (28) imply
$\vec{b}^{\,T}Q(A)A^{\ell}\vec{e}=\mathcal{L}(x^{\ell}Q(x))$ for all
$0\leq\ell\leq N$. On the other hand, $\vec{b}^{\,T}Q(A)=0$. Hence,
$\mathcal{L}(Q_{\ell}Q)=\zeta_{\ell}\alpha_{\ell}=0$ for all $0\leq\ell\leq
N$, so $\alpha_{j}=0$ for $j\leq N$ (since $\zeta_{j}\neq 0$). If $N\geq d$,
this would imply $\mathcal{L}(QQ_{d})=\zeta_{d}=0$, which is a contradiction.
Note: $p\leq 2d$ or $\lfloor\frac{p+1}{2}\rfloor\leq\dim(Y)$ when $p\geq 1$
and $d\geq 0$ are integers. ∎

Utilizing the the basis $(Q_{j})_{j\geq 0}$, we now work out an expression for
the stability function $R(z)$. We first expand $\vec{b}^{\,T}(I-zA)^{-1}$ in
the Krylov basis $\vec{b}^{\,T}Q_{0}(A)$,$\ldots$, $\vec{b}^{\,T}Q_{d-1}(A)$
with (to be determined) coefficients $\beta_{j}=\beta_{j}(z)$ functions of
$z$:

$\vec{b}^{\,T}(I-zA)^{-1}=\vec{b}^{\,T}\big{(}\beta_{d-1}Q_{d-1}(A)+\beta_{d-2}Q_{d-2}(A)+\ldots+\beta_{0}Q_{0}(A)\big{)}\;.$
(32)

Right-multiplying (32) by $I-zA$ and using the fact that
$\vec{b}^{\,T}Q(A)=0$, leads to the algebraic equation for the coefficients:

$(1-zx)\big{(}\beta_{d-1}Q_{d-1}(x)+\beta_{d-2}Q_{d-2}(x)+\ldots+\beta_{0}Q_{0}(x)\big{)}\equiv
1\mod Q(x)\;.$ (33)

To solve for
$\vec{\beta}=\left(\beta_{0},\beta_{1},\ldots,\beta_{d-1}\right)^{T}$, write
$\vec{Q}(x)=\left(Q_{0}(x),Q_{1}(x),\ldots,Q_{d-1}(x)\right)^{T}$. Using
relation (30), multiplication $x\vec{Q}$ can be written as a matrix
multiplication $\mathbf{S}\vec{Q}$:

$x\vec{Q}=\\!\begin{pmatrix}Q_{1}+\dfrac{1}{2}Q_{0}\\\
Q_{2}-\xi_{1}^{2}Q_{0}\\\ \vdots\\\
Q_{d}-\xi_{d-1}^{2}Q_{d-2}\end{pmatrix}\\!=\mathbf{S}\vec{Q},\;\textrm{where}\;\mathbf{S}:=\\!\begin{pmatrix}\frac{1}{2}&1&0&\cdots&0\\\
\\!-\xi_{1}^{2}&0&1&&\vdots\\\ 0&-\xi_{2}^{2}&0&&\\\ \vdots&&&\ddots&1\\\
0&\cdots&&-\xi_{d-1}^{2}&0\end{pmatrix}\\!-\vec{e}_{d}\vec{\alpha}^{T}.$

Here $\vec{e}_{j}$ is the $j$-th unit vector and
$\vec{\alpha}=(\alpha_{0},\ldots,\alpha_{d-1})^{T}$ are the coefficients of
$Q(x)$ in (31). Using these notations and relations, (33) becomes:

$\displaystyle(1-zx)\vec{\beta}^{T}\vec{Q}\equiv 1\mod Q(x)\;,$
$\displaystyle\left(\left(I-z\mathbf{S}^{T}\right)\vec{\beta}-\vec{e}_{1}\right)^{\\!T\\!}\vec{Q}\equiv
0\mod Q(x)\;.$

The vector $\vec{\beta}$ then has solution:

$\left(I-z\mathbf{S}^{T}\right)\vec{\beta}=\vec{e}_{1},\quad\textrm{and via
Cramer's
rule:}\quad\beta_{j-1}(z)=\frac{\det\big{(}(I-z\mathbf{S}^{T})_{j}\big{)}}{\det(I-z\mathbf{S})}\;,$
(34)

where $(I-z\mathbf{S}^{T})_{j}$ is the matrix $I-z\mathbf{S}^{T}$ with column
$j$ replaced by $\vec{e}_{1}$. Substituting the expression for
$\vec{b}^{\,T}(I-zA)^{-1}$ into $R(z)$ yields:

$R(z)=1+\sum_{j=0}^{d-1}\left(\vec{b}^{\,T}Q_{j}(A)\vec{e}\right)z\beta_{j}(z)\;.$
(35)

Finally, the orthogonality property Lemma 4.1 implies that
$\vec{b}^{\,T}Q_{j}(A)\vec{e}$ agrees with $\mathcal{L}(Q_{j})=\delta_{0j}$
for all $j\leq p-1$. This (significantly) simplifies the summation in (35):

######  Lemma 4.2.

(Stability function when $p\geq\dim Y$) If $R(z)=e^{z}+{\mathcal{O}}(z^{p+1})$
is a $p$-th order approximation as $z\rightarrow 0$ and $p\geq\dim Y$, then

$R(z)=1+z\beta_{0}(z)\;,$ (36)

where $\beta_{0}(z)$ is given by (34) and depends only on the expansion of
$Q(x)$ from (31), i.e., $R(z)$ is a function of $\alpha_{j}$ (for $j=p-\dim
Y,\ldots,\dim Y-1$).

To the authors’ knowledge, the orthogonal polynomials $(Q_{j})_{j\geq 0}$ have
not previously been introduced in the context of RK schemes. However, they are
connected with existing ideas in the (extensive) literature on RK methods. For
instance, the values $\xi_{j}$ defined in (30) appear in the $W$-transform and
in recurrence relations for the shifted Legendre polynomials [18, Chap. IV.5]
commonly used in the construction of high stage order schemes. Lemma 4.1 is
also similar in spirit to the orthogonality relation Lemma 5.15 in [18, Chap.
IV.5] which expands the characteristic polynomial of $C$ (when $n_{c}=s$) in
the shifted Legendre polynomial basis.

######  Remark 4.3.

(Orthogonal polynomials with respect to a linear functional) In the
(semi)classical theory of orthogonal polynomials, the Hankel (moment) matrices
are positive definite and $\mathcal{L}(\cdot)$ defines an inner product as an
integral with respect to a (non-negative) measure [16, 21]. In contrast, for
the matrices $\mathbf{H}_{d}$ considered here, both $\det\mathbf{H}_{d}\neq 0$
and $\zeta_{i}\neq 0$. However, $\zeta_{i}<0$ for some values of $i$, so
$\mathcal{L}(\cdot)$ does not have an inner product representation. While much
of the (semi)classical theory of orthogonal polynomials still holds [16]
(e.g., the three-term recurrence), some properties do not. For instance, the
basis polynomials $Q_{j}(x)$ have complex roots, whereas in the classical
theory of orthogonal polynomials the roots are real and simple. $\spadesuit$

######  Remark 4.4.

(Moment-generating function of $\mathcal{L}$) The moments $(\mu_{n})_{n\geq
0}$ do not arise by accident; their generator function is
$G(z)=(e^{z}-1)/z=\sum_{n\geq 0}\mu_{n}z^{n}$, whose Taylor expansion agrees
up to a suitable order with $(R(z)-1)/z$. The Hankel determinants, or
coefficients of the recurrence relation (30) appear in the J-fraction
(continued fraction) of the generating function. For $G(z)$, this continued
fraction appears in the RK literature through Padé approximations of $e^{z}$
(see e.g. Exercise 4 in Chap. IV.3 or Theorem 5.18 of [18]). $\spadesuit$

######  Remark 4.5.

(Stiff accuracy) Equation (36) can be extended to incorporate the structure
imposed by stiff accuracy (for invertible $A$) using an alternative polynomial
basis $(\tilde{Q})_{n\geq 0}$. Specifically, let $\tilde{\mu}_{0}=1$,
$\tilde{\mu}_{n+1}:=\mu_{n}$ for $n\geq 0$ and
$\tilde{\mathcal{L}}(x^{n}):=\tilde{\mu}_{n}$. The Hankel matrices with
moments $\tilde{\mu}_{n}$ are shown in (49) to have non-zero determinants and
define a basis orthogonal with respect to $\tilde{\mathcal{L}}$, i.e.,
$\tilde{Q}_{0}(x)=1$, $\tilde{Q}_{1}(x)=x-1$, etc. The stability function
(having one less degree of freedom), for stiffly accurate schemes can then be
obtained in this basis. $\spadesuit$

##  5\. Weak Stage Order for DIRK Schemes

This section provides stronger versions of the results from §3 in the setting
of diagonally implicit schemes. The DIRK structure constrains the solvability
of the WSO equations, tightens Theorem 3.6, and restricts how the spectrum of
$A$ is partitioned into the minimal polynomials $Q(x)$, $P(x)$. These results
are crucial in the construction of high-WSO DIRK schemes.

###  5.1. The Subspace $K_{m}$ for DIRK Schemes

We now turn to obtaining lower bounds on the dimension of $K_{m}$. To expedite
the proofs below, we define a truncation map $\left[\cdot\right]_{j}$ which
acts on matrices as $\left[\cdot\right]_{j}:\mathbb{R}^{s\times
s}\rightarrow\mathbb{R}^{j\times j}$ to isolate the upper $j\times j$ block of
$A$, and acts on vectors as
$\left[\cdot\right]_{j}:\mathbb{R}^{s}\rightarrow\mathbb{R}^{j}$ to isolate
the top $s$ components of $\vec{c}$:

$\left[A\right]_{j}:=\begin{pmatrix}{a_{11}}&{a_{12}}&{\cdots}&{a_{1j}}\\\
{a_{21}}&{\ddots}&{}&{\vdots}\\\ {\vdots}&{}&{}&{}\\\
{a_{j1}}&{a_{j2}}&{\cdots}&{a_{jj}}\end{pmatrix}\ \ \text{and}\ \
\left[\vec{c}\right]_{j}:=\begin{pmatrix}{c_{1}}\\\ {c_{2}}\\\ {\vdots}\\\
{c_{j}}\end{pmatrix}.$

For a lower-triangular matrix $A$, the first $j$ ($\leq s$) components of the
vectors $\vec{\tau}^{(k)}$ and $\vec{c}$ are functions of the upper $j\times
j$ sub-block of the matrix $A\in\mathbb{R}^{s\times s}$ only. Moreover:

$\left[A^{n}\right]_{j}=\left(\left[A\right]_{j}\right)^{n},\quad\textrm{and}\quad\left[\vec{c}^{\,n}\right]_{j}=\left(\left[\vec{c}\right]_{j}\right)^{n}.$
(37)

We will also refer to the vectors $\vec{\tau}^{(k)}$ that arise when $A$ is
replaced by $\left[A\right]_{j}$ in (8):

$\displaystyle\vec{\tau}^{(k)}(\left[A\right]_{j})$
$\displaystyle:=\left[A\right]_{j}\left[\vec{c}^{\,k-1}\right]_{j}-\tfrac{1}{k}\left[\vec{c}^{\,k}\right]_{j}\quad\textrm{for
}k\geq 1\;,$ (38)

which commute with the map in the sense that for every $1\leq j\leq s$, $n\geq
0$, and $k\geq 1$:

$\left[\vec{\tau}^{(k)}(A)\right]_{j}=\vec{\tau}^{(k)}(\left[A\right]_{j})\;,\quad\left[A^{n}\vec{\tau}^{(k)}(A)\right]_{j}=\left[A\right]_{j}^{n}\vec{\tau}^{(k)}(\left[A\right]_{j})\;.$
(39)

We also extend $\left[\cdot\right]$ to vector subspaces
$U\subseteq\mathbb{R}^{s}$ as
$\left[U\right]_{j}:=\textrm{span}\\{\left[\vec{u}\right]_{j}\;|\;\vec{u}\in
U\\}$, which implies $\dim U\geq\dim\left[U\right]_{j}$.

The following definition generalizes the concept of an EDIRK (i.e., DIRK
scheme with $a_{11}=0$), and will be useful in the next two lemmas which
improve Theorem 3.5 to the case when $A$ is lower triangular.

######  Definition 5.1.

(GEDIRK) We call a DIRK scheme a generalized EDIRK (short: GEDIRK), if
$\vec{c}$ contains at least one zero and $a_{\ell\ell}=0$, where
$\ell=\min\\{j\>|\;c_{j}=0\\}$ is the index of the first zero in $\vec{c}$.

######  Lemma 5.2.

For a DIRK scheme that is not a GEDIRK scheme, if $m\geq 2n_{c}$ then $\dim
K_{m}\geq n_{c}$ (where the space $K_{m}$ is defined in (15)).

###### Proof.

If $\vec{c}$ has no zero, then Theorem 3.5(b)(d) applies and we are done. Let
$\ell$ be the first occurrence of $c_{\ell}=0$ in $\vec{c}$ (i.e., $c_{j}\neq
0$ for $j<\ell$). Note that $\ell>1$ since $c_{1}=a_{11}\neq 0$ (because $A$
is not a GEDIRK). Theorem 3.5(c)(d) shows that $\mathbf{W}_{m}\subseteq K_{m}$
has dimension $n_{c}-1$, that the vector $\vec{c}\in K_{m}$, and that
$\mathbf{W}_{m}$ is orthogonal to $\vec{e}_{\ell}$, i.e.,
$\vec{e}^{\,T}_{\ell}\vec{u}=0$ for all $\vec{u}\in\mathbf{W}_{m}$. In
addition, $\mathbf{W}_{m}$ contains a vector of the form (since $\ell$ is the
first occurrence of $0$)

$\vec{y}:=(1,\ldots,1,0,\star)\in\mathbf{W}_{m}\;,$

where the first $\ell-1$ entries of $\vec{y}$ are $1$ and $\star$ is
unimportant. Since $K_{m}$ is $A$-invariant, we have $\vec{c}-A\vec{y}\in
K_{m}$; however $\vec{e}^{\,T}_{\ell}(\vec{c}-A\vec{y})=a_{\ell\ell}\neq 0$
(where $a_{\ell\ell}\neq 0$ since $A$ is not a GEDIRK), showing that $\dim
K_{m}\geq n_{c}$. ∎

######  Lemma 5.3.

(Dimension of $K_{m}$ for a DIRK scheme) Consider a DIRK scheme with $n_{c}$
distinct abscissas, and the corresponding space $K_{m}$ defined in (15). Then

$\dim
K_{m}\geq\min\left\\{\left\lfloor\frac{m+\kappa}{2}\right\rfloor,n_{c}\right\\}-\kappa\;,$
(40)

where $\kappa=1$ if $A$ is a GEDIRK scheme, and $\kappa=0$ otherwise.

The main idea in the proof is that the commutation property (39) implies that
the matrix $\left[A\right]_{j}$ and vectors
$\vec{\tau}^{(k)}(\left[A\right]_{j})$, define an
$\left[A\right]_{j}$-invariant space. Specifically, $\left[K_{m}\right]_{j}$
is exactly the smallest $\left[A\right]_{j}$-invariant space containing
$\vec{\tau}^{(k)}(\left[A\right]_{j})$ for $k=2,\ldots,m$. Hence, Theorem 3.5
and Lemma 5.2 apply666The new space $\mathbf{W}_{m}$ as well as the definition
$I[\cdot]$ must use the minimal polynomial of $\left[C\right]_{j}$. to
$\left[A\right]_{j}$ for each $1\leq j\leq s$, where $n_{c}$ is replaced with
the number of distinct values in $\\{c_{1},\ldots,c_{j}\\}$. The dimension of
the associated $\left[A\right]_{j}$-invariant space (for any $j$) bounds $\dim
K_{m}\geq\dim\left[K_{m}\right]_{j}$. We then pick the “worst case” $j$. For
simplicity let

$T(j):=\max\\{r\;|\>\\#\\{c_{1},\ldots,c_{r}\\}=j\\}\;,$

where $\\#$ denotes the number of distinct values in a set. Then $T(\cdot)$ is
strictly increasing, $T(1)\geq 1$, $T(n_{c})=s$, and the set
$\\{c_{1},\ldots,c_{T(j)}\\}$ contains $j$ distinct values.

###### Proof.

(of Lemma 5.3) Assume $m>1$ (it holds trivially for $m=1$). If $A$ is not a
GEDIRK, set $r:=\min\\{\lfloor\frac{m}{2}\rfloor,n_{c}\\}$. Apply Lemma 5.2 to
$\left[A\right]_{T(r)}$, which has $r$ distinct $\vec{c}$ values, to obtain
$\dim K_{m}\geq r$ (which is (40) when $\kappa=0$).

If $A$ is a GEDIRK we use the more general bound Theorem 3.5(a–c) applied to a
DIRK matrix $\left[A\right]_{T(r^{\prime})}$ with $r^{\prime}$ ($\leq m$)
distinct abscissas values:

$\displaystyle\dim K_{m}\geq\min\\{r^{\prime}-1,m-r^{\prime}\\}.$ (41)

Set $r^{\prime}:=\min\\{\frac{m}{2},n_{c}\\}$ if $m$ is even and
$r^{\prime}:=\min\\{\frac{m+1}{2},n_{c}\\}$ if $m$ is odd. Applying (41), and
observing that (in all 4 cases)
$\min\\{r^{\prime}-1,m-r^{\prime}\\}=r^{\prime}-1$ yields (40) with
$\kappa=1$. ∎

Combining Theorem 3.3 and Lemma 5.3 leads to a sharpening of Theorem 3.6 for
DIRK schemes. Theorem 3.6 still applies when $q\geq 2n_{c}-1$.

######  Theorem 5.4.

(Main Result for DIRK) Let a DIRK scheme be given with $s$ stages, $n_{c}$
distinct abscissas, order $p$, and weak stage order $q\leq 2n_{c}-1$. Then

$\left\lfloor\frac{q+\kappa}{2}\right\rfloor-\kappa+p\leq s+1-\sigma\;,$

where $\sigma=1$ if $A$ is invertible and the scheme is stiffly accurate, and
$\sigma=0$ otherwise; and $\kappa=1$ if $A$ is a GEDIRK scheme, and $\kappa=0$
otherwise.

Note that the bound (24) is used to obtain Theorem 5.4. However, in the
special case when $p=1$, (23) gives $\dim Y\geq 1$ which is tighter than (24)
($\dim Y\geq 0$) and yields a slight improvement to Theorem 5.4.

######  Remark 5.5.

(Order barriers) When $p=s+1-\sigma$, Theorem 5.4 implies that DIRK schemes
are limited to WSO $q\leq 3$, with the exception of GEDIRK schemes which are
limited to WSO $q\leq 4$. For a fixed $s$, decreasing $p$ by $1$ increases the
upper bound on $q$ by $2$. There exist schemes that satisfy these barriers
sharply (see §6). $\spadesuit$

###  5.2. Necessary Conditions on $P(x)$ for DIRK Schemes

Theorem 5.4 provides a bound on the WSO for a method with order $p$ and number
of stages $s$. However, the theorem does not explain how one might go about
constructing schemes with high WSO. We examine the solvability of the
equations $P(A)\vec{\tau}^{(k)}=0$ for the matrices $\left[A\right]_{j}$. The
results impose necessary conditions on $P(x)$ (e.g., constraints on the roots)
for the construction of high WSO schemes. We focus in this section on schemes
that are not GEDIRK schemes.

Here we write $p_{j}(x)$ to denote the minimal polynomial of
$\left[A\right]_{j}$ when $j\geq 1$ (we have $p_{j}(x)=1$ if $j=0$). The main
result is that $p_{r}$ must divide $P(x)$ when the first $r$ abscissas are
distinct.

######  Lemma 5.6.

(Necessary condition on $P(x)$ for a DIRK) Let a DIRK scheme that is not a
GEDIRK be given with WSO $q$. If the abscissas
$\\{c_{1},c_{2},\ldots,c_{r}\\}$ are distinct for
$r\leq\lfloor\frac{q}{2}\rfloor$ then $p_{r}(x)$ divides $P(x)$.

###### Proof.

Apply the map $\left[\cdot\right]_{r}$ to (20) and use (37) to obtain:

$P(\left[A\right]_{r})\;\tau^{(k)}(\left[A\right]_{r})=0,\;\textrm{for}\;k=2,\ldots,q\;.$
(42)

By Lemma 5.2 we have that $\left[K_{q}\right]_{r}=\mathbb{R}^{r}$, which when
combined with (42) implies $P(\left[A\right]_{r})=0$. Thus $p_{r}(x)$ divides
$P(x)$. ∎

######  Corollary 5.7.

For any non-confluent DIRK scheme (that is not a GEDIRK scheme), $p_{r}(x)$
divides $P(x)$ where $r=\lfloor\frac{q}{2}\rfloor$.

Example 2.4 demonstrates that the assumption of distinct abscissas in Lemma
5.6 is necessary. It provides a confluent scheme ($\vec{c}=(1,2,2,1)^{T}$)
that is $S$-irreducible and for which the minimal polynomial of $A$ (i.e.,
$p_{4}(x)=(x-1)^{3}(x-2)$) is equal to the characteristic polynomial and does
not divide $P(x)=(x-1)^{3}$.

One might wonder (I) whether the assumption of distinct abscissas in Lemma 5.6
can be relaxed; and (II) whether $p_{j}(x)$ must ever be the characteristic
polynomial of $\left[A\right]_{j}$. In general the answer to both questions is
no (see Example 2.4). However, due to the small number of variables in the top
block of $\left[A\right]_{j}$ the answer to both questions is yes when $j\leq
3$ in (I) and $j\leq 2$ in (II), provided that the scheme is $S$-irreducible.

######  Definition 5.8.

(S-reducible, Def. 12.17 in [18]) A RK method is S-reducible, if for some
partition $(S_{1},\ldots,S_{\tilde{r}})$ of $\\{1,\ldots,s\\}$ $(\tilde{r}<s)$
the indicator vectors777In a slight abuse of notation we use $\vec{S}_{j}$ to
denote both the set and the indicator vector. The notation $\vec{S}_{j}$ is
non-standard—the definition typically writes (43) as a summation.
$\vec{S}^{(m)}=\sum_{j\in S_{m}}\vec{e}_{j}$ satisfy for all $\ell$, $m$:

$(\vec{e}_{i}-\vec{e}_{j})^{T}A\vec{S}^{(m)}=0,\quad\textrm{if}\;i,j\in
S_{\ell}\;.$ (43)

An S-reducible scheme is equivalent to a smaller $\tilde{r}$-stage RK scheme
$(A^{*},\vec{b}^{*})$ where each partition $S_{j}$ is replaced by a single
stage (see Eq. (12.24) of [18, Chapt. IV.12]; the smaller scheme generates the
same stage solutions $g_{i}$ in (2a)). For S-reducible DIRKs, the new scheme
can also be made a DIRK888Formally, if we order the partitions so that
$\min\\{x\in S_{i}\\}<\min\\{y\in S_{j}\\}$ whenever $i<j$. Then for $1\leq
i\leq\tilde{r}$ define the new scheme (which is a DIRK) as
$a_{ij}^{*}=\sum_{k\in S_{j}}a_{ij}$, $b_{j}^{*}=\sum_{k\in S_{j}}b_{j}$..

######  Remark 5.9.

(A simple $S$-reducibility observation) Any DIRK scheme with $c_{1}=c_{2}$ is
$S$-reducible: apply Def. 12.17 in [18], where the partition of equivalent
stages (i.e., partition of the integers $\\{1,\ldots,s\\}$) is taken as
$S_{1}=\\{1,2\\}$ and $S_{2}=\\{3\\}$, $\ldots$, $S_{s-1}=\\{s\\}$. Then [20,
Thm. 2.2] implies that the first $2$ stages of $A$ yield the same intermediate
stage value solutions—and can be replaced by a single stage.

######  Lemma 5.10.

(A general necessary condition on $P(x)$) Let $A$ be an $S$-irreducible DIRK
scheme with invertible $\left[A\right]_{3}$ and WSO $q$.

  1. (a)

If $q>1$, then $\deg P\geq 1$ and $a_{11}$ is a root of $P$.

  2. (b)

If $q>3$, then $\deg P\geq 2$, and $a_{11},a_{22}$ are roots of $P$.

  3. (c)

If $q>5$, then $P(x)=p_{3}(x)\tilde{P}(x)$ where $p_{3}(x)$ is the minimal
polynomial of $\left[A\right]_{3}$ and is divisible by $(x-a_{11})(x-a_{22})$.

###### Proof.

(a) Since the set $\\{c_{1}\\}$ contains 1 distinct element (trivially), we
can apply Lemma 5.6 with $r=1$ when $q>1$, so that $p_{1}(x)=(x-a_{11})$
divides $P(x)$.

For (b), Remark 5.9 implies that $\\{c_{1},c_{2}\\}$ must contain two distinct
elements for $A$ to be $S$-irreducible. Applying Lemma 5.6 with $r=2$ when
$q>3$ implies that $p_{2}(x)$ divides $P(x)$. If $p_{2}(x)=(x-a_{11})$ has
only $1$ root (which must be $a_{11}$ by (a)), then
$p_{2}(\left[A\right]_{2})=0$ implies $\left[A\right]_{2}=a_{11}I$, whence
$c_{1}=c_{2}=a_{11}$ (which is a contradiction). Thus,
$p_{2}(x)=(x-a_{11})(x-a_{22})$ is the characteristic polynomial of
$\left[A\right]_{2}$.

For (c), suppose by way of contradiction that $p_{3}(x)$ does not divide
$P(x)$, so by (b), $p_{2}(x)$ is the minimal polynomial of
$\left[A\right]_{3}$ restricted to $\left[K_{q}\right]_{3}$. When $q>5$, we
claim that

$\displaystyle\left[K_{q}\right]_{3}=\textrm{span}\\{\left[\vec{e}\right]_{3},\left[\vec{c}\right]_{3}\\}\;,$
(44)

where either $c_{1}=c_{3}\neq c_{2}$ or $c_{2}=c_{3}\neq c_{1}$. The reasons
are: (i) Lemma 5.6 implies that $\vec{c}$ is confluent (with $c_{1}\neq
c_{2}$) so $n_{c}=2$ and (by the same argument in the proof)
$\dim\left[K_{q}\right]_{3}\leq 2$; (ii) if $c_{j}=0$, (44) follows
immediately from Theorem 3.5(b); (iii) If $c_{j}=0$, the construction in Lemma
5.2 implies
$\left[K_{q}\right]_{3}=\textrm{span}\\{\left[\vec{c}\right]_{3},\left[A\vec{c}\right]_{3}\\}$
which is (44) when $\left[A\right]_{3}$ is invertible (i.e., multiply by
$\left[A\right]_{3}^{-1}$). For $\left[K_{q}\right]_{3}$ to be $A$-invariant
requires

$A=\\!\begin{pmatrix}a_{11}&0&0\\\ a_{21}&a_{22}&0\\\
a_{11}\\!-\\!a_{33}\\!&0&a_{33}\end{pmatrix}\textrm{if}\;c_{1}\\!=\\!c_{3},\;A=\\!\begin{pmatrix}a_{11}&0&0\\\
a_{21}&a_{22}&0\\\
a_{21}&\\!a_{22}\\!-\\!a_{33}\\!&a_{33}\end{pmatrix}\textrm{if}\;c_{2}\\!=\\!c_{3}\;.$
(45)

In both cases $A$ is $S$-reducible with partitions $S_{1}=\\{1,3\\}$,
$S_{2}=\\{2\\}$ ($c_{1}=c_{3}$) or $S_{1}=\\{1\\}$, $S_{2}=\\{2,3\\}$
($c_{2}=c_{3}$). Hence $p_{2}(x)$ divides $p_{3}(x)$ which divides $P(x)$. ∎

######  Remark 5.11.

The proof in Lemma 5.10 follows from a more general fact (which holds for
arbitrary RK schemes) that an RK matrix is $S$-reducible if the column space
of the Vandermonde matrix
$V=(\vec{e}\,|\,\vec{c}\,|\,\cdots\,|\,\vec{c}^{\,n_{c}-1})$ is $A$-invariant.
$\spadesuit$

While this section has focused entirely on $P(x)$, the final remark concerns a
similar-in-spirit result for $Q(x)$.

######  Remark 5.12.

($a_{ss}$ is a root of $Q(x)$) A DJ-irreducible DIRK scheme has
$\vec{e}^{\,T}_{s}\vec{b}\neq 0$; otherwise, the scheme is independent of the
last stage. Since
$0=\vec{b}^{\,T}Q(A)\vec{e}_{s}=Q(a_{ss})\vec{b}^{\,T}\vec{e}_{s}$, the entry
$a_{ss}$ must be a root of $Q$, i.e., $Q(x)=(x-a_{ss})\tilde{Q}(x)$.
$\spadesuit$

##  6\. Examples

This section provides some concrete examples that highlight how the theory
established above can be used to reduce the number of degrees of freedom when
constructing RK schemes with high WSO. In Example 6.1 and Example 6.2 we
parameterize DIRK schemes with invertible $A$ and WSO $q=3$. Based on Remark
2.6 and Lemma 5.6, we can determine that certain diagonal entries of $A$
(those closest to the top left) are roots of $P(x)$, while those near the
bottom right are roots of $Q(x)$.

######  Example 6.1.

The theory enables a complete characterization of DIRK schemes with
$(s,p,q)=(2,2,3)$. Theorem 3.3 and Lemma 5.3 require: $\dim Y=1$, $\dim
K_{q}=1$. Via Lemma 5.6 the associated polynomials are $P(x)=x-a_{11}$ and
$Q(x)=x-a_{22}$. For a $p=2$ method, the orthogonality property from Lemma 4.1
implies that $Q(x)=Q_{1}(x)=x-\tfrac{1}{2}$, so $a_{22}=\tfrac{1}{2}$. The
first row of equation $P(A)\vec{\tau}^{(k)}=0$ for $k=2,3$ is automatically
satisfied; the second row yields two equations, the solution of which
determines $A$. The last order condition $\vec{b}^{\,T}\vec{e}=1$ determines
$\vec{b}$, resulting in two schemes:

$A=\begin{pmatrix}1\mp\tfrac{\sqrt{2}}{2}&0\\\
\tfrac{1}{2}\pm\tfrac{\sqrt{2}}{2}&\tfrac{1}{2}\end{pmatrix},\quad\vec{b}=\begin{pmatrix}\tfrac{1}{2}\pm\tfrac{\sqrt{2}}{4}\\\
\tfrac{1}{2}\mp\tfrac{\sqrt{2}}{4}\end{pmatrix},\quad
R(z)=\frac{1+\tfrac{z}{2}}{1-\tfrac{z}{2}}\;.$ (46)

Either choice of signs leads to a $(2,2,3)$ scheme. In both cases,
$\vec{b}^{\,T}A=\tfrac{1}{2}\vec{b}^{\,T}$ is a left eigenvector of $A$ and
$K_{3}=\textrm{span}\\{2\vec{c}-\vec{e}\,\\}=\mathbf{W}_{3}$.

######  Example 6.2.

(DIRK schemes with $(s,p,q)=(3,3,3)$) The chosen values of $(s,p,q)$ satisfy
Theorem 5.4 sharply and require that $\dim Y=2$ and $\deg P=\dim K_{q}=1$ with
$P(x)=x-a_{11}$, $Q(x)=(x-a_{22})(x-a_{33})$. Since $\left[A\right]_{2}$
satisfies (42) it must be a scalar multiple (call the ratio $a$) of the
Butcher matrix in (46). The polynomial $Q(x)=Q_{2}(x)+\alpha_{1}Q_{1}(x)$ then
has $a_{22}=a/2$ as a root so that
$\alpha_{1}=-Q_{2}(\tfrac{a}{2})/Q_{1}(\tfrac{a}{2})$, which uniquely
determines (i) $a_{33}$ as the second root of $Q(x)$; and (ii) the stability
function in terms of $a$. The matrix $A$ and $R(z)$ are constrained to be:

$A=\begin{pmatrix}(1\mp\tfrac{\sqrt{2}}{2})a&0&0\\\
(\tfrac{1}{2}\pm\tfrac{\sqrt{2}}{2})a&\tfrac{1}{2}a&0\\\
a_{31}&a_{32}&\tfrac{3a-2}{6(a-1)}\end{pmatrix},\quad
R(z)=\frac{(1+6\alpha_{1})z^{2}+(6+12\alpha_{1})z+12}{(1-6\alpha_{1})z^{2}+(12\alpha_{1}-6)z+12}\;.$

The values $(a_{31},a_{32})$ can be parameterized in terms of $a$ via the
$3$rd row equations of $P(A)\vec{\tau}^{(j)}=0$ for $j=2,3$. For each choice
of $\pm$ in $A$, we have the following. There are $3$ solution branches for
$(a_{31},a_{32})$, two of which yield S-reducible schemes with the structures
given in (45). The one irreducible solution branch has that $\vec{\tau}^{(2)}$
and $\vec{\tau}^{(3)}$ are parallel. The vector $\vec{b}$ is determined via
$\vec{b}^{\,T}(\vec{e},\vec{c},\vec{\tau}^{(2)})=(1,\tfrac{1}{2},0)$. By
construction, $\vec{b}^{\,T}Q(A)=0$ so that $\vec{b}^{\,T}\vec{e}=1$ and
$\vec{b}^{\,T}\vec{c}=\tfrac{1}{2}$ ensure
$\vec{b}^{\,T}A^{2}\vec{e}=\tfrac{1}{6}$. Finally, combining
$\vec{b}^{\,T}(A\vec{\tau}^{(1)}+\vec{\tau}^{(2)})=0$ with
$\vec{b}^{\,T}A^{2}\vec{e}=\tfrac{1}{6}$ yields
$\vec{b}^{\,T}\vec{c}^{\,2}=\tfrac{1}{3}$.

As expected, with $s=3$ there is no choice of $a$ that satisfies the
additional conditions imposed by order $p=4$ or $q=4$. Finally,
$0<a<\tfrac{2}{3}$ or $a>1$ is required for $A$ to have positive eigenvalues
and $R(z)$ to have all its poles in the right half plane.

######  Example 6.3.

The bound in Theorem 3.6 can be sharp for both stiffly accurate schemes and
EDIRKs. Specifically, stiffly accurate DIRK schemes with $(s,p,q)=(4,3,3)$
were constructed in an ad-hoc fashion in [22], where a-posteriorily $s$ is
observed to be sharp. Setting $\dim K_{q}=0$ in Lemma 5.3, the fact
$\tilde{q}\leq q$ recovers the known result that stage order is limited to
$\tilde{q}=2$ for EDIRKs and $\tilde{q}=1$ for DIRKs.

######  Example 6.4.

(Schemes with high stage order) The Gauss-Legendre RK methods satisfy the
bound in Theorem 3.6 sharply since $p=2s$, $q=p$ (so that $\dim(K_{q})=0$) and
$n_{c}=s$ (and are not stiffly accurate). As an example of Corollary 3.4, we
obtain the bound $p\leq 2s$, well known already from the theory of numerical
quadrature.

######  Remark 6.5.

(Guide to constructing DIRK schemes with high WSO) For schemes with $q=p$, one
can take $r=\lfloor q/2\rfloor$, set $P(x)=(x-a_{11})\cdots(x-a_{rr})$ and
then solve (20) as a sufficient condition (and use Theorem 5.4 to guide the
number of stages). $\spadesuit$

##  7\. Conclusions and Outlook

Weak stage order can be viewed as a geometric condition that, when satisfied,
can remove order reduction in Runge-Kutta schemes applied to linear problems
with time-independent operators. The general theory of WSO provided here
relates geometric objects (WSO invariant subspaces) to associated algebraic
objects (minimal polynomials). This relationship allows us to establish order
barriers for WSO that generalize the well-known bounds on the RK order $p$ in
terms of $s$ (both for general RK schemes and DIRKs). Along the way, we also
provide new formulas for the RK stability function in terms of a family of
polynomials which are “orthogonal” with respect to a linear functional. The
new necessary conditions show how one needs to split the spectrum of $A$ into
the roots of $P(x)$ and $Q(x)$—which is of practical use when constructing
high WSO schemes (i.e., beyond 3). Indeed, the necessary conditions were used
in the companion work [6] in the construction of new DIRK schemes with WSO 4
and 5 (and satisfy the full set of order conditions).

The results presented here give rise to several research directions. First, is
Theorem 5.4 sharp for all $p,q,s$? And, is it further constrained by the non-
tall tree order conditions? There is also the related (practical) issue of
constructing DIRK schemes with $q=p$ (or $q=p-1$) and the fewest stages
$s\sim\frac{2}{3}p$ satisfying Theorem 5.4 (or better yet, analytically
parameterizing such schemes). While the main results in this work do apply to
all (E)DIRKs, DIRKs with additional constraints such as SDIRKs or ERKs are
expected to further impact the bounds in Theorem 5.4.

## Acknowledgments

This material is based upon work supported by the National Science Foundation
under Grants No. DMS–2012271 (Biswas, Seibold), No. DMS–1952878 (Seibold), and
No. DMS–2012268 (Shirokoff). D. Shirokoff gratefully acknowledges the
hospitality of Robert Pego and the Dept. of Mathematical Sciences at Carnegie
Melon University where portions of this work were completed.

##  Appendix A Proofs of Theorems 2.2 and 2.3

We collect the linear algebra proofs here, which are adapted from various
materials (e.g., [31, Chap. 2.6]).

###### Proof.

(Theorem 2.2, $A$-invariant subspaces) Let $\mathbf{U}\in\mathbb{R}^{s\times
d}$ form an orthogonal basis for $U$. Since $U$ is $A$-invariant, there is a
square matrix $A_{UU}:=\mathbf{U}^{T}A\mathbf{U}\in\mathbb{R}^{d\times d}$
such that (each column of $\mathbf{U}$ is mapped back into the column space of
$\mathbf{U}$)

$A\mathbf{U}=\mathbf{U}A_{UU},$ (47)

For any polynomial $p(x)$, (47) then implies that
$p(A)\mathbf{U}=\mathbf{U}p(A_{UU})$. Thus, $p(x)$ in (18) is the minimal
polynomial of $A_{UU}$ (which has real coefficients since $A_{UU}$ is real)
because it satisfies $p(A)\mathbf{U}=0$ and any other smaller degree
polynomial fails to satisfy $p(A)\mathbf{U}=\mathbf{U}p(A_{UU})\neq 0$.

For (a), $p(x)$ divides the characteristic polynomial of $A_{UU}$, which
divides $\textrm{char}_{A}(x)$; (c) also follows since
$\textrm{deg}\;p\leq\textrm{deg}\;\textrm{char}_{A_{UU}}(x)=\textrm{dim}(U)$.

For (b), every root $\lambda$ of $p(x)$ is an eigenvalue of $A_{UU}$ and hence
of $A$.

For (d), set $\ell=\textrm{deg}\;p$. The condition $p(A)\vec{v}=0$ (where
$\vec{v}$ is the vector defining $U$) implies that $A^{\ell}\vec{v}$ is a
linear combination of $\vec{v},\ldots,A^{\ell-1}\vec{v}$, which is only
possible if $\textrm{dim}(U)\leq\textrm{deg}\;p$. Combined with part (c), this
shows $\dim U=\deg p$.

For (e), when $U=\\{\vec{v},A\vec{v},\ldots,A^{d-1}\vec{v}\\}$, condition (18)
implies that $p(A)\vec{v}=0$ since $\vec{v}\in U$. Conversely, any element
$\vec{u}\in U$ is a linear combination of $A^{j}\vec{v}$ for $j=0,\ldots,d-1$.
But then, $p(A)\vec{v}=0$ implies (since $A^{j}$ and $p(A)$ commute) that
$p(A)A^{j}\vec{v}=0$ for any $j$, and hence $p(A)\vec{u}=0$ for any
$\vec{u}\in U$. ∎

###### Proof.

(Theorem 2.3, Left and right orthogonal $A$-invariant subspaces). Similar to
the proof of Theorem 2.2, introduce the orthogonal matrix
$\mathbf{O}=(\mathbf{U}\ |\ \mathbf{W}\ |\ \mathbf{V})\in\mathbb{R}^{s\times
s}$, where the columns of $\mathbf{U}\in\mathbb{R}^{s\times d_{u}}$ and
$\mathbf{V}\in\mathbb{R}^{s\times d_{v}}$ form an orthogonal basis for $U$ and
$V$, respectively (the columns of $\mathbf{W}$ span the remaining space). Now
$A$ block-diagonalizes as

$\mathbf{O}^{T}A\mathbf{O}=\begin{pmatrix}A_{UU}&A_{UW}&A_{UV}\\\
0&A_{WW}&A_{WV}\\\ 0&0&A_{VV}\end{pmatrix},$

where $A_{\Sigma\Theta}=\Sigma^{T}A\Theta$ where
$\Sigma,\Theta\in\\{\mathbf{U},\mathbf{V},\mathbf{W}\\}$. Hence,

$\textrm{char}_{A}(x)=\textrm{char}_{A_{UU}}(x)\;\textrm{char}_{A_{WW}}(x)\;\textrm{char}_{A_{VV}}(x)\;.$
(48)

Via Theorem 2.2 the polynomial $p(x)$ is the minimal polynomial of $A_{UU}$.
Theorem 2.2 applies to $A^{T}$ with space $\mathbf{V}$, so that $q(x)$ is the
minimal polynomial of $A_{VV}^{T}$ (and is the same as the minimal polynomial
of $A_{VV}$). Hence, $p(x)q(x)$ divides
$\textrm{char}_{A_{WW}}(x)\textrm{char}_{A_{VV}}(x)$, which, from (48),
divides $\textrm{char}_{A}(x)$. ∎

It is generally not true that $p(x)q(x)$ divides the minimal polynomial of
$A$.

##  Appendix B Proofs of Hankel Matrix Determinant and Orthogonal Polynomial
Recurrence Relation

Here we provide details for the determinant computation ($m\in\mathbb{Z}_{\geq
0}$) of

$\mathbf{H}_{n}(m)=\begin{pmatrix}\tfrac{1}{m!}&\tfrac{1}{(m+1)!}&\cdots&\tfrac{1}{(m+n-1)!}\\\
\tfrac{1}{(m+1)!}&\ddots&&\tfrac{1}{(m+n)!}\\\ \vdots&&\ddots&\vdots\\\
\tfrac{1}{(m+n-1)!}&\cdots&\cdots&\tfrac{1}{(m+2n-2)!}\end{pmatrix}\in\mathbb{R}^{n\times
n}.$ (49)

Let $M_{n}(m):=\det(\mathbf{H}_{n}(m))$. A formula (without proof) for
$M_{n}(m)$ is given in [24]:

$\displaystyle M_{n}(m)$
$\displaystyle=\sigma(n)\frac{c(n)c(m+n-1)}{c(m+2n-1)},\quad\textrm{where}\;c(n):=\Pi_{i=1}^{n-1}i!\;,$
(50) $\displaystyle\textrm{and }\sigma(n)$
$\displaystyle:=\left\\{\begin{array}[]{cr}\phantom{-}1,&\textrm{if
}r=0\textrm{ or }1,\\\ -1,&\textrm{if }r=2\textrm{ or
}3,\end{array}\right.\textrm{ where }n\equiv r\mod 4\;.$ (53)

Since we could not find a proof of (50) in the literature, we provide a brief
one here. Following the approach in [23, Method 2 in §4], for any square
matrix $A$, a determinant formula due to Jacobi reads:

$\displaystyle\det A\cdot\det A_{1,n}^{1,n}=\det A_{1}^{1}\cdot\det
A_{n}^{n}-\det A_{1}^{n}\cdot\det A_{n}^{1}\;,$ (54)

where $A_{i}^{j}$ denotes the submatrix of $A$ with row $i$ and column $j$
removed. Applying (54) to the matrix $M_{n}(m)$ and using the symmetry of the
Hankel matrix yields the recursion:

$M_{n}(m)M_{n-2}(m+2)=M_{n-1}(m+2)M_{n-1}(m)-M_{n-1}(m+1)^{2}\;.$ (55)

Formula (55) relates $M_{n}(m)$ to matrices of size $n-1$ and $n-2$, and can
be used to prove (50) via induction. Note that for $M_{1}(m)=1/m!$ and
$M_{2}(m)=-[(m+1)!(m+2)!]^{-1}$, formula (50) is readily verified (and holds
$\forall m\geq 0$). Assume (50) holds for $1\leq k<n$. We can then divide (55)
through by $M_{n-1}(m+1)M_{n-2}(m+2)$ (which are non-zero by assumption) to
obtain:

$\tfrac{M_{n}(m)}{M_{n-1}(m+1)}=\tfrac{M_{n-1}(m+2)M_{n-1}(m)}{M_{n-2}(m+2)M_{n-1}(m+1)}-\tfrac{M_{n-1}(m+1)}{M_{n-2}(m+2)}\;.$
(56)

It is then a matter of substituting the formulas from (50) into (56) to verify
(after several lines of factorial cancellations) that the following
holds999Note that $\sigma(n)/\sigma(n-1)$ always has opposite sign to
$\sigma(n-1)/\sigma(n-2)$.:

$\tfrac{(n-1)!}{(m+2n-2)!}=-\left(\tfrac{(n-2)!(m+n-1)!(m+2n-3)!}{(m+2n-2)!(m+2n-3)!(m+n-2)!}-\tfrac{(n-2)!}{(m+2n-3)!}\right)\;.$

Since (for every $n$) the Hankel matrix determinants $M_{n}(m)$ do not vanish,
the entries of $\mathbf{H}_{n}(m)$ define a _quasi-definite linear functional_
which then have associated orthogonal polynomials. Monic orthogonal
polynomials (always) satisfy a three term linear recurrence of the general
form

$Q_{n+1}(x)=(x+\gamma_{n})Q_{n}(x)+\beta_{n}Q_{n-1}(x)\;,$ (57)

where $\mathbf{H}_{n}(m)$ is initialized via $Q_{0}(x)=1$,
$Q_{1}(x)=x-1/(m+1)$. We now obtain formulas for $\beta_{n}$, $\gamma_{n}$ via
determinant computations. From [23, Thm. 29] we have:

$\displaystyle\beta_{n}(m)=-\tfrac{M_{n-1}(m)}{M_{n}(m)}\cdot\tfrac{M_{n+1}(m)}{M_{n}(m)}=\tfrac{n(m+n-1)}{(m+2n)(m+2n-1)^{2}(m+2n-2)}\;,$
$\displaystyle\textrm{When }m=1:\quad\beta_{n}=\tfrac{1}{4(4n^{2}-1)}\;.$

The $\gamma_{n}$ depend only on the second leading term $\lambda_{n}(m)$ of
$Q_{n}(x)=x^{n}+\lambda_{n}(m)x^{n-1}+\textrm{low order terms}$. That is,
evaluating the coefficient of $x^{n}$ in equation (57) yields
$\gamma_{n}(m)=\lambda_{n+1}(m)-\lambda_{n}(m)$.

Orthogonality of $Q_{j}(x)$ demands that the coefficients
$\vec{\lambda}=(\lambda_{0}(m),\ldots,\lambda_{n}(m))^{T}$ satisfy
$\mathbf{H}_{n}(m)\vec{\lambda}=\vec{h}_{n}$ with
$\vec{h}_{n}=-(1/(m+n)!,\ldots 1/(m+2n+1)!)^{T}$. Using Cramer’s rule for
$\lambda_{n}(m)$ and (54), one obtains the two term recursion for
$\lambda_{n}(m)$

$\displaystyle\lambda_{n}(m)M_{n}(m)M_{n-2}(m+2)$
$\displaystyle=\lambda_{n-1}(m+2)M_{n-1}(m+2)M_{n-1}(m)$
$\displaystyle-\lambda_{n-1}(m+1)M_{n-1}(m+1)^{2}\;,$

which simplifies to

$\lambda_{n}(m)=\tfrac{(m+2n-2)}{(n-1)}\lambda_{n-1}(m+1)-\tfrac{(m+n-1)}{(n-1)}\lambda_{n-1}(m+2)\;.$

Induction then shows that

$\lambda_{n}(m)=-\tfrac{n}{m+2n-1},\;\forall n\geq 1,m\geq 0\;.$

The orthogonal polynomials in this paper use the values $m=0,1$. In the case
when $m=1$, the coefficient $\lambda_{n}(1)=-\frac{1}{2}$ is constant, in
which case $\gamma_{n}(1)=0$, $\forall n\geq 1$.

## References

  * [1] S. Abarbanel, D. Gottlieb, and M. H. Carpenter.  On the removal of boundary errors caused by Runge–Kutta integration of nonlinear partial differential equations.  SIAM Journal on Scientific Computing, 17(3):777–782, 1996.
  * [2] I. Alonso-Mallo.  Runge-Kutta methods without order reduction for linear initial boundary value problems.  Numer. Math., 91(4):577–603, 2002.
  * [3] I. Alonso-Mallo and B. Cano.  Avoiding order reduction of Runge–Kutta discretizations for linear time-dependent parabolic problems.  BIT Numerical Mathematics, 44(1):1–20, 2004.
  * [4] I Alonso-Mallo and C Palencia.  Optimal orders of convergence for runge–kutta methods and linear, initial boundary value problems.  Applied numerical mathematics, 44(1-2):1–19, 2003.
  * [5] Sheldon Axler.  Linear Algebra Done Right.  Undergraduate Texts in Mathematics. Springer, third edition, 2015.
  * [6] A. Biswas, D. Ketcheson, B. Seibold, and D. Shirokoff.  Design of DIRK schemes with high weak stage order, 2022.  Under review.
  * [7] Abhijit Biswas.  Structure-Optimized and Order-Preserving Methods for Advection and Time Stepping.  PhD thesis, Temple University, 2021.
  * [8] K. Burrage, W.H. Hundsdorfer, and J.G. Verwer.  A study of B-Convergence of Runge-Kutta methods.  Computing, 36:17–34, 1986.
  * [9] K. Burrage and L. Petzold.  On order reduction for Runge-Kutta methods applied to differential/algebraic systems and to stiff systems of ODEs.  SIAM J. Numer. Anal., 27(2):447–456, 1990.
  * [10] J. C. Butcher and N. Goodwin.  Numerical methods for ordinary differential equations, volume 2.  Wiley Online Library, 2008.
  * [11] M. Calvo and C. Palencia.  Avoiding the order reduction of Runge-Kutta methods for linear initial boundary value problems.  Math. Comp., 71(240):1529–1543, 2002.
  * [12] M.P. Calvo and C Palencia.  Avoiding the order reduction of Runge-Kutta methods for linear initial boundary value problems.  Math. Comput., 71(240):1529–1543, 2002.
  * [13] M. H. Carpenter, D. Gottlieb, S Abarbanel, and W-S Don.  The theoretical accuracy of Runge-Kutta time discretizations for the initial boundary value problem: a study of the boundary error.  SIAM J. Sci. Comput., 16(6):1241–1252, 1995.
  * [14] M. Crouzeix and P.-A. Raviart.  Approximation des problèmes d’évolution.  Lecture Notes, University Rennes, 1980.
  * [15] K. Dekker and J. G. Verwer.  Stability of Runge-Kutta methods for stiff nonlinear differential equations, volume 2 of CWI Monographs.  North-Holland Publishing Co., 1984.
  * [16] J. García-Ardila, F. Marcellán, and M. Marriaga.  From standard orthogonal polynomials to Sobolev orthogonal polynomials: The role of semiclassical linear functionals.  In Foupouagnigni M., Koepf W. (eds) Orthogonal Polynomials. AIMSVSW 2018., pages 245–292. Birkhäuser, Cham, 2020.
  * [17] E. Hairer, S. P. 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.  Solving ordinary differential equations II.  Springer Berlin Heidelberg, 1996.
  * [19] W. H. Hundsdorfer.  Stability and B-Convergence of linearly implicit Runge-Kutta methods.  Numer. Math., 50:83–95, 1986.
  * [20] WH Hundsdorfer and MN Spijker.  A note on B-stability of Runge-Kutta methods.  Numer. Math., 36(3):319–331, 1980.
  * [21] M. Ismail.  Classical and Quantum Orthogonal Polynomials in One Variable.  Cambridge University Press, New York, 2005.
  * [22] D. Ketcheson, B. Seibold, D. Shirokoff, and D. Zhou.  DIRK schemes with high weak stage order.  In Lecture Notes in Computational Science and Engineering, pages 453–463. Proceedings of ICOSAHOM 2018, Springer, 2020.
  * [23] C. Krattenthaler.  Advanced determinant calculus: A complement.  Linear Algebra Appl., 411:68–166, 2005.
  * [24] S. Landsburg.  Some Hankel determinants.  MathOverflow.  https://mathoverflow.net/q/249838 (version: 2017-04-13).
  * [25] Ch. Lubich and A. Ostermann.  Runge-Kutta approximation of quasi-linear parabolic equations.  Math. Comp., 64(210):601–627, 1995.
  * [26] JI Montijano.  Estudio de los metodos SIRK para la resolucion numérica de ecuaciones diferenciales de tipo stiff.  PhD thesis, Thesis, University of Zaragoza, 1983.
  * [27] S. P. Nørsett and A. Wolfbrandt.  Attainable order of rational approximations to the exponential function with only real poles.  BIT Numer. Math., 17:200–208, 1977.
  * [28] A. Ostermann and M. Roche.  Runge-Kutta methods for partial differential equations and fractional orders of convergence.  Math. Comp., 59(200):403–420, 1992.
  * [29] A. Ostermann and M. Roche.  Rosenbrock methods for partial differential equations and fractional orders of convergence.  SIAM J. Numer. Anal., 30(4):1084–1098, 1993.
  * [30] A. Ostermann and M. Thalhammer.  Convergence of Runge-Kutta methods for nonlinear parabolic equations.  Appl. Numer. Math., 42(1):367–380, 2002.
  * [31] P. Petersen.  Linear Algebra.  Springer, 2021.
  * [32] P.J. Prince and J.R. Dormand.  High order embedded Runge-Kutta formulae.  J. Comput. Appl. Math., 7(1):67–75, March 1981.
  * [33] A. Prothero and A. Robinson.  On the stability and accuracy of one-step methods for solving stiff systems of ordinary differential equations.  Math. Comput., 28(125):145–162, 1974.
  * [34] J. Rang.  An analysis of the Prothero–Robinson example for constructing new DIRK and ROW methods.  J. Comput. Appl. Math, 262:105–114, 2014.
  * [35] J. Rang.  The Prothero and Robinson example: Convergence studies for Runge-Kutta and Rosenbrock-Wanner methods.  Appl. Numer. Math., 108:37–56, 2016.
  * [36] S. Roberts.  Multimethods for the Efficient Solution of Multiscale Differential Equations.  PhD thesis, Virginia Polytechnic Institute and State University, 2021\.
  * [37] S. Roberts and A. Sandu.  Eliminating order reduction on linear, time-dependent ODEs with GARK methods, 2022.  arXiv:2201.07940.
  * [38] R. R. Rosales, B. Seibold, D. Shirokoff, and D. Zhou.  Spatial manifestations of order reduction in Runge-Kutta methods for initial boundary value problems, 2017.  arXiv:1712.00897.
  * [39] Steven J. Ruuth and Raymond J Spiteri.  High-order strong-stability-preserving Runge-Kutta methods with downwind-biased spatial discretizations.  SIAM J. Numer. Anal., 42:974–996, 2004.
  * [40] J. M. Sanz-Serna, J. G. Verwer, and W.-H. Hundsdorfer.  Convergence and order reduction of Runge-Kutta schemes applied to evolutionary problems in partial differential equations.  Numerische Mathematik, 50(4):405–418, 1986.
  * [41] S. Scholz.  Order barriers for the B-convergence of ROW methods.  Computing, 41(3):219–235, 1989.
  * [42] J.G. Verwer and J.M. Sanz-Serna.  Convergence of method of lines approximations to partial differential equations.  Computing, 33:297–313, 1984.


# A Proof that Coarse Correlated Equilibrium Implies Nash Equilibrium in Two-
Player Zero-Sum Games

Revan MacQueen

(2023)

###### Abstract

We give a simple proof of the well-known result that the marginal strategies
of a coarse correlated equilibrium form a Nash equilibrium in two-player zero-
sum games. A corollary of this fact is that no-external-regret learning
algorithms that converge to the set of coarse correlated equilibria will also
converge to Nash equilibria in two-player zero-sum games. We show an
approximate version: that $\epsilon$-coarse correlated equilibria imply
$2\epsilon$-Nash equilibria.

##  1 Background

It is well-known that the marginal strategies of a coarse correlated
equilibrium (CCE) are Nash equilibria in two-player zero-sum games. We present
a simple proof of this fact, which may be useful to others.

This fact is useful, since a no-external-regret learning algorithm’s average
strategy profile will approach the set of CCE (Greenwald and Jafari, 2003),
and so will compute an approximate Nash equilibrium in two-player zero-sum
games. Hence, the marginal strategies of a CCE in two-player zero-sum games
have the desirable properties of Nash equilibrium, including bounded
exploitability (Fudenberg and Tirole, 1991).

######  Definition 1.1.

A normal form game $G$ is a 3 tuple $G=(N,\mathrm{P},u)$ where $N$ is a set of
players, $\mathrm{P}=\bigtimes_{i\in N}\mathrm{P}_{i}$ is a joint pure
strategy space where $\mathrm{P}_{i}$ is a set of _pure strategies_ for player
$i$. $u=(u_{i})_{i\in N}$ is a set of utility functions where
$u_{i}:\mathrm{P}\to\mathbb{R}$.

Pure strategies are deterministic choices of actions in the game. In a normal-
form game, each agent simultaneously chooses their pure strategy
$\rho_{i}\in\mathrm{P}_{i}$; we call a joint selection of pure strategies
$\rho\in\mathrm{P}$ a _pure strategy profile_. The pure strategy profile
determines the payoff (or utility) to player via the utility function.

Players may also randomize their actions through the use of a _mixed strategy_
: a probability distribution $s_{i}$ over $i$’s pure strategies. Let
$S_{i}=\Delta(\mathrm{P}_{i})$ be the set of player $i$’s mixed strategies
(where $\Delta(X)$ denotes the set of probability distributions over a domain
$X$), and let $S=\bigtimes_{i\in N}S_{i}$ be the set of mixed strategy
profiles.We overload the definition of utility function to accept mixed
strategies as follows:

$\displaystyle u_{i}(s)=\sum_{\rho\in\mathrm{P}}\left(\prod_{i\in
N}s_{i}(\rho_{i})\right)u_{i}(\rho)$

In a two player game, when one player plays a pure strategy while the other
plays a mixed strategy, we write

$\displaystyle
u_{i}(\rho_{i},s_{-i})=\sum_{\rho_{-i}\in\mathrm{P}_{-i}}s_{-i}(\rho_{-i})u_{i}(\rho_{i},\rho_{-i})$

Where $-i$ is the other player from player $i$ in a two-player game.

######  Definition 1.2.

A two-player game $G=(\\{i,-i\\},\mathrm{P},u)$ is zero-sum111Note that two-
player zero-sum games are no less general than two-player constant-sum games,
since we can turn any two-player constant-sum game into a strategically
identical zero-sum game by either subtracting or adding a constant to each
player’s utility. if

$\displaystyle u_{i}(\rho)+u_{-i}(\rho)=0\quad\forall\rho\in\mathrm{P}$

######  Definition 1.3 ($\epsilon$-Nash Equilibrium).

A strategy profile $s$ is an $\epsilon$-Nash equilibrium if $\forall i\in N$,
we have

$\displaystyle
u_{i}(s^{\prime}_{i},s_{-i})-u_{i}(s_{i},s_{-i})\leq\epsilon\quad\forall
s^{\prime}_{i}\in S_{i}$

######  Definition 1.4 ($\epsilon$-CCE (Moulin and Vial, 1978)).

We say $\mu\in\Delta(\mathrm{P})$ is an $\epsilon$-CCE if $\forall i\in N$, we
have

$\displaystyle\mathbb{E}_{\rho\sim\mu}\left[u_{i}(s_{i}^{\prime},\rho_{-i})-u_{i}(\rho)\right]\leq\epsilon\quad\forall
s_{i}^{\prime}\in S_{i}$

######  Definition 1.5 (Marginal strategy).

Given $\mu\in\Delta(\mathrm{P})$, let $s_{i}^{\mu}$ be the _marginal strategy_
for $i$, where
$s_{i}^{\mu}(\rho_{i})\doteq\sum_{\rho_{-i}\in\mathrm{P}_{-i}}\mu(\rho_{i},\rho_{-i})$.
Let $s^{\mu}$ be a _marginal strategy profile_ , where each player plays their
marginal strategy from $\mu$.

##  2 The Proof

Here we show a simple proof that CCE imply Nash in two-player zero-sum games.
We begin showing that if a player $i$ deviates from their CCE recommendations,
if $-i$ continues to play their CCE recommendations, this is equivalent to
$-i$ playing their marginal strategy for the CCE.

######  Lemma 2.1.

If $\mu\in\Delta(\mathrm{P})$, then in a two-player game, for any deviation
$s_{i}$, we have
$\mathbb{E}_{\rho\sim\mu}\left[u_{i}(s_{i},\rho_{-i})\right]=u_{i}(s_{i},s^{\mu}_{-i})$
.

###### Proof.

We have

$\displaystyle\mathbb{E}_{\rho\sim\mu}\left[u_{i}(s_{i},\rho_{-i})\right]$
$\displaystyle=\sum_{\rho_{i}\in\mathrm{P}_{i}}\sum_{\rho_{-i}\in\mathrm{P}_{-i}}s_{i}(\rho_{i})\mu_{-i}(\rho_{-i})u_{i}(\rho_{i},\rho_{-i}),$

where
$\mu_{-i}(\rho_{-i})\doteq\sum_{\rho_{-i}\in\mathrm{P}_{-i}}\mu(\rho_{i},\rho_{-i})$.
Note, however, that $\mu_{-i}(\rho_{-i})=s^{\mu}_{-i}(\rho_{-i})$ so

$\displaystyle\mathbb{E}_{\rho\sim\mu}\left[u_{i}(s_{i},\rho_{-i})\right]$
$\displaystyle=\sum_{\rho_{i}\in\mathrm{P}_{i}}\sum_{\rho_{-i}\in\mathrm{P}_{-i}}s_{i}(\rho_{i})s^{\mu}_{-i}(\rho_{-i})u_{i}(\rho_{i},\rho_{-i})$
$\displaystyle=u_{i}(s_{i},s^{\mu}_{-i}).$

∎

######  Proposition 2.2.

If $\mu$ is an $\epsilon$-CCE of a two-player constant-sum game, then
$|\mathbb{E}_{\rho\sim\mu}\left[u_{i}(\rho)\right]-u_{i}(s^{\mu})|\leq\epsilon$

###### Proof.

Suppose not, then
$\left|\mathbb{E}_{\rho\sim\mu}\left[u_{i}(\rho)\right]-u_{i}(s^{\mu})\right|>\epsilon$
which means either

$\displaystyle
u_{i}(s^{\mu})-\mathbb{E}_{\rho\sim\mu}\left[u_{i}(\rho)\right]>\epsilon$ (1)

or

$\displaystyle\mathbb{E}_{\rho\sim\mu}\left[u_{i}(\rho)\right]-u_{i}(s^{\mu})>\epsilon.$
(2)

Consider (1). By Lemma 2.1 we have
$u_{i}(s^{\mu})=\mathbb{E}_{\rho\sim\mu}\left[u_{i}(s_{i}^{\mu},\rho_{-i})\right]$,
which means

$\displaystyle\mathbb{E}_{\rho\sim\mu}\left[u_{i}(s_{i}^{\mu},\rho_{-i})\right]-\mathbb{E}_{\rho\sim\mu}\left[u_{i}(\rho)\right]>\epsilon,$

which contradicts the fact that $\mu$ is an $\epsilon$-CCE, since
$s^{\mu}_{i}$ is a deviation that is more than $\epsilon$-profitable for
player $i$. Next, consider (2); since the game is two-player zero-sum, we
have:

$\displaystyle\mathbb{E}_{\rho\sim\mu}\left[u_{i}(\rho)\right]-u_{i}(s^{\mu})>\epsilon$
$\displaystyle\implies$
$\displaystyle-\mathbb{E}_{\rho\sim\mu}\left[u_{-i}(\rho)\right]-(-u_{-i}(s^{\mu}))>\epsilon$
$\displaystyle\implies$ $\displaystyle
u_{-i}(s^{\mu})-\mathbb{E}_{\rho\sim\mu}\left[u_{-i}(\rho)\right]>\epsilon.$

At this point we may repeat the steps above to show that $\mu$ is not an
$\epsilon$-CCE, since $s^{\mu}_{-i}$ is a deviation that is more than
$\epsilon$-profitable for player $-i$. ∎

######  Proposition 2.3.

If $\mu$ is an $\epsilon$-CCE of a two-player constant-sum game $G$, then
$s^{\mu}$ is a $2\epsilon$-Nash equilibrium.

###### Proof.

Choose player $i$ arbitrarily. Either
$u_{i}(s^{\mu})\geq\mathbb{E}_{\rho\sim\mu}\left[u_{i}(\rho)\right]$ or
$u_{i}(s^{\mu})<\mathbb{E}_{\rho\sim\mu}\left[u_{i}(\rho)\right]$. Consider
the first case. Starting from the definition of $\epsilon$-CCE:

$\displaystyle\mathbb{E}_{\rho\sim\mu}\left[u_{i}(\rho_{i}^{\prime},\rho_{-i})\right]-\mathbb{E}_{\rho\sim\mu}\left[u_{i}(\rho)\right]\leq\epsilon\quad\forall\rho_{i}^{\prime}\in\mathrm{P}_{i}.$

But since $u_{i}(s^{\mu})\geq\mathbb{E}_{\rho\sim\mu}\left[u_{i}(\rho)\right]$
we have

$\displaystyle\mathbb{E}_{\rho\sim\mu}\left[u_{i}(\rho_{i}^{\prime},\rho_{-i})\right]-u_{i}(s^{\mu})\leq\epsilon\quad\forall\rho_{i}^{\prime}\in\mathrm{P}_{i}.$

Which by Lemma 2.1 means

$\displaystyle
u_{i}(\rho_{i}^{\prime},s^{\mu}_{-i})-u_{i}(s^{\mu})\leq\epsilon\quad\forall\rho_{i}^{\prime}\in\mathrm{P}_{i}.$

Thus $s^{\mu}$ is an $\epsilon$-Nash. Next, suppose
$u_{i}(s^{\mu})<\mathbb{E}_{\rho\sim\mu}\left[u_{i}(\rho)\right]$. By
Proposition 2.2 we have,

$\displaystyle\mathbb{E}_{\rho\sim\mu}\left[u_{i}(\rho)\right]-u_{i}(s^{\mu})\leq\epsilon$
(3) $\displaystyle\implies$
$\displaystyle\mathbb{E}_{\rho\sim\mu}\left[u_{i}(\rho)\right]\leq
u_{i}(s^{\mu})+\epsilon.$ (4)

Then, starting from the definition of $\epsilon$-CCE and applying (4),

$\displaystyle\mathbb{E}_{\rho\sim\mu}\left[u_{i}(\rho_{i}^{\prime},\rho_{-i})\right]-\mathbb{E}_{\rho\sim\mu}\left[u_{i}(\rho)\right]\leq\epsilon\quad\forall\rho_{i}^{\prime}\in\mathrm{P}_{i}$
$\displaystyle\implies$
$\displaystyle\mathbb{E}_{\rho\sim\mu}\left[u_{i}(\rho_{i}^{\prime},\rho_{-i})\right]-(u_{i}(s^{\mu})+\epsilon)\leq\epsilon\quad\forall\rho_{i}^{\prime}\in\mathrm{P}_{i}$
$\displaystyle\implies$
$\displaystyle\mathbb{E}_{\rho\sim\mu}\left[u_{i}(\rho_{i}^{\prime},\rho_{-i})\right]-u_{i}(s^{\mu})\leq
2\epsilon\quad\forall\rho_{i}^{\prime}\in\mathrm{P}_{i}.$

By Lemma 2.1 we have

$\displaystyle u_{i}(\rho_{i}^{\prime},s^{\mu}_{-i})-u_{i}(s^{\mu})\leq
2\epsilon\quad\forall\rho_{i}^{\prime}\in\mathrm{P}_{i}.$

∎

## References

  * Fudenberg and Tirole [1991] Drew Fudenberg and Jean Tirole.  _Game Theory_.  MIT Press Books. The MIT Press, 1991.
  * Greenwald and Jafari [2003] Amy Greenwald and Amir Jafari.  A general class of no-regret learning algorithms and game-theoretic equilibria.  In _COLT_ , volume 3, pages 2–12, 2003.
  * Moulin and Vial [1978] H. Moulin and J.P. Vial.  Strategically zero-sum games: The class of games whose completely mixed equilibria cannot be improved upon.  _International Journal of Game Theory_ , 7(3):201–221, 1978.