# On the preservation of second integrals by Runge-Kutta methods

Benjamin K Tapley email<EMAIL_ADDRESS>Department of Mathematical
Sciences, The Norwegian University of Science and Technology, 7491 Trondheim,
Norway

###### Abstract

One can elucidate integrability properties of ordinary differential equations
(ODEs) by knowing the existence of second integrals (also known as weak
integrals or Darboux polynomials for polynomial ODEs). However, little is
known about how they are preserved, if at all, under numerical methods. Here,
we show that in general all Runge-Kutta methods will preserve all affine
second integrals but not (irreducible) quadratic second integrals. A number of
interesting corollaries are discussed, such as the preservation of certain
rational integrals by arbitrary Runge-Kutta methods. We also study the special
case of affine second integrals with constant cofactor and discuss the
preservation of affine higher integrals.

##  1 Introduction

Many ordinary differential equations (ODEs) possess structure that gives rise
to certain qualitative behaviors of the exact solution. For example, the
Hamiltonian function remains constant along the flow of Hamiltonian systems
[1] or the phase space volume is preserved along the flow of source-free ODEs
[2]. When solving such ODEs numerically, it is important, especially for long-
term simulations, that these features are preserved when possible [3]. The
preservation of first integrals is arguably one of the most important
properties for a numerical method to inherit and there exist many well
designed numerical methods for this purpose such as discrete gradient methods
[4]. Runge-Kutta methods, on the other hand, will always preserve all linear
first integrals [3] and they preserve all quadratic first integrals when their
stability matrix vanishes [5]. Runge-Kutta methods cannot preserve all cubic
(or higher) first integrals [3], however there exist particular Runge-Kutta
methods that can preserve a polynomial Hamiltonian [6]. The theory behind
preserving first integrals when solving ODEs numerically has been the center
of many influential studies [7, 8, 9, 10], however, many ODEs also possess
second integrals [11], which are also important for numerical methods to
preserve. Second integrals also go by the names weak integrals, special
integrals [12], algebraic invariant curves [13] or particular algebraic
solutions [14], for example, and will refer to polynomial second integrals of
polynomial vector fields as Darboux polynomials [15]. Second integrals are
important to study as they behave like first integrals on their zero level
sets. These level sets also divide phase space into regions that are
qualitatively different (e.g., bounded or unbounded) and represent barriers
for transport. Moreover, many first integrals, especially rational or time-
dependent first integrals, are constructed by taking certain products and/or
quotients of second integrals. Despite their importance, the preservation of
second integrals by numerical methods has not been addressed until only very
recently [16, 17]. In this paper we will discuss the second-integral-
preserving properties of Runge-Kutta methods.

Consider an autonomous ODE in $\mathbb{R}^{n}$

$\dot{\boldsymbol{x}}=\boldsymbol{f}(\boldsymbol{x}),$ (1)

then a second integral [11] of the ODE (1) is a function $p(\boldsymbol{x})$
that satisfies

$\dot{p}(\boldsymbol{x})=\boldsymbol{f}(\boldsymbol{x})^{T}\nabla
p(\boldsymbol{x})=c(\boldsymbol{x})p(\boldsymbol{x}),$ (2)

where the dot denotes $\frac{\mathrm{d}}{\mathrm{d}t}$ and $c(\boldsymbol{x})$
is called the cofactor of $p(\boldsymbol{x})$. If
$\boldsymbol{f}(\boldsymbol{x})$ is polynomial, then
$p(\boldsymbol{x})\in\mathbb{K}_{m}[\boldsymbol{x}]$ is referred to as a
Darboux polynomial, where $\mathbb{K}_{m}[\boldsymbol{x}]$ is the class of
polynomials over the field $\mathbb{K}$ of degree $m$ in the variables
$\boldsymbol{x}$. The discrete-time analogue of $p(\boldsymbol{x})$, is
referred to as a discrete second integral of a map
$\varphi_{h}:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ (or discrete Darboux
polynomial when $\varphi_{h}(\boldsymbol{x})$ is polynomial or rational),
which is a function $p(\boldsymbol{x})$ that satisfies

$p(\varphi_{h}(\boldsymbol{x}))=\tilde{c}(\boldsymbol{x})p(\boldsymbol{x}),$
(3)

where $\tilde{c}(\boldsymbol{x})$ is called the discrete cofactor of
$p(\boldsymbol{x})$. This is a discrete analogue of equation (2) and was
recently introduced in [16, 17]. In our context,
$\varphi_{h}(\boldsymbol{x})\approx\boldsymbol{x}(h)$ is a numerical method
designed to approximate the solution of the ODE (1) for some time step $h$.
The map $\varphi_{h}$ is said to preserve a Darboux polynomial
$p(\boldsymbol{x})$ if there exists a discrete cofactor
$\tilde{c}(\boldsymbol{x})$ that satisfies (3) and is independent of the
choice of basis for $p(\boldsymbol{x})$. The last point is important,
otherwise any polynomial $p(\boldsymbol{x})$ is a discrete Darboux polynomial
by letting
$\tilde{c}(\boldsymbol{x}):=p(\varphi_{h}(\boldsymbol{x}))/p(\boldsymbol{x})$.

In this paper we will mainly consider ODEs with one or more affine second
integrals of the form $p(\boldsymbol{x})=\boldsymbol{p}^{T}\boldsymbol{x}+r$
where $\boldsymbol{p}\in\mathbb{R}^{n}$ and $r\in\mathbb{R}$. Note that we can
take the constant $r=0$ without loss of generality as we will do throughout
the paper. In particular, we focus on the case where the map $\varphi_{h}$
comes about as a Runge-Kutta method applied to an ODE. That is, letting
$\varphi_{h}(\boldsymbol{x})$ denote one step of an $s$-stage Runge-Kutta
method applied to (1) with initial condition $\boldsymbol{x}$ and Butcher
tableau given by

$\begin{array}[]{c|c}\mathcal{C}&\mathcal{A}\\\
\hline\cr{}\hfil&\boldsymbol{b}^{T}\end{array}$ (4)

then this defines the method

$\displaystyle\boldsymbol{g}_{i}=$
$\displaystyle\boldsymbol{x}+h\sum_{j=1}^{s}a_{ij}\boldsymbol{f}(\boldsymbol{g}_{j}),\quad{\mathrm{for}~{}i=1,...,s}$
(5) $\displaystyle\varphi_{h}(\boldsymbol{x})=$
$\displaystyle\boldsymbol{x}+h\sum_{j=1}^{s}b_{j}\boldsymbol{f}(\boldsymbol{g}_{j}),$
(6)

where $\varphi_{h}(\boldsymbol{x})$ is the Runge-Kutta map with time-step $h$
applied to the initial point $\boldsymbol{x}$ and $\mathcal{A}=[a_{ij}]$. For
explicit Runge-Kutta maps the sum in equation (5) runs from 1 to $i-1$ as
$\mathcal{A}$ is strictly lower triangular.

The paper begins by proving that Runge-Kutta methods cannot preserve
irreducible quadratic second integrals. We then prove that Runge-Kutta methods
always preserve affine second integrals and discuss the implications of this
theorem including the fact that all Runge-Kutta methods preserve all rational
integrals that are affine in the numerator and denominator. We then focus on
the special case of affine second integrals with constant cofactor. The final
section is on the preservation of affine higher integrals.

##  2 The non-preservation of quadratic second integrals

######  Theorem 1.

Given an autonomous ODE with an irreducible quadratic Darboux polynomial
$p(\boldsymbol{x})\in\mathbb{K}_{2}[\boldsymbol{x}]$, then no Runge-Kutta
method can preserve $p(\boldsymbol{x})$.

###### Proof.

The proof begins in a similar way to that of [3, Thm. 2.2]. Consider a
quadratic second integral
$p(\boldsymbol{x})=\boldsymbol{x}^{T}Q\boldsymbol{x}$ for some symmetric
adjoint matrix $Q\in\mathbb{K}^{n\times n}$. Then

$p(\varphi_{h}(\boldsymbol{x}))=\boldsymbol{x}^{T}Q\boldsymbol{x}+h\sum_{i=1}^{s}b_{i}\boldsymbol{f}(\boldsymbol{g}_{i})^{T}Q\boldsymbol{x}+h\sum_{j=1}^{s}b_{j}\boldsymbol{x}^{T}Q\boldsymbol{f}(\boldsymbol{g}_{i})+h^{2}\sum_{i,j=1}^{s}b_{i}b_{j}\boldsymbol{f}(\boldsymbol{g}_{i})^{T}Q\boldsymbol{f}(\boldsymbol{g}_{j}).$
(7)

Inserting
$\boldsymbol{x}=\boldsymbol{g}_{i}-h\sum_{j=1}^{s}a_{ij}\boldsymbol{f}(\boldsymbol{g}_{j})$
in the two $O(h)$ terms yields

$p(\varphi_{h}(\boldsymbol{x}))=\boldsymbol{x}^{T}Q\boldsymbol{x}+2h\sum_{i=1}^{s}b_{i}\boldsymbol{g}_{i}^{T}Q\boldsymbol{f}(\boldsymbol{g}_{i})+h^{2}\sum_{i,j=1}^{s}(b_{i}b_{j}-b_{i}a_{ij}-b_{j}a_{ji})\boldsymbol{f}(\boldsymbol{g}_{i})^{T}Q\boldsymbol{f}(\boldsymbol{g}_{j}).$
(8)

Note that equation (2) implies
$\boldsymbol{g}_{i}^{T}Q\boldsymbol{f}(\boldsymbol{g}_{i})=c(\boldsymbol{g}_{i})p(\boldsymbol{g}_{i})$.
Computing $p(\boldsymbol{g}_{i})$ gives

$\displaystyle p(\boldsymbol{g}_{i})=$
$\displaystyle\boldsymbol{x}^{T}Q\boldsymbol{x}+h\sum_{i=1}^{s}a_{ij}\boldsymbol{f}(\boldsymbol{g}_{j})^{T}Q\boldsymbol{x}+h\sum_{j=1}^{s}a_{ik}\boldsymbol{x}^{T}Q\boldsymbol{f}(\boldsymbol{g}_{k})+h^{2}\sum_{j,k=1}^{s}a_{ij}a_{ik}\boldsymbol{f}(\boldsymbol{g}_{j})^{T}Q\boldsymbol{f}(\boldsymbol{g}_{k})$
(9) $\displaystyle=$ $\displaystyle
p(\boldsymbol{x})+h\sum_{j=1}^{s}a_{ij}c(\boldsymbol{g}_{j})p(\boldsymbol{g}_{j})+\overbrace{h^{2}\sum_{j,k=1}^{s}(a_{ij}a_{ik}-a_{ij}a_{jk}-a_{ik}a_{kj})\boldsymbol{f}(\boldsymbol{g}_{j})^{T}Q\boldsymbol{f}(\boldsymbol{g}_{k})}^{:=H_{i}},$
(10)

where we have again inserted the expression for $\boldsymbol{x}$ and used
equation (2) to arrive at (10). Note that this is now a linear system of $s$
equations for $p(\boldsymbol{g}_{j})$, which we can solve for. Denoting by
$\mathbf{P}:=(p(\boldsymbol{g}_{1}),...,p(\boldsymbol{g}_{s}))^{T}\in\mathbb{R}^{s}$
and $\mathbf{H}\in\mathbb{R}^{s}$ the vector whose $i$th component is defined
in equation (10) then

$\displaystyle\mathbf{P}$
$\displaystyle=\mathbb{1}p(\boldsymbol{x})+h\mathcal{A}D_{c}\mathbf{P}+\mathbf{H}$
(11)
$\displaystyle=(I-h\mathcal{A}D_{c})^{-1}\left(p(\boldsymbol{x})\mathbb{1}+\mathbf{H}\right),$
(12)

where
$D_{c}:=\mathrm{diag}([c(\boldsymbol{g}_{1}),...,c(\boldsymbol{g}_{s})])\in\mathbb{R}^{s\times
s}$ and $\mathbb{1}\in\mathbb{R}^{s}$ is the ones vector. Equation (8)
therefore reads

$\displaystyle p(\varphi_{h}(\boldsymbol{x}))=$ $\displaystyle
p(\boldsymbol{x})+h\boldsymbol{b}^{T}D_{c}\mathbf{P}+h^{2}\sum_{i,j=1}^{s}(b_{i}b_{j}-b_{i}a_{ij}-b_{j}a_{ji})\boldsymbol{f}(\boldsymbol{g}_{i})^{T}Q\boldsymbol{f}(\boldsymbol{g}_{j})$
(13) $\displaystyle=$ $\displaystyle
p(\boldsymbol{x})+h\boldsymbol{b}^{T}D_{c}(I-h\mathcal{A}D_{c})^{-1}\left(p(\boldsymbol{x})\mathbb{1}+\mathbf{H}\right)$
(14)
$\displaystyle\qquad+h^{2}\sum_{i,j=1}^{s}(b_{i}b_{j}-b_{i}a_{ij}-b_{j}a_{ji})\boldsymbol{f}(\boldsymbol{g}_{i})^{T}Q\boldsymbol{f}(\boldsymbol{g}_{j}).$
(15)

For the discrete cofactor
$\tilde{c}(\boldsymbol{x})=p(\varphi_{h}(\boldsymbol{x}))/p(\boldsymbol{x})$
to be independent of $p(\boldsymbol{x})$, we require that $p(\boldsymbol{x})$
divides equation (15). This is only true if the following conditions are
satisfied

$b_{i}b_{j}-b_{i}a_{ij}-b_{j}a_{ji}=0,\quad
a_{ij}a_{ik}-a_{ij}a_{jk}-a_{ik}a_{kj}=0,\quad\text{for}\quad i,j,k=1,...,s$
(16)

which is only true when $a_{i,j}=0$ and $b_{j}=0$ for $i,j=1,...,s$. Therefore
no Runge-Kutta method can, in general, preserve irreducible quadratic second
integrals. ∎

Note that above theorem only applies to irreducible second integrals as Runge-
Kutta methods can preserve quadratic second integrals that are the product of
affine second integrals. For example, when a quadratic second integral can be
expressed as
$p(\boldsymbol{x})=\boldsymbol{x}^{T}Q\boldsymbol{x}=(\boldsymbol{x}^{T}\boldsymbol{u})(\boldsymbol{v}^{T}\boldsymbol{x})$
where $Q:=\boldsymbol{u}\boldsymbol{v}^{T}$ then the preservation of
$p(\boldsymbol{x})$ is equivalent to the preservation of the two linear second
integrals $\boldsymbol{x}^{T}\boldsymbol{u}$ and
$\boldsymbol{x}^{T}\boldsymbol{v}$. Due to theorem 2 in the next section, this
is always true for Runge-Kutta methods.

Furthermore, in [16] many examples of Kahan’s method (which is a Runge-Kutta
method for quadratic vector fields [18]) preserving irreducible quadratic
Darboux polynomials are given. This suggests that despite the negative result
of theorem 1, there do exist cases where certain Runge-Kutta methods can
preserve certain quadratic second integrals. We will leave this for a future
study and instead focus now on the preservation of affine second integrals by
Runge-Kutta methods.

##  3 The preservation of affine second integrals

We begin with an example of a planar ODE and its discretisation by a second-
order Runge-Kutta method.

######  Example 1.

Consider the following ODE in two dimensions

$\dot{x}={x}^{2}+2\,xy+3\,{y}^{2},\quad\dot{y}=2\,y\left(2\,x+y\right),$ (17)

This ODE was studied in [19] and has the following three linear Darboux
polynomials

$p_{1}(\boldsymbol{x})=x+y,\quad p_{2}(\boldsymbol{x})=x-y,\quad
p_{3}(\boldsymbol{x})=y,$ (18)

that correspond to the cofactors

$c_{1}(\boldsymbol{x})=x+5y,\quad c_{2}(\boldsymbol{x})=x-y,\quad
c_{3}(\boldsymbol{x})=4x+2y.$ (19)

(a)

(b)

Figure 1: The phase portrait of the ODE (17) and the errors of the second
integrals $p_{i}(\boldsymbol{x})$ for initial conditions satisfying
$p_{i}(\boldsymbol{x}_{0})=0$. Note that
$p_{3}(\boldsymbol{x}_{n})-p_{3}(\boldsymbol{x}_{0})=0$ and therefore does not
show on the semi-log axis. The initial conditions are shown by black dots and
are located on the grid $(-10+i,-10+j)$ for $i,j=0,...,20$.

Consequently, the system also possesses the first integral
$p_{1}(\boldsymbol{x})p_{2}(\boldsymbol{x})^{3}p_{3}(\boldsymbol{x})^{-1}$.
The system is discretised using Ralston’s method with a time step of $h=0.001$
and the phase portrait on the square $[-10,10]^{2}$ is presented in figure
1(a). Here, the level sets $p_{i}(\boldsymbol{x})=0$ for $i=1,2,3$ are
represented by blue dashed lines. We see that numerical solutions starting on
one of these zero level sets remain on the level set. This is exemplified in
figure 1(b) which shows the errors
$p_{i}(\boldsymbol{x}_{n})-p_{i}(\boldsymbol{x}_{0})$ for the numerical
solutions starting from $p_{i}(\boldsymbol{x}_{0})=0$ and for $i=1,2,3$. Here,
we see that the errors are all within machine precision, implying that these
three second integrals behave like first integrals on their zero level sets as
the exact solution dictates. ∎

Second integrals are important as they divide phase space into sections with
qualitatively different behavior. We see that Ralston’s method has preserved
the Darboux polynomials $p_{i}(\boldsymbol{x})$, meaning that it produces a
qualitatively similar phase portrait. This fact is due to the following
theorem.

######  Theorem 2.

If an autonomous ODE $\dot{\boldsymbol{x}}=\boldsymbol{f}(\boldsymbol{x})$
possesses an affine second integral
$p(\boldsymbol{x})=\boldsymbol{p}^{T}\boldsymbol{x}$ for
$\boldsymbol{p}\in\mathbb{R}^{n}$ with cofactor $c(\boldsymbol{x})$ satisfying
$\boldsymbol{p}^{T}\boldsymbol{f}=c(\boldsymbol{x})\boldsymbol{p}^{T}\boldsymbol{x}$
then a Runge-Kutta map $\varphi_{h}$ of the ODE possesses the discrete second
integral $\boldsymbol{p}^{T}\boldsymbol{x}$ that satisfies
$\boldsymbol{p}^{T}\varphi_{h}(\boldsymbol{x})=\tilde{c}(\boldsymbol{x})\,\boldsymbol{p}^{T}\boldsymbol{x}$
where the discrete cofactor is given by

$\tilde{c}(\boldsymbol{x})=1+h\boldsymbol{b}^{T}D_{c}\,(I-h\mathcal{A}D_{c})^{-1}\mathbb{1},$
(20)

where
$D_{c}:=\mathrm{diag}([c(\boldsymbol{g}_{1}),...,c(\boldsymbol{g}_{s})])\in\mathbb{R}^{s\times
s}$ and $\mathbb{1}\in\mathbb{R}^{n}$ is the vector of ones.

###### Proof.

Let $\varphi_{h}$ denote the Runge-Kutta map defined by equations (4), (5) and
(6). Now let $G:=(\boldsymbol{g}_{1},...,\boldsymbol{g}_{s})^{T}$ and
$F:=(\boldsymbol{f}(\boldsymbol{g}_{i}),...,\boldsymbol{f}(\boldsymbol{g}_{s}))^{T}$
denote the $s\times n$ matrices whose $i$’th rows are $\boldsymbol{g}_{i}^{T}$
and $\boldsymbol{f}_{i}^{T}$, respectively. Then

$p(\varphi_{h}(\boldsymbol{x}))=\boldsymbol{p}^{T}\boldsymbol{x}+h\sum_{j=1}^{s}b_{j}\boldsymbol{p}^{T}\boldsymbol{f}(\boldsymbol{g}_{j})=\boldsymbol{p}^{T}\boldsymbol{x}+h\boldsymbol{b}^{T}F\boldsymbol{p}=\boldsymbol{p}^{T}\boldsymbol{x}+h\boldsymbol{b}^{T}D_{c}G\boldsymbol{p}.$
(21)

We have for $G\boldsymbol{p}$ the following

$G\boldsymbol{p}=\mathbb{1}_{s}\boldsymbol{p}^{T}\boldsymbol{x}+h\mathcal{A}F\boldsymbol{p}=\left(I-h\,\mathcal{A}D_{c}\right)^{-1}\mathbb{1}_{s}\boldsymbol{p}^{T}\boldsymbol{x}$
(22)

due to the fact that $F\boldsymbol{p}=D_{c}G\boldsymbol{p}$. Inserting (22)
into (21) and dividing by $\boldsymbol{p}^{T}\boldsymbol{x}$ we arrive at the
desired result

$\frac{\boldsymbol{p}^{T}\varphi_{h}(\boldsymbol{x})}{\boldsymbol{p}^{T}\boldsymbol{x}}=\tilde{c}=1+h\boldsymbol{b}^{T}D_{c}\,(I-h\mathcal{A}D_{c})^{-1}\mathbb{1}_{s}.$
(23)

∎

This is a generalisation of theorem 1 in [17], which shows that Kahan’s method
preserves all affine Darboux polynomials for quadratic vector fields. The
discrete cofactor $\tilde{c}(\boldsymbol{x})$ of theorem 2 depends only on the
Butcher table coefficients, the vector field $\boldsymbol{f}(\boldsymbol{x})$
and the (continuous) cofactor $c(\boldsymbol{x})$. Furthermore,
$\tilde{c}(\boldsymbol{x})$ is in general rational and implicitly defined due
the dependence of $D_{c}$ on the stage values $\boldsymbol{g}_{i}$. However,
when $\boldsymbol{f}(\boldsymbol{x})$ is polynomial then explicit Runge-Kutta
maps applied to $\boldsymbol{f}(\boldsymbol{x})$ yield polynomial maps. One
would therefore expect $\tilde{c}$ to be known explicitly and be polynomial.

######  Remark 1.

For all explicit Runge-Kutta methods applied to polynomial ODEs then
$\tilde{c}(\boldsymbol{x})$ is polynomial and can be written explicitly. This
can be shown by observing that the matrix $I-h\mathcal{A}D_{c}=I+L$ where
$L:=-h\mathcal{A}D_{c}$ is strictly lower triangular and therefore
$(I+L)^{-1}=I+\tilde{L}$ where $\tilde{L}$ is also strictly lower triangular.
Moreover, as $\det(I+L)=1$, its inverse is equal to its adjugate, that is
$(I+L)^{-1}=\mathrm{adj}(I+L)=I+\tilde{L}$ and therefore $\tilde{L}$ is
polynomial in the components of $L$. As the cofactor $c(\boldsymbol{x})$ is
polynomial for polynomial ODEs, it follows that $\tilde{c}(\boldsymbol{x})$ is
polynomial.

An interesting implication of theorem 2 is that if an ODE possesses two linear
second integrals with the same cofactor, then so does $\varphi_{h}$. This
leads to the following corollary about the preservation of rational integrals.

######  Corollary 1.

All Runge-Kutta methods preserve all rational first integrals of the form
$H(\boldsymbol{x})=Q(\boldsymbol{x})/R(\boldsymbol{x})$ for
$Q(\boldsymbol{x})$ and $R(\boldsymbol{x})$ affine.

###### Proof.

Any ODE with a first integral
$H(\boldsymbol{x})=Q(\boldsymbol{x})/R(\boldsymbol{x})$ can be written as the
following system

$\dot{\boldsymbol{x}}=\boldsymbol{f}(\boldsymbol{x})=S(\boldsymbol{x})\nabla\left(\frac{Q(\boldsymbol{x})}{R(\boldsymbol{x})}\right)$
(24)

for some skew-symmetric matrix $S(\boldsymbol{x})=-S(\boldsymbol{x})^{T}$.
Without loss of generality we can let
$Q(\boldsymbol{x})=\boldsymbol{q}^{T}\boldsymbol{x}$ and
$R(\boldsymbol{x})=\boldsymbol{r}^{T}\boldsymbol{x}$ for constant vectors
$\boldsymbol{q},\boldsymbol{r}\in\mathbb{R}^{n}$. We now show that these two
functions are second integrals of the ODE with the same cofactor by computing
their time derivatives

$\frac{\mathrm{d}}{\mathrm{d}t}Q(\boldsymbol{x})=\boldsymbol{q}^{T}\boldsymbol{f}=\frac{\boldsymbol{q}^{T}S(\boldsymbol{x})\left(\boldsymbol{q}^{T}\boldsymbol{x}\,\boldsymbol{r}-\boldsymbol{r}^{T}\boldsymbol{x}\,\boldsymbol{q}\right)}{{R(\boldsymbol{x})^{2}}}=\frac{\boldsymbol{q}^{T}S(\boldsymbol{x})\boldsymbol{r}}{{R(\boldsymbol{x})^{2}}}Q(\boldsymbol{x})$
(25)

due to skew-symmetry of $S(\boldsymbol{x})$. Similarly,

$\frac{\mathrm{d}}{\mathrm{d}t}R(\boldsymbol{x})=\frac{\boldsymbol{q}^{T}S(\boldsymbol{x})\boldsymbol{r}}{{R(\boldsymbol{x})^{2}}}R(\boldsymbol{x}),$
(26)

that is, $Q(\boldsymbol{x})$ and $R(\boldsymbol{x})$ are second integrals with
cofactor
$\frac{\boldsymbol{q}^{T}S(\boldsymbol{x})\boldsymbol{r}}{{R(\boldsymbol{x})^{2}}}$.
Due to theorem 2, all Runge-Kutta methods preserve these second integrals and
they both correspond to the same discrete cofactor. Therefore their quotient
is an integral of the Runge-Kutta map. ∎

Note that corollary 1 applies to general ODEs. The same statement can me made
for polynomial ODEs by scaling (24) by $R(x)^{a}$, $a\geq 2$ and
$a\in\mathbb{N}$. We now give two examples of explicit Runge-Kutta methods
preserving a rational integral. The first is of a simple Lotka-Volterra
system, the second is of a non-polynomial vector field.

######  Example 2 (A Lotka-Volterra system with a rational integral).

Consider the following 2D Lotka-Volterra system

$\dot{x}=x(x-y),\quad\dot{y}=y(x-y).$ (27)

Here $x$ and $y$ are clearly Darboux polynomials with cofactor $c=x-y$
implying that $H(\boldsymbol{x})=\frac{x}{y}$ is a first integral of the ODE.
Now consider the generic second order explicit Runge-Kutta map $\varphi_{h}$
defined by the Butcher Tableau

$\begin{array}[]{c|cc}0&0&0\\\ \theta&\theta&0\\\
\hline\cr{}\hfil&1-\frac{1}{2\theta}&\frac{1}{2\theta}\\\ \end{array}$ (28)

Note that setting $\theta=1/2,2/3,1$ yields the explicit midpoint, Ralston’s
and Heun’s methods, respectively. Then according to theorem 2 the discrete
cofactor is

$\tilde{c}(\boldsymbol{x})=1+{{\left(\,x_{{1}}-\,x_{{2}}\right)h}}+\left(x_{{1}}-x_{{2}}\right)^{2}{h}^{2}+{\frac{\theta}{2}\,\left(x_{{1}}-x_{{2}}\right)^{3}{h}^{3}}$
(29)

and is polynomial due to remark 1. Indeed, we can show that the Runge-Kutta
map satisfies

$\varphi_{h}(\boldsymbol{x})=\tilde{c}(\boldsymbol{x})\boldsymbol{x}$ (30)

meaning that $x$ and $y$ are discrete Darboux polynomials of $\varphi_{h}$
with cofactor $\tilde{c}$. This implies that $H(\boldsymbol{x})=\frac{x}{y}$
is a first integral of $\varphi_{h}$. ∎

######  Example 3 (A non-polynomial ODE with a rational integral).

Consider the following ODE

$\dot{x}=\frac{2\,x+\alpha}{\left(x-y\right)^{{\frac{3}{2}}}},\quad\dot{y}=\frac{2\,y+\alpha}{\left(x-y\right)^{{\frac{3}{2}}}}$
(31)

which has the following two second integrals

$p_{1}(\boldsymbol{x})=x+y+\alpha,\quad p_{2}(\boldsymbol{x})=x-y$ (32)

that both correspond to the cofactor

$c(\boldsymbol{x})=\frac{2}{\left(x-y\right)^{3/2}}.$ (33)

This means that $H(\boldsymbol{x})=(x+y+\alpha)/(x-y)$ is a first integral of
the ODE. Now apply to this ODE the generic second-order explicit Runge-Kutta
map $\varphi_{h}$ from example 2. Then $\varphi_{h}$ possesses the discrete
second integrals $p_{1}(\boldsymbol{x})$ and $p_{2}(\boldsymbol{x})$ with
cofactor

$\tilde{c}(\boldsymbol{x})=\frac{\left(h\sqrt{x-y}+a(\boldsymbol{x})\left(\left(\left(x-y\right)^{{\frac{3}{2}}}+2\,h\right)\theta-h\right)\right)}{a(\boldsymbol{x})\theta\left(x-y\right)^{{\frac{3}{2}}}}$
(34)

where

$a(\boldsymbol{x})=\sqrt{{\left(\left(x-y\right)^{{\frac{1}{2}}}+2\,h\theta\left(x-y\right)^{-{\frac{1}{2}}}\right)}},$
(35)

in agreement with theorem 2, hence $H(\boldsymbol{x})$ is an integral of
$\varphi_{h}$. We can also verify by direct computation that $\varphi_{h}$
satisfies $H(\varphi_{h}(\boldsymbol{x}))=H(\boldsymbol{x})$. ∎

If an ODE possesses a rational integral of the form prescribed by corollary 1,
then it is not necessarily a straightforward task to determine its form
without knowledge of the cofactor. However, due to the following theorem, one
can immediately determine the discrete cofactor corresponding to the numerator
and denominator of the rational integral by computing the Jacobian determinant
of a Runge-Kutta map applied to the ODE. From this one can infer the
(continuous) cofactor due to the fact that $\lim\limits_{h\rightarrow
0}((\tilde{c}(\boldsymbol{x})-1)/h)=c(\boldsymbol{x})$

######  Theorem 3.

If a Runge-Kutta map $\varphi_{h}$ possesses two or more affine second
integrals with the same cofactor $\tilde{c}(\boldsymbol{x})$, then
$\tilde{c}(\boldsymbol{x})$ divides the Jacobian determinant of the map
$\det\left(\frac{\partial\varphi_{h}(\boldsymbol{x})}{\partial\boldsymbol{x}}\right)$.

###### Proof.

Without loss of generality we can assume the two affine second integrals are
$p_{1}(\boldsymbol{x})=x_{1}$ and $p_{2}(\boldsymbol{x})=x_{2}$ that both
correspond to the cofactor $\tilde{c}(\boldsymbol{x})$. Letting
$\varphi_{h}:(x_{1},...,x_{n})\mapsto\varphi_{h}(x_{1},...,x_{n}):=(x_{1}^{\prime},...,x_{n}^{\prime})$
then the map can be written as

$\varphi_{h}\left(\begin{array}[]{c}x_{1}\\\ x_{2}\\\ x_{3}\\\ \vdots\\\
x_{n}\end{array}\right)=\left(\begin{array}[]{c}x_{1}\tilde{c}(\boldsymbol{x})\\\
x_{2}\tilde{c}(\boldsymbol{x})\\\ x_{3}^{\prime}\\\ \vdots\\\
x_{n}^{\prime}\end{array}\right).$ (36)

Taking the Jacobian of $\varphi_{h}(\boldsymbol{x})$ with respect to
$\boldsymbol{x}$ gives

${J_{\varphi_{h}}:=\frac{\partial\varphi_{h}(\boldsymbol{x})}{\partial\boldsymbol{x}}={\left(\begin{array}[]{cccc}x_{1}\frac{\partial\tilde{c}}{\partial
x_{1}}&\dots&x_{1}\frac{\partial\tilde{c}}{\partial x_{n}}\\\
x_{2}\frac{\partial\tilde{c}}{\partial
x_{1}}&\dots&x_{2}\frac{\partial\tilde{c}}{\partial x_{n}}\\\ \frac{\partial
x^{\prime}_{3}}{\partial x_{1}}&\dots&\frac{\partial x^{\prime}_{3}}{\partial
x_{n}}\\\ \vdots&&\vdots\\\ \frac{\partial x^{\prime}_{n}}{\partial
x_{1}}&\dots&\frac{\partial x^{\prime}_{n}}{\partial x_{n}}\\\
\end{array}\right)}+\left(\begin{array}[]{ccccc}\tilde{c}(\boldsymbol{x})&0&0&\dots&0\\\
0&\tilde{c}(\boldsymbol{x})&0&&\\\ 0&0&0&&\\\ \vdots&&&\ddots&\vdots\\\
0&&&\dots&0\\\ \end{array}\right)}.$ (37)

As we are taking the determinant of $J_{\varphi_{h}}$ we can perform
elementary row and column operations to manipulate this matrix. Denoting by
$\boldsymbol{r}_{i}$ and $\boldsymbol{c}_{i}$ the $i$th row and column of
$J_{\varphi_{h}}$ then letting
$\boldsymbol{c}_{1}\rightarrow\sum\boldsymbol{c}_{i}$, then
$\boldsymbol{r}_{2}\rightarrow\boldsymbol{r}_{2}-\boldsymbol{r}_{1}$, we get
for row two $\boldsymbol{r}_{2}=(0,\tilde{c}(\boldsymbol{x}),0,...,0)$. Then
by Laplace expansion across this row, we see that
$\det(J_{\varphi_{h}})=\tilde{c}(\boldsymbol{x})\sum M_{2,j}$, where $M_{2,j}$
are the determinantal minors across the second row. ∎

We also remark that theorem 3 implies that if a Runge-Kutta map has $m$ second
integrals corresponding to the cofactor $\tilde{c}(\boldsymbol{x})$, then
$\tilde{c}(\boldsymbol{x})^{m-1}$ will divide the Jacobian. Furthermore, for
the proof to hold, we only need the map to be locally invertible, such that
$\det(J_{\varphi_{h}})\neq 0$ and therefore should hold for other affinely
equivalent maps that possess a rational integral.

Using theorem 3 one can easily detect the existence of a rational integral of
the form given in corollary 1 by factoring the Jacobian determinant of any
Runge-Kutta map applied to an ODE. This immediately provides the form of the
cofactor, if it exists as given by the following example.

######  Example 4 (A rational ODE with an unknown rational integral).

Consider the following ODE

$\displaystyle\dot{x}=$ $\displaystyle x-y+z,$ (38) $\displaystyle\dot{y}=$
$\displaystyle{\frac{2\,{y}^{2}-xz-{z}^{2}}{y+z}},$ (39)
$\displaystyle\dot{z}=$
$\displaystyle{\frac{{z}^{2}+z\left(x+y\right)-{y}^{2}}{y+z}}.$ (40)

It is not immediately obvious if this ODE possesses a rational integral.
However, we apply to this ODE the forward Euler method, denoted by
$\varphi_{h}(\boldsymbol{x})$, and take its Jacobian determinant. This yields

$\det\left(J_{\varphi_{h}}\right)=\frac{K_{1}K_{2}}{D^{2}}$ (41)

where

$\displaystyle K_{1}=$ $\displaystyle yh+y+z,$ (42) $\displaystyle K_{2}=$
$\displaystyle\left(x+2\,y+z\right){h}^{2}+\left(x+4\,y+3\,z\right)h+y+z,$
(43) $\displaystyle D=$ $\displaystyle y+z.$ (44)

Using $K_{1}/D$ as a discrete cofactor and setting
$p(\boldsymbol{x})=\alpha_{1}x+\alpha_{2}y+\alpha_{3}z+\alpha_{4}$ we solve

$p(\varphi_{h}(\boldsymbol{x}))=\frac{K_{1}}{D}p(\boldsymbol{x})$ (45)

for $\alpha_{i}$, which gives the following two discrete second integral
solutions

$p_{1}(\boldsymbol{x})=y+z,\quad\text{and}\quad p_{2}(\boldsymbol{x})=x+y.$
(46)

This implies that $H(\boldsymbol{x})=(y+z)/(x+y)$ is a first integral of
$\varphi_{h}(\boldsymbol{x})$ and therefore also of the original ODE.
Computing

$c(\boldsymbol{x})=\lim\limits_{h\rightarrow
0}\left(\frac{\tilde{c}(\boldsymbol{x})-1}{h}\right)=\frac{y}{y+z},$ (47)

yields the form of the rational (continuous) cofactor. ∎

We remark that the above method for finding rational affine integrals involves
only solving linear systems and hence is fast and scalable to ODE in higher
dimensions with free parameters.

###  3.1 Constant cofactor case

We now restrict our discussion to the special case of affine Darboux
polynomials with constant cofactor

$\dot{p}(\boldsymbol{x})=\lambda p(\boldsymbol{x}),$ (48)

where $\lambda$ is constant. If such a Darboux polynomial exists then we can
solve for $p(\boldsymbol{x})$

$p(\boldsymbol{x})=Ke^{\lambda t},$ (49)

for some arbitrary constant $K$, and therefore

$H=p(\boldsymbol{x})e^{-\lambda t}$ (50)

is a time-dependent integral of the ODE (1). This has a natural discrete-time
analogue as follows.

######  Proposition 1.

Consider a map $\varphi_{h}$ with a discrete Darboux polynomial satisfying

$p(\varphi_{h}(\boldsymbol{x}))=\tilde{c}p(\boldsymbol{x}),$

where $\tilde{c}$ is constant. This implies that

$p(\boldsymbol{x})=K\tilde{c}^{k},$ (51)

where $k$ is the iteration index and $K$ is an arbitrary constant, moreover
the function

$\tilde{H}=\tilde{c}^{-k}p(\boldsymbol{x}),$ (52)

is an integral of the map $\varphi_{h}^{\circ k}$.

We now consider the case where the map $\varphi_{h}$ comes about as a Runge-
Kutta method applied to an ODE that possesses one or more Darboux polynomials
of the form (48). We begin with an example of a Lotka-Volterra system with a
time-dependent integral and its discretisation under an arbitrary Runge-Kutta
method.

######  Example 5 (A Lotka-Volterra system with one time-dependent integral).

Consider the following Lotka-Volterra system

$\frac{\mathrm{d}}{\mathrm{d}t}\left(\begin{array}[]{c}x_{1}\\\ x_{2}\\\
x_{3}\\\
\end{array}\right)=\left(\begin{array}[]{c}x_{{1}}\left(a_{{1}}x_{{2}}-a_{{2}}x_{{3}}+b\right)\\\
x_{{2}}\left(-a_{{1}}x_{{1}}+a_{{3}}x_{{3}}+b\right)\\\
x_{{3}}\left(a_{{2}}x_{{1}}-a_{{3}}x_{{2}}+b\right)\\\ \end{array}\right).$
(53)

Then $p_{1}(x)=x_{1}+x_{2}+x_{3}$ is a Darboux polynomial corresponding to the
cofactor $c_{1}=b$, hence

$H=p_{1}(x)e^{-bt}$ (54)

is a time-dependent integral of the ODE (53). Now consider discretisation of
the above ODE by a the forward Euler method $\varphi_{h}$. Then, $p_{1}(x)$ is
a discrete Darboux polynomial of $\varphi_{h}$ with the cofactor
$\tilde{c}_{1}=1+bh$, hence

$\tilde{H}=\tilde{c}_{1}^{k}p_{1}(x),$ (55)

where $k$ is the iteration index, is an integral of the map
$\varphi_{h}^{\circ k}$. ∎

Note that $\tilde{c}_{1}$ in the above example also corresponds to the
stability function $R(bh)$ of the forward Euler method. This is no coincidence
and is due to the following corollary of theorem 2.

######  Corollary 2.

If an autonomous polynomial ODE
$\dot{\boldsymbol{x}}=\boldsymbol{f}(\boldsymbol{x})$ possesses an affine
Darboux polynomial $p(\boldsymbol{x})$ with a constant cofactor $\lambda$,
then a Runge-Kutta map $\varphi_{h}$ with stability function $R(z)$ satisfies

$p(\varphi_{h}(\boldsymbol{x}))=R(\lambda h)p(\boldsymbol{x}),$ (56)

where

$R(z)=1+z\boldsymbol{b}^{T}(I-z\mathcal{A})^{-1}\mathbb{1}$ (57)

is the (constant) discrete cofactor, $\mathbb{1}$ is the ones vector and $I$
is the $s\times s$ identity matrix.

###### Proof.

This can be seen by setting $c(\boldsymbol{x})=\lambda$ in theorem 2.
Alternatively, we note that as Runge-Kutta methods are affinely equivariant
[20], we can take $p(x)=x_{1}$ without loss of generality, then equation (48)
is identical to Dahlquist’s famous test ODE [21] and $R(z)$ is identical to
the the stability function in the context of A-stability of one-step methods.
The rest follows from e.g., [21, Proposition 3.1]. See appendix A for details.
∎

The fact that the cofactor is the stability function is intuitive. Equation
(2) implies that $p(\boldsymbol{x})=A_{0}\exp(\int
c(\boldsymbol{x})\mathrm{d}t)$ for some integration constant $A_{0}$. Letting
$c(\boldsymbol{x})=\lambda$ and $p(\boldsymbol{x})=x_{1}$ implies that
$x_{1}=A_{0}\exp(\lambda t)$ as expected for Dalquist’s linear test equation,
while an order-$p$ Runge-Kutta map yields $\varphi_{h}(x_{1})=R(\lambda
h)x_{1}$, where $R(z)\approx\exp(z)$ is an order-$p$ Padé approximation to the
exponential function.

This result implies that Runge-Kutta methods can generally preserve integrals
that are the product of Darboux polynomials. This can be seen by considering
an ODE that has the first integral
$H(\boldsymbol{x})=p_{1}(\boldsymbol{x})^{\alpha_{1}}p_{2}(\boldsymbol{x})^{\alpha_{2}}$
where the $p_{i}(\boldsymbol{x})$ are Darboux polynomials with constant
cofactors $\lambda_{i}$. The fact that $H(\boldsymbol{x})$ is a first integral
implies that $\alpha_{1}\lambda_{1}+\alpha_{2}\lambda_{2}=0$. Then
$H(\boldsymbol{x})$ is only a first integral of the Runge-Kutta map
$\varphi_{h}$ if and only if
$R(\lambda_{1})^{\alpha_{1}}R(\lambda_{2})^{\alpha_{2}}=1$, which is only true
only when $R(z)=\exp(z)$, which is impossible. The exception of course, is
when $\alpha_{1}=-\alpha_{2}=\pm 1$ and $\lambda_{1}=\lambda_{2}$. This
results in the integral
$H(\boldsymbol{x})=p_{1}(\boldsymbol{x})/p_{2}(\boldsymbol{x})$ being
preserved by all Runge-Kutta maps, as shown in corollary 1. However, despite
this result, there are some instances where Runge-Kutta methods can preserve a
modified integral as we will now discuss.

The existence of Darboux polynomials with constant cofactors implies the
existence of time-dependent integrals of an ODE and therefore iteration-index-
dependent integrals of a Runge-Kutta map in the discrete-time case. One can
eliminate this time dependence (and iteration index dependence) by taking
quotients. This leads to the preservation of some non-rational modified
integrals by Runge-Kutta methods by the following corollary.

######  Corollary 3.

Given an ODE with a (non-rational) first integral given by

$H(\boldsymbol{x})=\frac{p_{1}(\boldsymbol{x})^{\sigma}}{p_{2}(\boldsymbol{x})},$
(58)

where $p_{1}(\boldsymbol{x})$ and $p_{2}(\boldsymbol{x})$ are affine Darboux
polynomials with constant cofactors $c_{1}$ and $c_{2}$ and

$\sigma=\frac{c_{2}}{c_{1}},$ (59)

then any Runge-Kutta map $\varphi_{h}$ with stability function $R(z)$
preserves the modified integral

$\tilde{H}(\boldsymbol{x})=\frac{p_{1}(\boldsymbol{x})^{\tilde{\sigma}}}{p_{2}(\boldsymbol{x})}$
(60)

with the modified exponent

$\tilde{\sigma}=\frac{\ln(R(hc_{2}))}{\ln(R(hc_{1}))}.$ (61)

###### Proof.

According to corollary 2, $\varphi_{h}$ preserves the Darboux polynomials
$p_{1}(\boldsymbol{x})$ and $p_{2}(\boldsymbol{x})$ with the modified
cofactors $R(hc_{1})$ and $R(hc_{2})$. It follows that
$\tilde{H}(\boldsymbol{x})$ is also a Darboux polynomial of $\varphi_{h}$ with
cofactor 1, and hence is a first integral of $\varphi_{h}$. ∎

This is demonstrated by the following example.

######  Example 6 (An ODE with a non-rational integral).

Consider the following ODE in three dimensions

$\frac{\mathrm{d}}{\mathrm{d}t}\left(\begin{array}[]{c}x_{1}\\\ x_{2}\\\
x_{3}\\\
\end{array}\right)=\left(\begin{array}[]{c}x_{{1}}x_{{2}}+x_{{2}}x_{{3}}+x_{{1}}\\\
-x_{{1}}x_{{2}}-x_{{2}}x_{{3}}+x_{{2}}\\\
\left(\sigma-1\right)\left(x_{{1}}+x_{{2}}\right)+\sigma\,x_{{3}}\\\
\end{array}\right),$ (62)

where $\sigma\in\mathbb{R}$. Discretising the above ODE by the forward Euler
method $\varphi_{h}$ with stability function $R(z)=1+z$, we find the following
affine Darboux polynomials and their corresponding discrete constant cofactors

$\displaystyle p_{1}(\boldsymbol{x})=x_{1}+x_{2},\quad\tilde{c}_{1}=1+h,$ (63)
$\displaystyle
p_{2}(\boldsymbol{x})=x_{1}+x_{2}+x_{3},\quad\tilde{c}_{2}=1+\sigma h.$ (64)

Therefore

$\tilde{H}=\frac{(x_{1}+x_{2})^{{\tilde{\sigma}}}}{x_{1}+x_{2}+x_{3}}$ (65)

defines an $h$-dependent integral of $\varphi_{h}$ with

$\tilde{\sigma}=\frac{\ln\left(1+\sigma h\right)}{\ln\left(1+h\right)}.$ (66)

We can take the continuum limit $\lim\limits_{h\rightarrow
0}\left(\tilde{\sigma}\right)=\sigma$, which implies that
$\lim\limits_{h\rightarrow 0}\tilde{H}(x)=H(x)$, where

$H(x)=\frac{(x_{1}+x_{2})^{\sigma}}{x_{1}+x_{2}+x_{3}}.$ (67)

Indeed, the irrational integral $H(x)$ is preserved by the flow of the
original ODE. ∎

##  4 ODEs with affine higher integrals

The notion of first, second, third111A third integral is a function
$K(\boldsymbol{x})$ that is preserved on a particular level set of a first
integral $H(\boldsymbol{x})$, e.g.,
$\dot{K}(\boldsymbol{x})=c(\boldsymbol{x})H(\boldsymbol{x})$ and higher
integrals can be generalized by considering solutions of the linear system for
$\underline{{\boldsymbol{p}}}(\boldsymbol{x})=(p_{1}(\boldsymbol{x}),p_{2}(\boldsymbol{x}),...,p_{m}(\boldsymbol{x}))^{T}\in\mathbb{R}^{m}$
satisfying [11]

$\dot{\underline{\boldsymbol{p}}}(\boldsymbol{x})=L(\boldsymbol{x})\underline{\boldsymbol{p}}(\boldsymbol{x}),$
(68)

where $L(\boldsymbol{x})\in\mathbb{R}^{n\times n}$. If $L_{ij}=0$ for
$j=1,...,n$ then $p_{i}(\boldsymbol{x})$ is a first integral and is preserved
from arbitrary initial conditions. If $L_{ii}\neq 0$ and $L_{ij}=0$ for $j\neq
i$ then $p_{i}(\boldsymbol{x})$ is a Darboux polynomial with cofactor $L_{ii}$
and is preserved for initial conditions that begin on its zero level sets. If,
for example, $L_{ik}\neq 0$ and $L_{ij}=0$ for $j\neq i,k$ then
$p_{i}(\boldsymbol{x})$ is preserved for initial conditions starting on the
zero level sets of $p_{k}(\boldsymbol{x})$. Here, if $p_{k}(\boldsymbol{x})$
is a first integral, then $p_{i}(\boldsymbol{x})$ is called a third integral.
In this sense, one can define the notion of higher integrals for ODE systems.
In this section we consider Runge-Kutta methods applied to ODEs with a sub-
system of affine higher integrals with constant coefficient matrix $L$.

######  Corollary 4.

Consider an autonomous ODE in $n$ dimensions that possesses a system of $m\leq
n$ linearly-independent affine polynomials
$\underline{{\boldsymbol{p}}}(\boldsymbol{x})=(p_{1}(\boldsymbol{x}),p_{2}(\boldsymbol{x}),...,p_{m}(\boldsymbol{x}))^{T}\in\mathbb{R}^{m}$
satisfying

$\dot{\underline{\boldsymbol{p}}}(\boldsymbol{x})=L\underline{\boldsymbol{p}}(\boldsymbol{x}),$
(69)

where $L$ is an $m\times m$ matrix that is independent of $\boldsymbol{x}$.
Then a Runge-Kutta map $\varphi_{h}$ satisfies

$\underline{\boldsymbol{p}}(\varphi_{h}(\boldsymbol{x}))=R(hL)\underline{\boldsymbol{p}}(\boldsymbol{x}),$
(70)

where $R(z)$ is the stability function of $\varphi_{h}$.

###### Proof.

Let $\underline{\boldsymbol{p}}(\boldsymbol{x})=D\boldsymbol{x}$, where $D$ is
an $m\times n$ matrix of rank $m\leq n$. The numerical solution of
$\varphi_{h}$ applied to a linear ODE in $n$ dimensions
$\dot{\underline{\boldsymbol{p}}}=L\underline{\boldsymbol{p}}$ is given by
(see for example, [3, p. 194])

$\varphi_{h}(\underline{\boldsymbol{p}})=R(hL)\underline{\boldsymbol{p}}.$
(71)

As Runge-Kutta methods commute with linear transformations we have
$\varphi_{h}(\underline{\boldsymbol{p}}(\boldsymbol{x}))=\underline{\boldsymbol{p}}(\varphi_{h}(\boldsymbol{x}))$,
which yields the desired result. ∎

We remark that an ODE in $n$ dimensions cannot possess more than $n$
functionally independent affine Darboux polynomials with constant cofactor. We
demonstrate corollary 4 in an example where $L$ is given in Jordan form. But
first we will briefly introduce a particular class of Runge-Kutta methods
called diagonal Padé Runge-Kutta methods.

######  Definition 1.

A diagonal Padé Runge-Kutta map is an $s$-stage Runge-Kutta map $\varphi_{h}$
whose stability function $R(z)=\frac{P(z)}{Q(z)}$ has equal degree in the
numerator and denominator, that is $\deg(P(z))=\deg(Q(z))$. Such a stability
function is an order $s$ approximation to $e^{z}$ with numerator and
denominator given by

$\displaystyle P(z)$ $\displaystyle=\sum_{i=0}^{s}a_{i}z^{i}$ (72)
$\displaystyle Q(z)$ $\displaystyle=\sum_{i=0}^{s}(-1)^{i}a_{i}z^{i}=P(-z),$
(73)

where $a_{i}$ are the constants

$a_{i}=\frac{s!(2s-i)!}{(2s)!i!(s-i)!}.$ (74)

See, for example, [22] and references therein. Such a map has a stability
function that satisfies $R(-z)R(z)=1$.

As an example of some well known Diagonal Padé Runge-Kutta methods, consider
the Runge-Kutta map $\varphi_{h}(\boldsymbol{x})=\boldsymbol{x}^{\prime}$
defined by

$\frac{\boldsymbol{x}^{\prime}-\boldsymbol{x}}{h}=(1-2\theta)\boldsymbol{f}(\frac{\boldsymbol{x}+\boldsymbol{x}^{\prime}}{2})+\theta\boldsymbol{f}(\boldsymbol{x}^{\prime})+\theta\boldsymbol{f}(\boldsymbol{x})$
(75)

then its stability function is the diagonal Padé approximation

$R(\lambda h)=\frac{1+\frac{1}{2}\lambda h}{1-\frac{1}{2}\lambda h}.$ (76)

Note that the three cases $\theta=0$, $\theta=\frac{1}{2}$ and
$\theta=-\frac{1}{2}$ respectively correspond to the midpoint rule,
trapezoidal rule and Kahan’s method, the latter being when
$\boldsymbol{f}(\boldsymbol{x})$ is quadratic [18].

######  Example 7 (A 3D ODE with a 2D linear subsystem in Jordan form).

Consider the following quadratic ODE in 3 dimensions that was considered in
example (6)

$\frac{\mathrm{d}}{\mathrm{d}t}\left(\begin{array}[]{c}x_{1}\\\ x_{2}\\\
x_{3}\\\
\end{array}\right)=\left(\begin{array}[]{c}x_{{1}}x_{{2}}+x_{{2}}x_{{3}}+\sigma
x_{{2}}\\\ -x_{{1}}x_{{2}}-x_{{2}}x_{{3}}+\sigma x_{{1}}\\\
\left(x_{{1}}+x_{{2}}\right)+\sigma\,x_{{3}}\\\ \end{array}\right).$ (77)

Setting
$\underline{\boldsymbol{p}}(\boldsymbol{x})=(x_{1}+x_{2}+x_{3},x_{1}+x_{2})^{T}$
then the above ODE has the following linear subsystem

$\dot{\underline{\boldsymbol{p}}}(\boldsymbol{x})=\left(\begin{array}[]{cc}\sigma&1\\\
0&\sigma\\\
\end{array}\right)\underline{\boldsymbol{p}}(\boldsymbol{x}):=L\underline{\boldsymbol{p}}(\boldsymbol{x}).$
(78)

If we apply Kahan’s method then we get

$\displaystyle\underline{\boldsymbol{p}}(\varphi_{h}(\boldsymbol{x}))=$
$\displaystyle\left(\begin{array}[]{cc}{\frac{-h\sigma-2}{h\sigma-2}}&4\,{\frac{h}{\left(h\sigma-2\right)^{2}}}\\\
\vskip 6.0pt plus 2.0pt minus 2.0pt\cr
0&{\frac{-h\sigma-2}{h\sigma-2}}\end{array}\right)\underline{\boldsymbol{p}}(\boldsymbol{x})$
(81) $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cc}R(h\sigma)&R^{\prime}(h\sigma)\\\
0&R(h\sigma)\\\
\end{array}\right)\underline{\boldsymbol{p}}(\boldsymbol{x})=R(hL)\underline{\boldsymbol{p}}(\boldsymbol{x}),$
(84)

which is the form prescribed in corollary 4.

∎

It is known that when the Kahan map is applied to a Hamiltonian ODE that it
preserves a modified Hamiltonian [18]. Here we show that there exists some
cases where the Kahan map (as well as other diagonal Padé Runge-Kutta maps)
preserves a quadratic first integral exactly.

######  Example 8 (An ODE with a quadratic integral that is preserved
exactly).

Consider the quadratic ODE in three dimensions

$\frac{\mathrm{d}}{\mathrm{d}t}\left(\begin{array}[]{c}x_{1}\\\ x_{2}\\\
x_{3}\\\
\end{array}\right)=\left(\begin{array}[]{c}2\,x_{{2}}+x_{{1}}+x_{{3}}+x_{{3}}\left(x_{{1}}+x_{{2}}\right)+{x_{{1}}}^{2}\\\
-x_{{1}}-x_{{2}}-x_{{3}}\left(x_{{1}}+x_{{2}}\right)-{x_{{1}}}^{2}\\\
x_{{3}}\left(x_{{1}}+x_{{2}}\right)+{x_{{1}}}^{2}\end{array}\right)$ (85)

Setting
$\underline{\boldsymbol{p}}(\boldsymbol{x})=(x_{1}+x_{2},x_{2}+x_{3})^{T}$
then this satisfies

$\dot{\underline{\boldsymbol{p}}}(\boldsymbol{x})=\left(\begin{array}[]{cc}0&1\\\
-1&0\\\ \end{array}\right)\underline{\boldsymbol{p}}(\boldsymbol{x})$ (86)

If we apply the Kahan discretisation to the above ODE we find that
$p_{1}(\boldsymbol{x})=\|\underline{\boldsymbol{p}}(\boldsymbol{x})\|^{2}=(x_{1}+x_{2})^{2}+(x_{2}+x_{3})^{2}$
is a quadratic Darboux polynomial of $\varphi_{h}$ with cofactor
$\tilde{c}=1$. Therefore
$H(\boldsymbol{x})=\|\underline{\boldsymbol{p}}(\boldsymbol{x})\|^{2}$ is an
integral of the map, i.e., $H(\varphi_{h}(\boldsymbol{x}))=H(\boldsymbol{x})$.
Moreover, we find that any diagonal Padé Runge-Kutta map preserves this
integral. This is due to the following corollary. ∎

######  Corollary 5.

Given an ODE $\dot{\boldsymbol{x}}=\boldsymbol{f}(\boldsymbol{x})$, if there
exists $m$ linear polynomials
$\underline{\boldsymbol{p}}(x)=(p_{1}(x),p_{2}(x),...,p_{m}(x))^{T}$ that
satisfy

$\dot{\underline{\boldsymbol{p}}}=AS\underline{\boldsymbol{p}}(\boldsymbol{x}),\quad\mathrm{where}\quad
A=-A^{T}\quad\mathrm{and}\quad S=S^{T}$ (87)

then
$H(x)=\underline{\boldsymbol{p}}(\boldsymbol{x})^{T}S\underline{\boldsymbol{p}}(\boldsymbol{x})$
is a first integral of the ODE and any diagonal Padé Runge-Kutta map (e.g.,
one of the form (75)) preserves this integral exactly.

###### Proof.

Using the fact that $R(z)=P(z)P(-z)^{-1}=P(-z)^{-1}P(z)$, $P(hAS)^{T}=P(-hSA)$
and $P(hSA)S=S\,P(hAS)$ then it’s straight forward to show that $H$ is
preserved under $\varphi_{h}$

$\displaystyle H(\varphi_{h}(\boldsymbol{x}))=$
$\displaystyle\underline{\boldsymbol{p}}(\boldsymbol{x})^{T}R(hAS)^{T}SR(hAS)\underline{\boldsymbol{p}}(\boldsymbol{x})$
(88) $\displaystyle=$
$\displaystyle\underline{\boldsymbol{p}}(\boldsymbol{x})^{T}P(-hAS)^{-T}P(hAS)^{T}SP(hAS)P(-hAS)^{-1}\underline{\boldsymbol{p}}(\boldsymbol{x})$
(89) $\displaystyle=$
$\displaystyle\underline{\boldsymbol{p}}(\boldsymbol{x})^{T}P(hSA)^{-1}P(-hSA)SP(hAS)P(-hAS)^{-1}\underline{\boldsymbol{p}}(\boldsymbol{x})$
(90) $\displaystyle=$
$\displaystyle\underline{\boldsymbol{p}}(\boldsymbol{x})^{T}SP(hAS)^{-1}P(-hAS)P(-hAS)^{-1}P(hAS)\underline{\boldsymbol{p}}(\boldsymbol{x})$
(91) $\displaystyle=$
$\displaystyle\underline{\boldsymbol{p}}(\boldsymbol{x})^{T}S\underline{\boldsymbol{p}}(\boldsymbol{x})$
(92) $\displaystyle=$ $\displaystyle H(\boldsymbol{x})$ (93)

∎

We can therefore make the following statement about linear ODEs with quadratic
integrals.

######  Corollary 6.

For all linear ODEs with quadratic first integrals, all diagonal Padé Runge-
Kutta maps preserve the integral exactly.

###### Proof.

A linear ODE with a quadratic integral
$H=\frac{1}{2}\boldsymbol{x}^{T}S\boldsymbol{x}$ can be written in the form

$\dot{\boldsymbol{x}}=A\nabla H=AS\boldsymbol{x}$ (94)

for $A=-A^{T}$ and $S=S^{T}$, which is in the form of (87) and according to
corollary 5, all diagonal Padé Runge-Kutta maps preserve the integral $H$ ∎

This is a slight generalisation of proposition 5 in [6].

Corollary 5 also applies to ODEs with multiple quadratic integrals as shown in
the following example.

######  Example 9 (A 5D system with 2 quadratic integrals preserved exactly).

Consider the following ODE

$\frac{\mathrm{d}}{\mathrm{d}t}\left(\begin{array}[]{c}x_{1}\\\ x_{2}\\\
x_{3}\\\ x_{4}\\\ x_{5}\\\
\end{array}\right)=\left(\begin{array}[]{c}{x_{{3}}}^{2}+2\,x_{{3}}x_{{4}}+2\,x_{{3}}x_{{5}}+2\,{x_{{4}}}^{2}+2\,x_{{4}}x_{{5}}+2\,{x_{{5}}}^{2}-10\,x_{{1}}-17\,x_{{2}}-7\,x_{{3}}-5\,x_{{4}}+4\,x_{{5}}\par\\\
-{x_{{3}}}^{2}-2\,x_{{3}}x_{{4}}-2\,x_{{3}}x_{{5}}-2\,{x_{{4}}}^{2}-2\,x_{{4}}x_{{5}}-2\,{x_{{5}}}^{2}+6\,x_{{1}}+11\,x_{{2}}+5\,x_{{3}}+5\,x_{{4}}-3\,x_{{5}}\par\\\
{x_{{3}}}^{2}+2\,x_{{3}}x_{{4}}+2\,x_{{3}}x_{{5}}+2\,{x_{{4}}}^{2}+2\,x_{{4}}x_{{5}}+2\,{x_{{5}}}^{2}-x_{{2}}-x_{{3}}-4\,x_{{4}}+2\,x_{{5}}\par\\\
-x_{{1}}-2\,x_{{2}}-x_{{3}}+2\,x_{{4}}\par\\\
x_{{1}}+3\,x_{{2}}+2\,x_{{3}}+2\,x_{{4}}-2\,x_{{5}}\par\\\ \end{array}\right)$
(95)

Let $\lambda=2+i$, where $i^{2}=-1$ and discretise the ODE using Kahan’s
method, which is a diagonal Padé Runge-Kutta map. We find the following affine
discrete Darboux polynomials with constant cofactor.

$\begin{array}[]{c|c|c}j&p_{j}(x)&\tilde{c}_{j}\\\ \hline\cr
1&-ix_{{4}}+x_{{1}}+2\,x_{{2}}+x_{{3}}&R(\lambda h)\\\
2&\left(106+52\,i\right)x_{{2}}+41\,x_{{3}}-\left(52-24\,i\right)x_{{5}}+\left(65+52\,i\right)x_{{1}}+\left(12-15\,i\right)x_{{4}}&R(-\lambda
h)\\\ 3&ix_{{4}}+x_{{1}}+2\,x_{{2}}+x_{{3}}&R(\lambda^{*}h)\\\
4&\left(106-52\,i\right)x_{{2}}+41\,x_{{3}}-\left(52+24\,i\right)x_{{5}}+\left(65-52\,i\right)x_{{1}}+\left(12+15\,i\right)x_{{4}}&R(-\lambda^{*}h)\\\
\end{array}$ (96)

From the above four linear Darboux polynomials, we can construct the following
complex quadratic integrals

$K_{1}=p_{1}(\boldsymbol{x})p_{2}(\boldsymbol{x}),\quad
K_{2}=p_{3}(\boldsymbol{x})p_{4}(\boldsymbol{x}),$ (97)

however if we observe that $K_{1}=K_{2}^{*}$ then the following two
independent real integrals can be constructed

$H_{1}=K_{1}+K_{2},\quad H_{2}=i(K_{1}-K_{2}).$ (98)

As these two integrals are independent of $h$, it follows that the ODE also
possesses these integrals. ∎

###  4.1 Detecting the existence of affine higher integrals

Many of our results in this section so far rely on the fact that there exist a
change of coordinates $\underline{\boldsymbol{p}}(x)=Dx$ that allow us to find
a linear ODE subsystem of higher integrals of the form
$\dot{\underline{\boldsymbol{p}}}=L\underline{\boldsymbol{p}}$. A logical
question to ask is what is the most general class of ODEs where this sub-
system exists and how to calculate this transformation. This is now addressed.

######  Theorem 4.

If an ODE in $n$ dimensions can be expressed in the following form

$\dot{\boldsymbol{x}}=A\boldsymbol{x}+\sum_{i=1}^{k}b_{i}(\boldsymbol{x})\boldsymbol{v}_{i}$
(99)

for the scalar functions $b_{i}(\boldsymbol{x})$ and the matrix
$V:=[\boldsymbol{v}_{1},...,\boldsymbol{v}_{k}]\in\mathbb{R}^{n\times k}$ has
rank $n-m$, where $m<n$, then there exists the linear transformations
$\underline{\boldsymbol{p}}=Q\boldsymbol{x}\in\mathbb{R}^{m}$ and
$\boldsymbol{y}=R\boldsymbol{x}\in\mathbb{R}^{n-m}$ that decouple the ODE into
an $m$ dimensional linear ODE for $\underline{\boldsymbol{p}}$ and an $n-m$
dimensional non-linear ODE for $\boldsymbol{y}$

$\displaystyle\dot{\underline{\boldsymbol{p}}}=$ $\displaystyle
L\underline{\boldsymbol{p}}$ (100) $\displaystyle\dot{\boldsymbol{y}}=$
$\displaystyle\boldsymbol{g}(\underline{\boldsymbol{p}},\boldsymbol{y})$ (101)

for some matrix $L$.

###### Proof.

Given an ODE in the form (99), then choose the linear transformation
$\underline{\boldsymbol{p}}=Q\boldsymbol{x}$ such that
$\boldsymbol{v}_{i}\in\ker(Q)$. Multiplying the ODE (99) by $Q$ gives

$Q\dot{\boldsymbol{x}}=QA\boldsymbol{x}+\sum_{i=1}^{k}b_{i}(\boldsymbol{x})Q\boldsymbol{v}_{i}=QA\boldsymbol{x}=\dot{\underline{\boldsymbol{p}}},$
(103)

which gives the form of the decoupled linear part of the ODE (100). Now we
need to construct the $n-m$ non-linear part of the ODE (101). To do this, we
simply choose some matrix $R\in\mathbb{R}^{(n-m)\times n}$ whose rows are
independent to the rows in $Q$. This then defines the one-to-one
transformation

$\boldsymbol{z}:=\left(\begin{array}[]{c}\underline{\boldsymbol{p}}\\\
\boldsymbol{y}\\\ \end{array}\right)=\left(\begin{array}[]{c}Q\\\ R\\\
\end{array}\right)\boldsymbol{x}:=G.\boldsymbol{x}$ (104)

Then the transformed ODE

$\dot{\boldsymbol{z}}=G\boldsymbol{f}(\boldsymbol{z})$

is in the desired form. ∎

A logical question now is how does one find such a linear transformation. We
will address this now by presenting a systematic algorithm to transform an ODE
given in the form (99) into its decoupled form (100) and (101).

######  Algorithm 1.

Given an ODE $\dot{\boldsymbol{x}}=\boldsymbol{f}(\boldsymbol{x})$

  1. 1.

Let $i=1$.

  2. 2.

Solve the condition $\nabla
p_{i}(\boldsymbol{x})^{T}\boldsymbol{f}(\boldsymbol{x})=\lambda
p_{i}(\boldsymbol{x})$ for affine $p_{i}(\boldsymbol{x})$ and constant
$\lambda_{i}$.

(a) If there is no solutions, or they have all been used in successive steps,
then end the algorithm.

(b) If there are one or more solutions, pick one, set $i\rightarrow i+1$ and
move to the next step.

  3. 3.

Set $p_{i-1}(\boldsymbol{x})=0$ in $\boldsymbol{f}$ and solve the condition
$\left.\nabla
p_{i}(\boldsymbol{x})^{T}\boldsymbol{f}(\boldsymbol{x})\right|_{p_{i-1}(\boldsymbol{x})=0}=\lambda_{i}p_{i}(\boldsymbol{x})$
for $p_{i}(\boldsymbol{x})$ and $\lambda_{i}$.

(a) If there is a solution, then set $i\rightarrow i+1$ and repeat step 3

(b) If there is no solution, go back to step 2 and pick a different Darboux
polynomial solution

  4. 4.

Calculate $Q=\nabla\underline{\boldsymbol{p}}(\boldsymbol{x})$

  5. 5.

Calculate $L$ from $\nabla(Q.f)=LQ$

  6. 6.

Choose a matrix $R$ s.t.

$G=\left(\begin{array}[]{c}Q\\\ R\\\ \end{array}\right)$ (105)

has full rank. Then the ODE

$\dot{\boldsymbol{z}}=G\boldsymbol{f}(\boldsymbol{z})$ (106)

is in the desired form. End algorithm.

We will now demonstrate this algorithm with an example.

######  Example 10 (Example of algorithm 1).

Consider the following vector field in five dimensions

$\frac{\mathrm{d}}{\mathrm{d}t}\left(\begin{array}[]{c}x_{1}\\\ x_{2}\\\
x_{3}\\\ x_{4}\\\ x_{5}\\\
\end{array}\right)=\left(\begin{array}[]{c}-2\,x_{{1}}x_{{3}}+5\,{x_{{2}}}^{2}+{x_{{4}}}^{2}-2\,{x_{{5}}}^{2}+x_{{2}}+x_{{5}}\\\
\left(-2\,x_{{3}}+15\right)x_{{1}}+5\,{x_{{2}}}^{2}+{x_{{4}}}^{2}-2\,{x_{{5}}}^{2}+6\,x_{{2}}+12\,x_{{3}}+8\,x_{{4}}+x_{{5}}\\\
\left(2\,x_{{3}}+2\right)x_{{1}}-5\,{x_{{2}}}^{2}-{x_{{4}}}^{2}+2\,{x_{{5}}}^{2}+x_{{3}}+2\,x_{{4}}-x_{{5}}\\\
\left(2\,x_{{3}}-7\right)x_{{1}}-5\,{x_{{2}}}^{2}-{x_{{4}}}^{2}+2\,{x_{{5}}}^{2}-3\,x_{{2}}-6\,x_{{3}}-3\,x_{{4}}-x_{{5}}\\\
x_{{1}}x_{{3}}-2\,{x_{{2}}}^{2}+{x_{{5}}}^{2}\\\ \end{array}\right)$ (107)

We will now implement algorithm 1 to decouple this ODE into a set of linear
and non-linear equations. First look for an affine Darboux polynomial that has
a constant cofactor. We find the following Darboux polynomial

$p_{1}=x_{{1}}+x_{{2}}+2\,x_{{4}}$ (108)

with cofactor 1 is the only affine Darboux polynomial that has constant
cofactor. Now eliminate a variable from the vector field by setting $p_{1}=0$,
for example, by substituting $x_{1}=-x_{2}-2x_{4}$ into $\boldsymbol{f}$.

$\frac{\mathrm{d}}{\mathrm{d}t}\left(\begin{array}[]{c}x_{2}\\\ x_{3}\\\
x_{4}\\\ x_{5}\\\
\end{array}\right)=\left(\begin{array}[]{c}\left(-2\,x_{{3}}+15\right)\left(-x_{{2}}-2\,x_{{4}}\right)+5\,{x_{{2}}}^{2}+{x_{{4}}}^{2}-2\,{x_{{5}}}^{2}+6\,x_{{2}}+12\,x_{{3}}+8\,x_{{4}}+x_{{5}}\\\
\left(2\,x_{{3}}+2\right)\left(-x_{{2}}-2\,x_{{4}}\right)-5\,{x_{{2}}}^{2}-{x_{{4}}}^{2}+2\,{x_{{5}}}^{2}+x_{{3}}+2\,x_{{4}}-x_{{5}}\\\
\left(2\,x_{{3}}-7\right)\left(-x_{{2}}-2\,x_{{4}}\right)-5\,{x_{{2}}}^{2}-{x_{{4}}}^{2}+2\,{x_{{5}}}^{2}-3\,x_{{2}}-6\,x_{{3}}-3\,x_{{4}}-x_{{5}}\\\
\left(-x_{{2}}-2\,x_{{4}}\right)x_{{3}}-2\,{x_{{2}}}^{2}+{x_{{5}}}^{2}\\\
\end{array}\right)$ (109)

We now look for Darboux polynomials with constant cofactor on the resulting
four dimensional vector field and find

$p_{2}=-x_{{2}}+x_{{3}}-2\,x_{{4}}$ (110)

also with cofactor 1. We now substitute $x_{2}=x_{3}-2x_{4}$ and compute
constant cofactor Darboux polynomials on the resulting three dimensional
system

$\frac{\mathrm{d}}{\mathrm{d}t}\left(\begin{array}[]{c}x_{3}\\\ x_{4}\\\
x_{5}\\\
\end{array}\right)=\left(\begin{array}[]{c}-\left(2\,x_{{3}}+2\right)x_{{3}}-5\,\left(x_{{3}}-2\,x_{{4}}\right)^{2}-{x_{{4}}}^{2}+2\,{x_{{5}}}^{2}+x_{{3}}+2\,x_{{4}}-x_{{5}}\\\
-\left(2\,x_{{3}}-7\right)x_{{3}}-5\,\left(x_{{3}}-2\,x_{{4}}\right)^{2}-{x_{{4}}}^{2}+2\,{x_{{5}}}^{2}-9\,x_{{3}}+3\,x_{{4}}-x_{{5}}\\\
-{x_{{3}}}^{2}-2\,\left(x_{{3}}-2\,x_{{4}}\right)^{2}+{x_{{5}}}^{2}\\\
\end{array}\right)$ (111)

this vector field has the polynomial $p_{3}=x_{4}-x_{3}$ with cofactor $1$.
Substituting $x_{3}=x_{4}$ into the above three dimensional vector field, we
find that there are no more affine Darboux polynomials with constant cofactor.
Setting
$\underline{\boldsymbol{p}}(\boldsymbol{x})=(p_{1}(x),p_{2}(x),p_{3}(x))^{T}$,
we can write $\underline{\boldsymbol{p}}(\boldsymbol{x})=Q\boldsymbol{x}$,
where

$Q=\left[\begin{array}[]{ccccc}1&1&0&2&0\\\ \vskip 6.0pt plus 2.0pt minus
2.0pt\cr 0&-1&1&-2&0\\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr
0&0&-1&1&0\end{array}\right].$ (112)

It must therefore be possible to write the original ODE as the following
decoupled set of ODEs

$\displaystyle\dot{\underline{\boldsymbol{p}}}=$ $\displaystyle
L\underline{\boldsymbol{p}},$ (113) $\displaystyle\dot{\boldsymbol{y}}=$
$\displaystyle\boldsymbol{g}(\boldsymbol{y},\underline{\boldsymbol{p}}),$
(114)

where $L$ is triangular. Multiplying the ODE
$\dot{\boldsymbol{x}}=\boldsymbol{f}(\boldsymbol{x})$ by $Q$, we get the
identity $Q\boldsymbol{f}(\boldsymbol{x})=L\underline{\boldsymbol{p}}$. In
other words, $Q\boldsymbol{f}(\boldsymbol{x})=B\boldsymbol{x}$ must be linear.
The coefficients of the matrix $L$ can therefore be easily found by solving
the linear problem $B=LQ$. Doing so, we find

$L=\left[\begin{array}[]{ccc}1&0&0\\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr
1&1&0\\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr-9&-6&1\end{array}\right].$
(115)

We can now define the linear transformation

$\left(\begin{array}[]{c}p\\\ y\\\
\end{array}\right)=\left[\begin{array}[]{c}Q\\\ R\\\
\end{array}\right]\boldsymbol{x}$ (116)

where $R$ is any $2\times n$ matrix whose rows are independent to the rows of
$Q$. We will choose $R=[0,I_{2}]$. Hence by the above linear transformation,
the ODE reads

$\displaystyle\dot{\underline{\boldsymbol{p}}}$
$\displaystyle=L\underline{\boldsymbol{p}},$ (117) $\displaystyle\dot{y_{1}}$
$\displaystyle=2\,p_{{1}}p_{{3}}-5\,{p_{{2}}}^{2}-{y_{{1}}}^{2}+2\,{y_{{2}}}^{2}-7\,p_{{1}}-3\,p_{{2}}-6\,p_{{3}}-3\,y_{{4}}-y_{{5}},$
(118) $\displaystyle\dot{y_{2}}$
$\displaystyle=p_{{1}}p_{{3}}-2\,{p_{{2}}}^{2}+{y_{{2}}}^{2}.$ (119)

∎

##  5 Conclusion

In this paper we have shown that all Runge-Kutta methods will preserve all
affine second integrals of ODEs. This implies that all Runge-Kutta methods
will preserve all rational first integrals where the numerator and denominator
are affine. Furthermore, we have demonstrated a method for detecting the
existence of such integrals in an ODE, if it exists. We have also shown that
Runge-Kutta methods cannot preserve irreducible quadratic second integrals, in
general.

The constant cofactor case is also studied, which gives rise to certain
iteration-index dependent first integrals. Due to this, one can eliminate this
iteration-index-dependence by taking quotients, which leads to a the
preservation of certain modified first integrals by Runge-Kutta methods
(corollary 3). This is the first example where an arbitrary Runge-Kutta method
preserves such an integral. For ODE systems that possess affine higher
integrals, we show that in some cases, the diagonal Padé Runge-Kutta methods
can preserve quadratic integrals exactly. Lastly, we present an algorithm that
detects the existence of such affine integrals as well as determines the
linear transformation that decouples the affine subsystem of higher integrals
from the nonlinear components.

## Acknowledgments

BK Tapley would like to thank GRW Quispel for many useful comments and
suggestions, especially with regards to the constant cofactor case.

## References

  * [1] Benedict Leimkuhler and Sebastian Reich.  Simulating hamiltonian dynamics.  Number 14. Cambridge university press, 2004.
  * [2] Feng Kang and Shang Zai-jiu.  Volume-preserving algorithms for source-free dynamical systems.  Numerische Mathematik, 71(4):451–463, Oct 1995.
  * [3] Ernst Hairer, Christian Lubich, and Gerhard Wanner.  Geometric numerical integration: structure-preserving algorithms for ordinary differential equations, volume 31.  Springer Science & Business Media, 2006.
  * [4] Morten Dahlby, Brynjulf Owren, and Takaharu Yaguchi.  Preserving multiple first integrals by discrete gradients.  Journal of Physics A: Mathematical and Theoretical, 44(30):305205, 2011.
  * [5] GJ967831 Cooper.  Stability of runge-kutta methods for trajectory problems.  IMA journal of numerical analysis, 7(1):1–13, 1987.
  * [6] Elena Celledoni, Robert I McLachlan, David I McLaren, Brynjulf Owren, G Reinout W Quispel, and William M Wright.  Energy-preserving runge-kutta methods.  ESAIM: Mathematical Modelling and Numerical Analysis, 43(4):645–649, 2009.
  * [7] GRW Quispel, HW Capel, et al.  Solving ode’s numerically while preserving all first integrals.  Preprint, La Trobe University, Melbourne, 1997.
  * [8] GRW Quispel and David Ian McLaren.  A new class of energy-preserving numerical integration methods.  Journal of Physics A: Mathematical and Theoretical, 41(4):045206, 2008.
  * [9] Ernst Hairer.  Energy-preserving variant of collocation methods.  Journal of Numerical Analysis, Industrial and Applied Mathematics, 5:73–84, 2010.
  * [10] Xinyuan Wu, Bin Wang, and Wei Shi.  Efficient energy-preserving integrators for oscillatory hamiltonian systems.  Journal of Computational Physics, 235:587–605, 2013.
  * [11] Alain Goriely.  Integrability and nonintegrability of dynamical systems, volume 19.  World Scientific, 2001.
  * [12] David W Albrecht, Elizabeth L Mansfield, and Alice E Milne.  Algorithms for special integrals of ordinary differential equations.  Journal of Physics A: Mathematical and General, 29(5):973, 1996\.
  * [13] CG Hewitt.  Algebraic invariant curves in cosmological dynamical systems and exact solutions.  General relativity and gravitation, 23(12):1363–1383, 1991.
  * [14] Jean-Pierre Jouanolou.  Equations de Pfaff algébriques, volume 708.  Springer, 2006.
  * [15] Jean Moulin Ollagnier, Andrzej Nowicki, and Jean-Marie Strelcyn.  On the non-existence of constants of derivations: the proof of a theorem of Jouanolou and its development.  Laboratoire analyse, géométrie et applications, Unité associée au CNRS 742, 1993.
  * [16] Elena Celledoni, Charalambos Evripidou, David McLaren, Brynjulf Owren, Reinout Quispel, and Benjamin Tapley.  Discrete darboux polynomials and the search for preserved measures and integrals of rational maps.  arXiv preprint arXiv:1902.04685, 2019.
  * [17] E Celledoni, C Evripidou, DI McLaren, B Owren, GRW Quispel, BK Tapley, and PH van der Kamp.  Using discrete darboux polynomials to detect and determine preserved measures and integrals of rational maps.  Journal of Physics A: Mathematical and Theoretical, 52(31):31LT01, 2019.
  * [18] Elena Celledoni, Robert I McLachlan, Brynjulf Owren, and G Reinout W Quispel.  Geometric properties of kahan’s method.  Journal of Physics A: Mathematical and Theoretical, 46(2):025201, 2012.
  * [19] CB Collins.  Algebraic conditions for a centre or a focus in some simple systems of arbitrary degree.  Journal of mathematical analysis and applications, 195(3):719–735, 1995.
  * [20] Robert I McLachlan, GRW Quispel, et al.  Six lectures on the geometric integration of ODEs.  Citeseer, 1998.
  * [21] Gerhard Wanner and Ernst Hairer.  Solving ordinary differential equations II.  Springer Berlin Heidelberg, 1996.
  * [22] Michael F Reusch, Lee Ratzan, Neil Pomphrey, and Wonchull Park.  Diagonal padé approximations for initial value problems.  SIAM journal on scientific and statistical computing, 9(5):829–838, 1988.

##  Appendix A Derivation of stability matrix for corollary 2

Let $\varphi_{h}$ denote the Runge-Kutta map defined by equations (4), (5) and
(6). Now let $G=(g_{1},g_{2},...,g_{s})^{T}$ denote the $s\times n$ matrix
whose $i$’th row is $g_{i}^{T}$. Then if $f(x)=\lambda x$ as in the case of
our test equation (5) can be written as

$\displaystyle
G=\mathbb{1}x^{T}+h\lambda\mathcal{A}G=(I-h\lambda\mathcal{A})^{-1}\mathbb{1}x^{T}$
(120)

We therefore have

$g_{i}=G^{T}\hat{e_{i}}=\left(\hat{e_{i}}^{T}(I-h\lambda\mathcal{A})^{-1}\mathbb{1}\right)x$
(121)

where $\hat{e_{i}}\in\mathbb{R}^{s}$ is the $i$’th canonical unit basis vector
with a 1 in the $i$’th component and 0 elsewhere. Inserting this into equation
(6) gives

$\displaystyle\varphi_{h}(x)=$ $\displaystyle
x+h\lambda\sum_{j=1}^{s}b_{j}\hat{e_{i}}^{T}(I-h\lambda\mathcal{A})^{-1}\mathbb{1}x$
(122) $\displaystyle=$ $\displaystyle(1+h\lambda
b^{T}(I-h\lambda\mathcal{A})^{-1}\mathbb{1})x$ (123) $\displaystyle:=$
$\displaystyle R(\lambda h)x$ (124)


elaboration

# Reading and Acting while Blindfolded:
The Need for Semantics in Text Game Agents

Shunyu Yao† Karthik Narasimhan† Matthew Hausknecht‡
†Princeton University ‡Microsoft Research
{shunyuy<EMAIL_ADDRESS>
<EMAIL_ADDRESS>
Work partly done during internship at Microsoft Research. Project page:
https://blindfolded.cs.princeton.edu

###### Abstract

Text-based games simulate worlds and interact with players using natural
language. Recent work has used them as a testbed for autonomous language-
understanding agents, with the motivation being that understanding the
meanings of words or semantics is a key component of how humans understand,
reason, and act in these worlds. However, it remains unclear to what extent
artificial agents utilize semantic understanding of the text. To this end, we
perform experiments to systematically reduce the amount of semantic
information available to a learning agent. Surprisingly, we find that an agent
is capable of achieving high scores even in the complete absence of language
semantics, indicating that the currently popular experimental setup and models
may be poorly designed to understand and leverage game texts. To remedy this
deficiency, we propose an inverse dynamics decoder to regularize the
representation space and encourage exploration, which shows improved
performance on several games including Zork I. We discuss the implications of
our findings for designing future agents with stronger semantic understanding.

##  1 Introduction

(a) Zork I _Observation 21:_ You are in the living room. There is a doorway to
the east, a wooden door with strange gothic lettering to the west, which
appears to be nailed shut, a trophy case, and a large oriental rug in the
center of the room. You are carrying: A brass lantern … _Action 21_ : move rug
_Observation 22_ : With a great effort, the rug is moved to one side of the
room, revealing the dusty cover of a closed trap door… Living room… You are
carrying: … _Action 22_ : open trap

(b) min-ob _Observation 21:_ Living Room _Action 21_ : move rug _Observation
22_ : Living Room _Action 22_ : open trap

(c) hash _Observation 21:_ 0x6fc _Action 21_ : 0x3a04 _Observation 22_ :
0x103b _Action 22_ : 0x16bb

Figure 1: (a): Sample original gameplay from Zork I. (b) (c): Our proposed
semantic ablations. (b) min-ob reduces observations to only the current
location name, and (c) hash replaces observation and action texts by their
string hash values.

Text adventure games such as Zork I (Figure 1 (a)) have been a testbed for
developing autonomous agents that operate using natural language. Since
interactions in these games (input observations, action commands) are through
text, the ability to understand and use language is deemed necessary and
critical to progress through such games. Previous work has deployed a spectrum
of methods for language processing in this domain, including word vectors
Fulda et al. (2017), recurrent neural networks Narasimhan et al. (2015);
Hausknecht et al. (2020), pre-trained language models Yao et al. (2020), open-
domain question answering systems Ammanabrolu et al. (2020), knowledge graphs
Ammanabrolu and Hausknecht (2020); Ammanabrolu et al. (2020); Adhikari et al.
(2020), and reading comprehension systems Guo et al. (2020).

Meanwhile, most of these models operate under the reinforcement learning (RL)
framework, where the agent explores the same environment in repeated episodes,
learning a value function or policy to maximize game score. From this
perspective, text games are just special instances of a partially observable
Markov decision process (POMDP) $(S,T,A,O,R,\gamma)$, where players issue text
actions $a\in A$, receive text observations $o\in O$ and scalar rewards
$r=R(s,a)$, and the underlying game state $s\in S$ is updated by transition
$s^{\prime}=T(s,a)$.

However, what distinguishes these games from other POMDPs is the fact that the
actions and observations are in language space $L$. Therefore, a certain level
of decipherable _semantics_ is attached to text observations $o\in O\subset L$
and actions $a\in A\subset L$. Ideally, these texts not only serve as
observation or action _identifiers_ , but also provide clues about the latent
transition function $T$ and reward function $R$. For example, issuing an
action “jump” based on an observation “on the cliff” would likely yield a
subsequent observation such as “you are killed” along with a negative reward.
Human players often rely on their understanding of language semantics to
inform their choices, even on games they have never played before, while
replacing texts with non-semantic identifiers such as their corresponding
hashes (Figure 1 (c)) would likely render games unplayable for people.
However, would this type of transformation affect current RL agents for such
games? In this paper, we ask the following question: _To what extent do
current reinforcement learning agents leverage semantics in text-based games?_

To shed light on this question, we investigate the Deep Reinforcement
Relevance Network (DRRN) He et al. (2016), a top-performing RL model that uses
gated recurrent units (GRU) Cho et al. (2014) to encode texts. We conduct
three experiments on a set of interactive fiction games from the Jericho
benchmark Hausknecht et al. (2020) to probe the effect of different semantic
representations on the functioning of DRRN. These include (1) using just a
location phrase as the input observation (Figure 1 (b)), (2) hashing text
observations and actions (Figure 1 (c)), and (3) regularizing vector
representations using an auxiliary inverse dynamics loss. While reducing
observations to location phrases leads to decreased scores and enforcing
inverse dynamics decoding leads to increased scores on some games, hashing
texts to break semantics surprisingly matches or even outperforms the baseline
DRRN on almost all games considered. This implies current RL agents for text-
based games might not be sufficiently leveraging the semantic structure of
game texts to learn good policies, and points to the need for developing
better experiment setups and agents that have a finer grasp of natural
language.

##  2 Models

### DRRN Baseline

Our baseline RL agent DRRN He et al. (2016) learns a Q-network $Q_{\phi}(o,a)$
parametrized by $\phi$. The model encodes the observation $o$ and each action
candidate $a$ using two separate GRU encoders $f_{o}$ and $f_{a}$, and then
aggregates the representations to derive the Q-value through a MLP decoder
$g$:

$Q_{\phi}(o,a)=g(\mathrm{concat}(f_{o}(o),f_{a}(a)))\vspace{-1.5mm}$ (1)

For learning $\phi$, tuples $(o,a,r,o^{\prime})$ of observation, action,
reward and the next observation are sampled from an experience replay buffer
and the following temporal difference (TD) loss is minimized:

$\mathcal{L}_{\text{{TD}}}(\phi)=(r+\gamma\max_{a^{\prime}\in
A}Q_{\phi}(o^{\prime},a^{\prime})-Q_{\phi}(o,a))^{2}\vspace{-1.5mm}$ (2)

During gameplay, a softmax exploration policy is used to sample an action:

$\pi_{\phi}(a|o)=\frac{\exp(Q_{\phi}(o,a))}{\sum_{a^{\prime}\in
A}\exp(Q_{\phi}(o,a^{\prime}))}\vspace{-1.5mm}$ (3)

Note that when the action space $A$ is large, (2) and (3) become intractable.
A valid action handicap Hausknecht et al. (2020) or a language model Yao et
al. (2020) can be used to generate a reduced action space for efficient
exploration. For all the modifications below, we use the DRRN with the valid
action handicap as our base model.

### Reducing Semantics via Minimizing Observation (min-ob)

Unlike other RL domains such as video games or robotics control, at each step
of text games the (valid) action space is constantly changing, and it reveals
useful information about the current state. For example, knowing “unlock box”
is valid leaks the existence of a locked box. Also, sometimes action semantics
indicate its value even unconditional on the state, e.g. “pick gold” usually
seems good. Given these, we minimize the observation to only a location phrase
$o\mapsto\mathrm{loc}(o)$ (Figure 1 (b)) to isolate the action semantics:
$Q^{\mathrm{loc}}_{\phi}(o,a)=g(f_{o}(\mathrm{loc}(o)),f_{a}(a)))$.

### Breaking Semantics via Hashing (hash)

GRU encoders $f_{o}$ and $f_{a}$ in the Q-network (1) generally ensure that
similar texts (e.g. a single word change) are given similar representations,
and therefore similar values. To study whether such a semantics continuity is
useful, we break it by hashing observation and action texts. Specifically,
given a hash function from strings to integers $h:L\to\mathbb{Z}$, and a
pseudo-random generator $G:\mathbb{Z}\to\mathbb{R}^{d}$ that turns an integer
seed to a random Gaussian vector, a hashing encoder $\hat{f}=G\circ
h:L\to\mathbb{R}^{d}$ can be composed. While $f_{o}$ and $f_{a}$ in (1) are
trainable, $\hat{f}$ is fixed throughout RL, and ensures two texts that only
differ by a word would have completely different representations. In this
sense, hashing breaks semantics and only serves to identify different
observations and actions in an abstract MDP problem (Figure 1 (c)):
$Q^{\mathrm{hash}}_{\phi}(o,a)=g(\hat{f}(o),\hat{f}(a))$.

Game | DRRN | min-ob | hash | inv-dy | Max
---|---|---|---|---|---
balances | 10 / 10 | 10 / 10 | 10 / 10 | 10 / 10 | 51
deephome | 57 / 66 | 8.5 / 27 | 58 / 67 | 57.6 / 67 | 300
detective | 290 / 337 | 86.3 / 350 | 290 / 317 | 290 / 323 | 360
dragon | -5.0 / 6 | -5.4 / 3 | -5.0 / 7 |  -2.7 / 8 | 25
enchanter | 20 / 20 | 20 / 40 | 20 / 30 | 20 / 30 | 400
inhumane | 21.1 / 45 | 12.4 / 40 | 21.9 / 45 | 19.6 / 45 | 90
library | 15.7 / 21 | 12.8 / 21 |  17 / 21 | 16.2 / 21 | 30
ludicorp | 12.7 / 23 | 11.6 / 21 |  14.8 / 23 | 13.5 / 23 | 150
omniquest | 4.9 / 5 | 4.9 / 5 | 4.9 / 5 | 5.3 / 10 | 50
pentari | 26.5 / 45 | 21.7 / 45 |  51.9 / 60 | 37.2 / 50 | 70
zork1 | 39.4 / 53 | 29 / 46 | 35.5 / 50 |  43.1 / 87 | 350
zork3 | 0.4 / 4.5 | 0.0 / 4 | 0.4 / 4 | 0.4 / 4 | 7
Avg. Norm | .21 / .38 | .12 / .35 |  .25 / .39 | .23 / .40 |

Table 1: Final/maximum score of different models.

### Regularizing Semantics via Inverse Dynamics Decoding (inv-dy)

The GRU representations in DRRN $f_{o}(o),f_{a}(a)$ are only optimized for the
TD loss (2). As a result, text semantics can degenerate during encoding, and
the text representation might arbitrarily overfit to the Q-values. To
regularize and encourage more game-related semantics to be encoded, we take
inspiration from Pathak et al. (2017) and propose an inverse dynamics
auxiliary task during RL. Given representations of current and next
observations $f_{o}(o),f_{o}(o^{\prime})$, we use a MLP $g_{inv}$ to predict
the action representation, and a GRU decoder $d$ to decode the action back to
text***Directly defining an L1/L2 loss between $f_{a}(a)$ and
$g_{inv}(\mathrm{concat}(f_{o}(o),f_{o}(o^{\prime})))$ in the representation
space will collapse text representations together.. The inverse dynamics loss
is defined as

$\mathcal{L}_{\text{{inv}}}(\phi,\theta)=-\log
p_{d}(a|g_{inv}(\mathrm{concat}(f_{o}(o),f_{o}(o^{\prime})))\vspace{-1.5mm}$

where $\theta$ denote weights of $g_{inv}$ and $d$, and $p_{d}(a|x)$ is the
probability of decoding token sequence $a$ using GRU decoder $d$ with initial
hidden state as $x$. To also regularize the action encoding, action
reconstruction from $f_{a}$ is also used as a loss term:

$\mathcal{L}_{\text{{dec}}}(\phi,\theta)=-\log
p_{d}(a|f_{a}(a))\vspace{-1.5mm}$

And during experience replay, these two losses are optimized along with the TD
loss:

$\mathcal{L}(\phi,\theta)=\mathcal{L}_{\text{{TD}}}(\phi)+\lambda_{1}\mathcal{L}_{\text{{inv}}}(\phi,\theta)+\lambda_{2}\mathcal{L}_{\text{{dec}}}(\phi,\theta)\vspace{-1.5mm}$

Figure 2: Transfer results from Zork I.

An intrinsic reward $r^{+}=\mathcal{L}_{\text{{inv}}}(\phi,\theta)$ is also
used to explore toward where the inverse dynamics is not learned well yet. All
in all, the aim of inv-dy is threefold: (1) regularize both action and
observation representations to avoid degeneration by decoding back to the
textual domain, (2) encourage $f_{o}$ to encode action-relevant parts of
observations, and (3) provide intrinsic motivation for exploration.

##  3 Results

### Setup

We train on 12 games†††We omit games where DRRN cannot score. from the Jericho
benchmark Hausknecht et al. (2020). These human-written interactive fictions
are rich, complex, and diverse in semantics.‡‡‡Please refer to Hausknecht et
al. (2020) for more details about these games. For each game, we train DRRN
asynchronously on 8 parallel instances of the game environment for $10^{5}$
steps, using a prioritized replay buffer. Following prior practice Hausknecht
et al. (2020), we augment observations with location and inventory
descriptions by issuing the ‘look’ and ‘inventory’ commands. We train three
independent runs for each game and report their average score. For hash, we
use the Python built-in hash function to process text as a tuple of token IDs,
and implement the random vector generator $G$ by seeding PyTorch with the hash
value. For inv-dy, we use $\lambda_{1}=\lambda_{2}=1$.

Figure 3: t-SNE visualization of seen and unseen state observations of Zork I.
DRRN (base) represents unseen states separated from seen states while inv-dy
mixes them by semantic similarity.

### Scores

Table 1 reports the final score (the average score of the final 100 episodes
during training), and the maximum score seen in each game for different
models. Average normalized score (raw score divided by game total score) over
all games is also reported. Compared to the base DRRN, min-ob turns out to
explore similar maximum scores on most games (except Deephome and Dragon), but
fails to memorize the good experience and reach high episodic scores, which
suggests the importance of identifying different observations using language
details. Most surprisingly, hash has a lower final score than DRRN on only one
game (Zork I), while on Pentari it almost doubles the DRRN final score. It is
also the model with the best average normalized final score across games,
which indicates that the DRRN model can perform as well without leveraging any
language semantics, but instead simply by identifying different observations
and actions with random vectors and memorizing the Q-values. Lastly, we
observe on some games (Dragon, Omniquest, Zork I) inv-dy can explore high
scores that other models cannot. Most notably, on Zork I the maximum score
seen is 87 (average of $54,94,113$), while any run of other models does not
explore a score more than 55. This might indicate potential benefit of
developing RL agents with more semantic representations.

### Transfer

We also investigate if representations of different models can transfer to a
new language environment, which is a potential benefit of learning natural
language semantics. So we consider the two most similar games in Jericho, Zork
I and Zork III, fix the language encoders of different Zork I models, and re-
train the Q-network on Zork III for 10,000 steps. As shown in Figure 2, inv-dy
representations can achieve a score around 1, which surpasses the best result
of models trained from scratch on Zork III for 100,000 steps (around 0.4),
showing great promise in better gameplay by leveraging language understanding
from other games. hash transfer is equivalent to training from scratch as the
representations are not learnt, and a score around 0.3 is achieved. Finally,
DRRN representations transfer worse than hash, possibly due to overfitting to
the TD loss (2).

### Visualizations

Finally, we use t-SNE Maaten and Hinton (2008) to visualize representations of
some Zork I walkthrough states in Figure 3. The first 30 walkthrough states
(red, score 0-45) are well experienced by the models during exploration,
whereas the last 170 states (blue, score 157-350) are unseen§§§The rest 150
states in the middle (score 45-157) are omitted as they might be seen by some
model but not others.. We also encircle the subset of states at location
‘living room’ for their shared semantics.

First, we note that the hash representations for living room states are
scattered randomly, unlike the other two models with GRU language encoders.
Further, the base DRRN overfits to the TD loss (2), representing unseen states
(blue) in a different subspace to seen states (red) without regarding their
semantic similarity. ind-dy is able to extrapolate to unseen states and
represent them similarly to seen states for their shared semantics, which may
explain its better performance on this game.

### Game stochasticity

All the above experiments were performed using a fixed game random seed for
each game, following prior work Hausknecht et al. (2020). To investigate if
randomness in games affects our conclusions, we run one trial of each game
with episode-varying random seeds¶¶¶Randomness includes transition uncertainty
(e.g. thief showing up randomly in Zork I) and occasional paraphrasing of text
observations.. We find the average normalized score for the base DRRN, hash,
inv-dy to be all around 17%, with performance drop mainly on three stochastic
games (Dragon, Zork I, Zork III). Notably, the core finding that the base DRRN
and hash perform similarly still holds. Intuitively, even though the Q-values
would be lower overall with unexpected transitions, RL would still memorize
observations and actions that lead to high Q-values.

##  4 Discussion

At a high level, RL agents for text-based games succeed by (1) exploring
trajectories that lead to high scores, and (2) learning representations to
stably reach high scores. Our experiments show that a semantics-regularized
inv-dy model manages to explore higher scores on some games (Dragon,
Omniquest, Zork I), while the hash model manages to memorize scores better on
other games (Library, Ludicorp, Pentari) using just a fixed, random, non-
semantic representation. This leads us to hypothesize two things. First,
fixed, stable representations might make Q-learning easier. Second, it might
be desirable to represent similar texts very differently for better gameplay,
e.g. the Q-value can be much higher when a key object is mentioned, even if it
only adds a few words to a long observation text. This motivates future
thought into the structural vs. functional use of language semantics in these
games.

Our findings also urge a re-thinking of the popular ‘RL + valid action
handicap’ setup for these games. On one hand, RL sets training and evaluation
in the same environment, with limited text corpora, and sparse, mostly
deterministic rewards as the only optimization objective. Such a combination
easily results in overfitting to the reward system of a specific game (Figure
2), or even just a specific stage of the game (Figure 3). On the other hand,
the valid action handicap reduces the action set to a small size tractable for
memorization, and reduces the language understanding challenge for the RL
agent. Thus for future research on text-based games, we advocate for more
attention towards alternative setups without RL or handicaps Hausknecht et al.
(2019); Yao et al. (2020); Yu et al. (2020). Particularly, in a ‘RL + no valid
action handicap’ setting, _generating_ action candidates rather than simply
_choosing_ from a set entails more opportunities and challenges with respect
to learning grounded language semantics Yao et al. (2020). Additionally,
training agents on a distribution of games and evaluating them on a separate
set of unseen games would require more general semantic understanding.
Semantic evaluation of these proposed paradigms is outside the scope of this
paper, but we hope it will spark a productive discussion on the next steps
toward building agents with stronger semantic understanding.

## Ethical Considerations

Autonomous decision-making agents are potentially impactful in our society,
and it is of great ethical consideration to make sure their understanding of
the world and their objectives align with humans. Humans use natural language
to convey and understand concepts as well as inform decisions, and in this
work we investigate whether autonomous agents leverage language semantics
similarly to humans in the environment of text-based games. Our findings
suggest that the current generation of agents optimized for reinforcement
learning objectives might not exhibit human-like language understanding, a
phenomenon we should pay attention to and further study.

## Acknowledgements

We appreciate helpful suggestions from anonymous reviewers as well as members
of the Princeton NLP Group and MSR RL Group.

## References

  * Adhikari et al. (2020) Ashutosh Adhikari, Xingdi (Eric) Yuan, Marc-Alexandre Côté, Mikulas Zelinka, Marc-Antoine Rondeau, Romain Laroche, Pascal Poupart, Jian Tang, Adam Trischler, and William L. Hamilton . 2020.  Learning dynamic knowledge graphs to generalize on text-based games.  In _NeurIPS 2020_.
  * Ammanabrolu and Hausknecht (2020) Prithviraj Ammanabrolu and Matthew J. Hausknecht. 2020.  Graph constrained reinforcement learning for natural language action spaces.  In _8th International Conference on Learning Representations, ICLR 2020, Addis Ababa, Ethiopia, April 26-30, 2020_. OpenReview.net.
  * Ammanabrolu et al. (2020) Prithviraj Ammanabrolu, Ethan Tien, Zhaochen Luo, and Mark O Riedl. 2020.  How to avoid being eaten by a grue: Exploration strategies for text-adventure agents.  _arXiv preprint arXiv:2002.08795_.
  * Cho et al. (2014) Kyunghyun Cho, Bart van Merriënboer, Dzmitry Bahdanau, and Yoshua Bengio. 2014\.  On the properties of neural machine translation: Encoder–decoder approaches.  In _Proceedings of SSST-8, Eighth Workshop on Syntax, Semantics and Structure in Statistical Translation_ , pages 103–111, Doha, Qatar. Association for Computational Linguistics.
  * Fulda et al. (2017) Nancy Fulda, Daniel Ricks, Ben Murdoch, and David Wingate. 2017.  What can you do with a rock? affordance extraction via word embeddings.  In _Proceedings of the Twenty-Sixth International Joint Conference on Artificial Intelligence, IJCAI 2017, Melbourne, Australia, August 19-25, 2017_ , pages 1039–1045. ijcai.org.
  * Guo et al. (2020) Xiaoxiao Guo, Mo Yu, Yupeng Gao, Chuang Gan, Murray Campbell, and Shiyu Chang. 2020\.  Interactive fiction game playing as multi-paragraph reading comprehension with reinforcement learning.  In _Proceedings of the 2020 Conference on Empirical Methods in Natural Language Processing (EMNLP)_ , pages 7755–7765, Online. Association for Computational Linguistics.
  * Hausknecht et al. (2019) Matthew Hausknecht, Ricky Loynd, Greg Yang, Adith Swaminathan, and Jason D Williams. 2019.  Nail: A general interactive fiction agent.  _arXiv preprint arXiv:1902.04259_.
  * Hausknecht et al. (2020) Matthew J. Hausknecht, Prithviraj Ammanabrolu, Marc-Alexandre Côté, and Xingdi Yuan. 2020.  Interactive fiction games: A colossal adventure.  In _The Thirty-Fourth AAAI Conference on Artificial Intelligence, AAAI 2020, The Thirty-Second Innovative Applications of Artificial Intelligence Conference, IAAI 2020, The Tenth AAAI Symposium on Educational Advances in Artificial Intelligence, EAAI 2020, New York, NY, USA, February 7-12, 2020_ , pages 7903–7910. AAAI Press.
  * He et al. (2016) Ji He, Jianshu Chen, Xiaodong He, Jianfeng Gao, Lihong Li, Li Deng, and Mari Ostendorf. 2016.  Deep reinforcement learning with a natural language action space.  In _Proceedings of the 54th Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers)_ , pages 1621–1630, Berlin, Germany. Association for Computational Linguistics.
  * Maaten and Hinton (2008) Laurens van der Maaten and Geoffrey Hinton. 2008.  Visualizing data using t-sne.  _Journal of machine learning research_ , 9(Nov):2579–2605.
  * Narasimhan et al. (2015) Karthik Narasimhan, Tejas Kulkarni, and Regina Barzilay. 2015.  Language understanding for text-based games using deep reinforcement learning.  In _Proceedings of the 2015 Conference on Empirical Methods in Natural Language Processing_ , pages 1–11, Lisbon, Portugal. Association for Computational Linguistics.
  * Pathak et al. (2017) Deepak Pathak, Pulkit Agrawal, Alexei A. Efros, and Trevor Darrell. 2017.  Curiosity-driven exploration by self-supervised prediction.  In _Proceedings of the 34th International Conference on Machine Learning, ICML 2017, Sydney, NSW, Australia, 6-11 August 2017_ , volume 70 of _Proceedings of Machine Learning Research_ , pages 2778–2787. PMLR.
  * Yao et al. (2020) Shunyu Yao, Rohan Rao, Matthew Hausknecht, and Karthik Narasimhan. 2020.  Keep CALM and explore: Language models for action generation in text-based games.  In _Proceedings of the 2020 Conference on Empirical Methods in Natural Language Processing (EMNLP)_ , pages 8736–8754, Online. Association for Computational Linguistics.
  * Yu et al. (2020) Mo Yu, Xiaoxiao Guo, Yufei Feng, Xiaodan Zhu, Michael Greenspan, and Murray Campbell. 2020.  Deriving commonsense inference tasks from interactive fictions.  _arXiv preprint arXiv:2010.09788_.