\section{Order conditions} \label{order_conditions}

\hskip\parindent{}The element-wise product of column vectors $\VEC{x}$ 
and $\VEC{y}$ will be denoted as $\VEC{x} \mkern0.5mu . \mkern1mu 
\VEC{y}$, \IE, $(\VEC{x} \mkern0.5mu . \mkern1mu \VEC{y})_i = x_i 
\mkern1mu y_i$. The element-wise product of $n$ copies of vector 
$\VEC{x}$ will be written as $\VEC{x}^{\mkern1mu n}$. Let $\VEC{1}$ be 
the $s$-dimensional column vector with all components being equal to 
$1$; $\MAT{A} = \bigl[ {a_{\mkern0.85mu i \mkern-0.85mu j}} \bigr]$ be the $s \times s$ matrix with 
${a_{\mkern0.85mu i \mkern-0.85mu j}}$ as its matrix element in the $i^{\mkern1mu \textrm{th}}$ row and 
$j^{\mkern1mu \textrm{th}}$ column; $\VEC{b} = \bigl[ {b_{\mkern-1.25mu j}} \bigr]$ be the 
weights row vector; and $\VEC{c} = \bigl[ {c_{\mkern0.5mu i}} \bigr]$ be the nodes 
vector. Let $\VEC{q}_n = \MAT{A} \mkern1mu \VEC{c}^{\mkern0.5mu n} - 
\frac{1}{n + 1} \mkern1mu \VEC{c}^{\mkern0.5mu n + 1}$, \EG, $\VEC{q}_0 
= \VEC{0}$ as $\MAT{A} \VEC{1} = \VEC{c}$.

Given rooted trees $\textrm{t}_1$, $\textrm{t}_2$, {...}} \def\DS{\displaystyle, 
$\textrm{t}_n$, a new tree $[ \mkern1mu \textrm{t}_1 \mkern2mu 
\textrm{t}_2 \mkern2mu {...}} \def\DS{\displaystyle \mkern2mu \textrm{t}_n \mkern1mu ]$ is 
obtained by connecting with $n$ edges their roots to a new vertex, the 
latter becomes a new root \citep[s.~301]{But16}, \citep[p.~152]{HNW93}, 
\citep[p.~44]{But21}, \citep[p.~53]{HLW06}. Consider a vector function 
$\VEC{\Phi} : \mathrm{T} \to \pmb{\mathbf{R}}^s$ on the set of rooted 
trees that is recursively defined as $\VEC{\Phi}( 
\begin{picture}(4,6)(0,0) \put(3.35,2.5){\makebox(0,0){\circle*{3}}} 
\end{picture} ) = \VEC{1}$ and $\VEC{\Phi} \bigl( [ \mkern1mu 
\textrm{t}_1 \mkern2mu \textrm{t}_2 \mkern2mu {...}} \def\DS{\displaystyle \mkern2mu 
\textrm{t}_n \mkern1mu ] \bigr) = \prod_{\mkern1mu m = 1}^{\mkern1mu n} 
{\MAT{A} \mkern0.5mu \VEC{\Phi}}(\textrm{t}_m)$, where the product of vectors is taken 
element-wise. This function coincides with \emph{derivative weights} 
\citep[def.~312A]{But16}, \citep[pp.~148 and 151]{HNW93}, 
\citep[def.~2.2]{AuCh98a}, it is closely related to \emph{internal} or 
\emph{stage weights} ${\MAT{A} \mkern0.5mu \VEC{\Phi}}(\textrm{t})$ \citep[p.~125]{But21} and 
\emph{elementary weights} ${\VEC{b} \mkern2mu \VEC{\Phi}}(\textrm{t})$ \citep[p.~55]{HLW06}.

Representation of analytical expressions by pictures or diagrams is 
often convenient and illustrative, and is used in many different 
contexts: countless diagram techniques in theoretical physics; group 
theory (see, \EG, \citep[]{Cvi08}); graphical representation of Fredholm 
series \citep[s.~11.3]{MaWa70}; and connection of derivatives with 
rooted trees \citep[]{Cay57}, which was heavily developed in the context 
of Runge--Kutta methods \citep[]{Mer57}, \citep[]{But72}, 
\citep[ch.~3]{But16}, \citep[s.~II.2]{HNW93}, \citep[ch.~2]{But21}, 
\citep[ch.~III]{HLW06}, with occasional usage of picturesque trees in 
diagrams \citep[pp.~188 and 189]{But16}, \citep[pp.~175 and 328]{HNW93}.

Similarly to \citep[]{Pen71} graphical notation, an object that is a 
tensor of type $(m ,\mkern2mu} \def\DWS{\mkern3.5mu n)$ will be depicted by a symbol of the object 
connected by lines to $m$ and $n$ black dots below and above it, 
respectively:$\mkern1mu$\footnote{~Such an 
up$\mkern2mu$/$\mkern0.75mu$down orientation is the reverse of the one 
in \citep[]{Pen71}, and is used here so that a tree has its root at the 
bottom. Not drawing the dots allows an elegant representation of 
Kronecker delta $\delta^{\mkern1mu a}_{\mkern1mu b}$ and metric tensors 
$g^{\mkern1mu ab}$ and $g_{\mkern1mu ab}$ \citep[pp.~226 and 
231]{Pen71}, while here the dots correspond to the vertices of rooted 
trees.}${}^{,\mkern2mu} \def\DWS{\mkern3.5mu}$\footnote{~The drawings of a leaf, a root, of trees 
\emph{Betula pendula} and \emph{Larix sibirica} are by Olga Stepanova, 
used with permission.}
 \begin{align*}
   \begin{tabular}{r|l}
   \begin{picture}(17,23)(0,0)
     \put(0.9,2.5){\includegraphics[height=19pt]{leaf.pdf}}
     \put(8,1.5){\makebox(0,0){\circle*{3}}}
   \end{picture} & \raisebox{3.67pt}{vector $\VEC{1}$} \\
   \hline
   \begin{picture}(60,55)(0,0)
     \put(1.5,44){\makebox(0,0){$\textrm{t}$}}
     \put(2.55,1.95){\includegraphics[height=49pt]{birch.pdf}}
     \put(15,1.5){\makebox(0,0){\circle*{3}}}
     \put(34,1.5){\includegraphics[height=49pt]{larch.pdf}}
     \put(48,1.5){\makebox(0,0){\circle*{3}}}
   \end{picture}
     & \raisebox{12pt}{$\VEC{\Phi}(\textrm{t})$ or an arbitrary vector} \\
   \hline
   \begin{picture}(21,24)(0,0)
     \put(1.5,10){\makebox(0,0){\rotatebox{90}{\rule[3pt]{17pt}{1.44pt}}}}
     \put(1.5,1.5){\makebox(0,0){\circle*{3}}}
     \put(1.5,18.5){\makebox(0,0){\circle*{3}}}
     \put(16.5,10){\makebox(0,0){\rotatebox{90}{\rule[3pt]{17pt}{1.44pt}}}}
     \put(16.5,10){\makebox(0,0){\rotatebox{90}{\color{white}\rule[3pt]{17pt}{0.72pt}}}}
     \put(16.5,1.5){\makebox(0,0){\circle*{3}}}
     \put(16.5,18.5){\makebox(0,0){\circle*{3}}}
   \end{picture}
     & \raisebox{3.67pt}{matrix $\MAT{A}$ and identity matrix $\MAT{I}$} \\
   \hline
   \begin{picture}(16,29)(0,0)
    \put(0.4,0.5){\includegraphics[height=22pt]{root.pdf}}
     \put(5.85,23.5){\makebox(0,0){\circle*{3}}}
   \end{picture}
     & \raisebox{6.67pt}{weights row vector $\VEC{b}$} \\
   \hline
   \begin{picture}(22,17)(0,0)
     \put(11,1.5){\makebox(0,0){\circle*{3}}}
     \put(6,7){\makebox(0,0){\rotatebox{135}{\rule[3pt]{14pt}{0.72pt}}}}
     \put(15.5,5.5){\makebox(0,0){\rotatebox{45}{\rule[3pt]{14.5pt}{0.72pt}}}}
     \put(8.75,7.4){\makebox(0,0){\rotatebox{112.5}{\rule[3pt]{12pt}{0.72pt}}}}
     \put(11.75,10.5){\makebox(0,0){\footnotesize$...$}}
   \end{picture}
     & \raisebox{1.33pt}{element-wise product of vectors}
   \end{tabular} \end{align*} Lines coming both from below and from 
above to a dot produce a tensor contraction, \IE, the summation over the 
values of the corresponding tensor index. Outer products are obtained by 
simply drawing objects next to each other \citep[p.~224]{Pen71}.

In bra-ket notation \citep{Dir39} vectors, linear forms, scalar and 
outer products (the latter is an example of an operator) are denoted 
through $\vert u \rangle$, $\langle f \vert$, $\langle f \vert u 
\rangle$, and $\vert u \rangle \langle f \vert$, respectively, with 
clear distinction between vectors and forms. In quantum mechanics the 
evolution is \emph{linear}, so there is no need for branching. 
Similarly, the stability function $R(z) = 1 + z \mkern2mu \VEC{b} 
\mkern0.5mu (\MAT{I} - z \mkern2mu \MAT{A})^{-1} \VEC{1} = 1 + z 
\mkern2mu \VEC{b} \mkern0.5mu \VEC{1} + z^{\mkern1mu 2} \VEC{b} 
\mkern2mu \MAT{A} \VEC{1} + z^{\mkern1mu 3} \VEC{b} \mkern2mu \MAT{A}^2 
\VEC{1} + {...}} \def\DS{\displaystyle = 1 + \sum_{n = 0}^\infty z^{\mkern1mu n + 1} \VEC{b} 
\mkern2mu \MAT{A}^n \VEC{1}$ (see, \EG, \citep[s.~238]{But16}, 
\citep[s.~5.3]{But21}, \citep[s.~4.4]{AsPe98}), which describes the 
behavior of a Runge--Kutta method applied to a \emph{linear} equation 
$\d \mkern0.5mu x / \d \mkern1mu t = \lambda \mkern1mu x$ with $z = \lambda 
\mkern1mu h$ (\IE, $\partial^2 \mkern-1mu f_{\mkern1mu i} / \partial 
x_{\mkern-1mu j} \mkern1.5mu \partial x_k = 0$), contains only the terms 
that correspond to rooting trees with no branching.

For the system $\d \mkern0.5mu \VEC{x} / \d \mkern1mu t = \VEC{f}(t ,\mkern2mu} \def\DWS{\mkern3.5mu \VEC{x})$ 
its integral form

\vspace{-11.3pt}

 \begin{gather*}
   \VEC{x}(t + h) = \, \VEC{x}( \mkern0.25mu t \mkern-0.5mu ) \, + \int\limits_t^{t + h} \d \mkern1mu t' 
\DWS \VEC{f} \bigl( t' \! ,\mkern2mu} \def\DWS{\mkern3.5mu \VEC{x}(t') \bigr) = \, \VEC{x}( \mkern0.25mu t \mkern-0.5mu ) \, + 
\, h \int\limits_0^1 \d \theta \DWS \VEC{f} \bigl( t + \theta h ,\mkern2mu} \def\DWS{\mkern3.5mu 
\VEC{x}(t + \theta h) \bigr)
 \end{gather*}

\vspace{-3.3pt}

\noindent can be interpreted as the application of an idealized Runge--Kutta 
method with stages being indexed by an interval $[ \mkern1mu 0, 1 \mkern0.5mu ]${} (Plato's Form 
of a Runge--Kutta method, it is called ``the Picard method'' in 
\citep[p.~97]{But72}, \citep[pp.~153 and 174]{But21}). In the upper half 
of the following table:

\vspace{-10.45pt}

 \begin{gather*}
   \begin{tabular}{r|l}
     $\VEC{1}$ & function $1(\theta) \equiv 1$ \\
     $\MAT{A} \vphantom{\vert_\vert^\vert}$ & operator $u(\theta) \mapsto \int_0^\theta \d 
\theta' \DWS u(\theta')$ \\
     $\VEC{b}$ & functional $u(\theta) \mapsto \int_0^1 \d \theta \DWS
u(\theta)$ \\
     $\VEC{u} . \VEC{v} \vphantom{\vert_{\vert_{\vert_\vert}}} \vphantom{\vert_\vert^\vert}$ & 
point-wise product $u(\theta) \mkern0.5mu v(\theta)$ \\
     \hline
     $\VEC{c} \vphantom{\vert^{\vert^{\vert^{-}}}}$ & $\theta$ \\
     $\VEC{q}_n \vphantom{\vert_\vert} $ & $0$ \\
     $\VEC{\Phi}(\textrm{t}) \vphantom{\vert_\vert^\vert}$ & $\vert \textrm{t} \vert \mkern2mu 
\theta^{\mkern1mu \vert \textrm{t} \vert - 1} \mkern-2mu / \mkern1mu 
\textrm{t}!$ \\
     ${\VEC{b} \mkern2mu \VEC{\Phi}}(\textrm{t}) \vphantom{\vert^{\vert^{\vert}}}$ & $1 / 
\mkern1mu \textrm{t}!$
   \end{tabular} \end{gather*}

\vspace{-2.35pt}

\noindent the process$\mkern1.25mu$/$\mkern0.5mu$result of the substitution of 
entries in the left column with what is on the right will be called a 
\emph{transfiguration}. It can be viewed in two ways: 1) considering 
quantities or statements in the idealized method; and 2) checking for 
desired properties of a Runge--Kutta method, such as simplifying 
assumptions, see, \EG, \citep[p.~52]{But64a}, \citep[s.~321]{But16}, 
\citep[pp.~175, 182, and 208]{HNW93}. Here $\vert \textrm{t} \vert$ is 
the order of tree $\textrm{t}$, \IE, the number of vertices in 
$\textrm{t}$. The factorial $\textrm{t}!$ is recursively defined as 
$\begin{picture}(4.2,6)(0,0) \put(3.35,2.8){\makebox(0,0){\circle*{3}}} 
\end{picture} ! = 1$ and if $\mathrm{t} = [ \mkern1mu \textrm{t}_1 
\mkern2mu \textrm{t}_2 \mkern2mu {...}} \def\DS{\displaystyle \mkern2mu \textrm{t}_n \mkern1mu 
]$, then $\mathrm{t}! = \vert \textrm{t} \vert \mkern1mu 
\prod_{\mkern1mu m = 1}^{\mkern1mu n} \bigl( \textrm{t}_m \bigr) 
\mkern-0.5mu ! \mkern1mu$.

A Runge--Kutta method of order $4$ should satisfy the order conditions 
${\VEC{b} \mkern2mu \VEC{\Phi}}(\textrm{t}) = 1 / \textrm{t}!$ for all rooted trees $\textrm{t}$ 
with $\vert \textrm{t} \vert \le 4$:

\vspace{-10pt}

 \begin{gather*}
   \begin{picture}(312,150)(0,0)
     \put(16.9,131){\includegraphics[height=19pt]{leaf.pdf}}
     \put(18.55,107){\includegraphics[height=22pt]{root.pdf}}
     \put(24,130){\makebox(0,0){\circle*{3}}}
     \put(40,130){\makebox(0,0){\scalebox{1.4}{$=$}}}
     \put(49,131.4){\makebox(0,0)[l]{\scalebox{1.25}{$1$}}}
     \put(94.9,131){\includegraphics[height=19pt]{leaf.pdf}}
     \put(96.55,90){\includegraphics[height=22pt]{root.pdf}}
     \put(102,121.5){\makebox(0,0){\rotatebox{90}{\rule[3pt]{17pt}{1.44pt}}}}
     \put(102,130){\makebox(0,0){\circle*{3}}}
     \put(102,113){\makebox(0,0){\circle*{3}}}
     \put(118,113){\makebox(0,0){\scalebox{1.4}{$=$}}}
     \put(127,114){\makebox(0,0)[l]{$\DS\frac{1}{2}$}}
     \put(165.2,126.5){\rotatebox{22.5}{\includegraphics[height=19pt]{leaf.pdf}}}
     \put(190.2,130.8){\rotatebox{-22.5}{\includegraphics[height=19pt]{leaf.pdf}}}
     \put(182.55,90){\includegraphics[height=22pt]{root.pdf}}
     \put(183.75,120.35){\makebox(0,0){\rotatebox{120}{\rule[2.78pt]{17pt}{1.44pt}}}}
     \put(192.25,120.35){\makebox(0,0){\rotatebox{60}{\rule[2.6pt]{17pt}{1.44pt}}}}
     \put(179.5,127.7){\makebox(0,0){\circle*{3}}}
     \put(196.5,127.7){\makebox(0,0){\circle*{3}}}
     \put(188,113){\makebox(0,0){\circle*{3}}}
     \put(208,113){\makebox(0,0){\scalebox{1.4}{$=$}}}
     \put(217,114){\makebox(0,0)[l]{$\DS\frac{1}{3}$}}
     \put(257.9,131){\includegraphics[height=19pt]{leaf.pdf}}
     \put(259.55,73){\includegraphics[height=22pt]{root.pdf}}
     \put(265,121.5){\makebox(0,0){\rotatebox{90}{\rule[3pt]{17pt}{1.44pt}}}}
     \put(265,104.5){\makebox(0,0){\rotatebox{90}{\rule[3pt]{17pt}{1.44pt}}}}
     \put(265,130){\makebox(0,0){\circle*{3}}}
     \put(265,113){\makebox(0,0){\circle*{3}}}
     \put(265,96){\makebox(0,0){\circle*{3}}}
     \put(281,96){\makebox(0,0){\scalebox{1.4}{$=$}}}
     \put(290,97){\makebox(0,0)[l]{$\DS\frac{1}{6}$}}
     \put(0.9,75){\includegraphics[height=19pt]{leaf.pdf}}
     \put(2.55,0){\includegraphics[height=22pt]{root.pdf}}
     \put(8,65.5){\makebox(0,0){\rotatebox{90}{\rule[3pt]{17pt}{1.44pt}}}}
     \put(8,48.5){\makebox(0,0){\rotatebox{90}{\rule[3pt]{17pt}{1.44pt}}}}
     \put(8,31.5){\makebox(0,0){\rotatebox{90}{\rule[3pt]{17pt}{1.44pt}}}}
     \put(8,74){\makebox(0,0){\circle*{3}}}
     \put(8,57){\makebox(0,0){\circle*{3}}}
     \put(8,40){\makebox(0,0){\circle*{3}}}
     \put(8,23){\makebox(0,0){\circle*{3}}}
     \put(24,23){\makebox(0,0){\scalebox{1.4}{$=$}}}
     \put(33,24){\makebox(0,0)[l]{$\DS\frac{1}{24}$}}
     \put(72.2,53.5){\rotatebox{22.5}{\includegraphics[height=19pt]{leaf.pdf}}}
     \put(97.2,57.8){\rotatebox{-22.5}{\includegraphics[height=19pt]{leaf.pdf}}}
     \put(89.55,0){\includegraphics[height=22pt]{root.pdf}}
     \put(95,31.5){\makebox(0,0){\rotatebox{90}{\rule[3pt]{17pt}{1.44pt}}}}
     \put(90.75,47.35){\makebox(0,0){\rotatebox{120}{\rule[2.78pt]{17pt}{1.44pt}}}}
     \put(99.25,47.35){\makebox(0,0){\rotatebox{60}{\rule[2.6pt]{17pt}{1.44pt}}}}
     \put(86.5,54.7){\makebox(0,0){\circle*{3}}}
     \put(103.5,54.7){\makebox(0,0){\circle*{3}}}
     \put(95,40){\makebox(0,0){\circle*{3}}}
     \put(95,23){\makebox(0,0){\circle*{3}}}
     \put(114,23){\makebox(0,0){\scalebox{1.4}{$=$}}}
     \put(123,24){\makebox(0,0)[l]{$\DS\frac{1}{12}$}}
     \put(162.2,36.5){\rotatebox{22.5}{\includegraphics[height=19pt]{leaf.pdf}}}
     \put(186.4,55.7){\includegraphics[height=19pt]{leaf.pdf}}
     \put(179.55,0){\includegraphics[height=22pt]{root.pdf}}
     \put(193.5,46.2){\makebox(0,0){\rotatebox{90}{\rule[3pt]{17pt}{1.44pt}}}}
     \put(180.75,30.35){\makebox(0,0){\rotatebox{120}{\rule[2.78pt]{17pt}{1.44pt}}}}
     \put(189.25,30.35){\makebox(0,0){\rotatebox{60}{\rule[2.6pt]{17pt}{1.44pt}}}}
     \put(193.5,54.7){\makebox(0,0){\circle*{3}}}
     \put(176.5,37.7){\makebox(0,0){\circle*{3}}}
     \put(193.5,37.7){\makebox(0,0){\circle*{3}}}
     \put(185,23){\makebox(0,0){\circle*{3}}}
     \put(205,23){\makebox(0,0){\scalebox{1.4}{$=$}}}
     \put(214,24){\makebox(0,0)[l]{$\DS\frac{1}{8}$}}
     \put(242.7,35.8){\rotatebox{22.5}{\includegraphics[height=19pt]{leaf.pdf}}}
     \put(263.9,38){\includegraphics[height=19pt]{leaf.pdf}}
     \put(278.7,40.1){\rotatebox{-22.5}{\includegraphics[height=19pt]{leaf.pdf}}}
     \put(265.55,0){\includegraphics[height=22pt]{root.pdf}}
     \put(264,30){\makebox(0,0){\rotatebox{135}{\rule[2.12pt]{20pt}{1.44pt}}}}
     \put(271,30){\makebox(0,0){\rotatebox{90}{\rule[3pt]{14pt}{1.44pt}}}}
     \put(278,30){\makebox(0,0){\rotatebox{45}{\rule[2.12pt]{20pt}{1.44pt}}}}
     \put(257,37){\makebox(0,0){\circle*{3}}}
     \put(271,37){\makebox(0,0){\circle*{3}}}
     \put(285,37){\makebox(0,0){\circle*{3}}}
     \put(271,23){\makebox(0,0){\circle*{3}}}
     \put(296,23){\makebox(0,0){\scalebox{1.4}{$=$}}}
     \put(305,24){\makebox(0,0)[l]{$\DS\frac{1}{4}$}}
   \end{picture} \end{gather*}

\vspace{-2.75pt}

\noindent or $\VEC{b} \mkern0.5mu \VEC{1} = 1$, $\VEC{b} \mkern1mu \VEC{c} = 
\frac12$, $\VEC{b} \mkern1mu \VEC{c}^2 = \frac13$, $\VEC{b} \mkern2mu 
\MAT{A} \mkern0.5mu \VEC{c} = \frac{1}{6}$, $\VEC{b} \mkern1mu \VEC{c}^3 
= \frac{1}{4}$, $\VEC{b} \bigl( \VEC{c} . (\MAT{A} \mkern0.5mu \VEC{c}) 
\bigr) = \frac{1}{8}$, $\VEC{b} \mkern2mu \MAT{A} \mkern0.5mu \VEC{c}^2 
= \frac{1}{12}$, $\VEC{b} \mkern2mu \MAT{A}^2 \VEC{c} = \frac{1}{24}$. 
Consider the diagram $\VEC{b} \mkern2mu \MAT{A} \mkern0.5mu \VEC{c}$ in 
the top right statement $\VEC{b} \mkern2mu \MAT{A} \mkern0.5mu \VEC{c} = 
\frac{1}{6}$. If the root or the leaf are removed, then the diagram 
becomes the vector $\MAT{A} \mkern0.5mu \VEC{c} = \MAT{A}^2 \VEC{1}$ or 
the row vector $\VEC{b} \mkern2mu \MAT{A}^2$, respectively. If both the 
root and the leaf are removed, then the diagram becomes the matrix 
$\MAT{A}^2$. The transfiguration of, \EG, the statement $\VEC{b} 
\mkern2mu \MAT{A} \mkern0.5mu \VEC{c}^2 = \frac{1}{12}$ is

 \vspace{-12pt}

 \begin{align*}
  & \int\limits_0^1 \d \theta \int\limits_0^\theta \d \theta' \Biggl( 
\int\limits_0^{\theta'} \d \theta'' \Biggr) \Biggl( 
\int\limits_0^{\theta'} \d \theta''' \Biggr) = \int\limits_0^1 \d \theta 
\int\limits_0^\theta \d \theta' \DWS \theta'^2 = \int\limits_0^1 \d 
\theta \DWS \frac{\theta^3}{3} = \frac{1}{12} \\
  & \mkern20mu \VEC{b} \mkern30mu \MAT{A} \mkern15mu \bigl( \mkern9mu 
(\MAT{A} \VEC{1}) \mkern18mu . \mkern18mu (\MAT{A} \VEC{1}) \mkern10mu 
\bigr) \mkern3mu = \mkern17mu \VEC{b} \mkern29mu \MAT{A} \mkern20mu 
\VEC{c}^2
 \end{align*}

 \vspace{-4pt}

Let $D \bigl( \VEC{\Phi}(\textrm{t}_1) ,\mkern2mu} \def\DWS{\mkern3.5mu \VEC{\Phi}(\textrm{t}_2) 
\bigr)$ denote the following property (see \citep[eq.~(2.5)]{AuCh98a}):

\vspace{-14pt}

 \begin{gather*}
   \begin{picture}(260,81)(-2,0)
     \put(-0.65,18){\rotatebox{30}{\includegraphics[height=49pt]{birch.pdf}}}
     \put(34,39.7){\rotatebox{-22.5}{\includegraphics[height=49pt]{larch.pdf}}}
     \put(29.55,0){\includegraphics[height=22pt]{root.pdf}}
     \put(40.28,29.72){\makebox(0,0){\rotatebox{45}{\rule{17pt}{1.44pt}}}}
     \put(47,35){\makebox(0,0){\circle*{3}}}
     \put(35,23){\makebox(0,0){\circle*{3}}}
     \put(-0.8,59){\makebox(0,0){$\textrm{t}_1$}}
     \put(46.5,77){\makebox(0,0){$\textrm{t}_2$}}
     \put(85,22.5){\makebox(0,0){\scalebox{2}{$+$}}}
     \put(92.85,31.35){\rotatebox{22.5}{\includegraphics[height=49pt]{birch.pdf}}}
     \put(122.7,29.3){\rotatebox{-30}{\includegraphics[height=49pt]{larch.pdf}}}
     \put(129.55,0){\includegraphics[height=22pt]{root.pdf}}
     \put(126.72,29.72){\makebox(0,0){\rotatebox{-45}{\rule{17pt}{1.44pt}}}}
     \put(123,35){\makebox(0,0){\circle*{3}}}
     \put(135,23){\makebox(0,0){\circle*{3}}}
     \put(182,21.5){\makebox(0,0){\scalebox{2}{$=$}}}
     \put(204.55,23.45){\includegraphics[height=49pt]{birch.pdf}}
     \put(211.55,0){\includegraphics[height=22pt]{root.pdf}}
     \put(217,23){\makebox(0,0){\circle*{3}}}
     \put(230,23){\includegraphics[height=49pt]{larch.pdf}}
     \put(238.55,0){\includegraphics[height=22pt]{root.pdf}}
     \put(244,23){\makebox(0,0){\circle*{3}}}
   \end{picture} \\
   \VEC{b} \bigl(\VEC{\Phi}(\textrm{t}_1) \mkern1mu . \mkern-1mu \bigl( 
\MAT{A} \VEC{\Phi}(\textrm{t}_2) \bigr) \bigr) + \VEC{b} \bigl( \bigl( 
\MAT{A} \VEC{\Phi}(\textrm{t}_1) \bigr) . \mkern1mu 
\VEC{\Phi}(\textrm{t}_2) \bigr) = \bigl( \VEC{b} \mkern2mu 
\VEC{\Phi}(\textrm{t}_1) \bigr) \bigl( \VEC{b} \mkern2mu 
\VEC{\Phi}(\textrm{t}_2) \bigr)
 \end{gather*}

\vspace{-4pt}

\noindent whose transfiguration or continuous analog $D(u ,\mkern2mu} \def\DWS{\mkern3.5mu \mkern-1mu v)$ 
would be

\vspace{-13pt}

\begin{gather*}
   \int\limits_0^1 \d \theta \DWS u(\theta) \int\limits_0^\theta \d 
\theta' \DWS v(\theta') \; + \, \int\limits_0^1 \d \theta \Biggl( 
\int\limits_0^\theta \d \theta' \DWS u(\theta') \Biggr) v(\theta) \, = 
\, \Biggl( \int\limits_0^1 \d \theta \DWS u(\theta) \Biggr) \Biggl( 
\int\limits_0^1 \d \theta \DWS v(\theta) \Biggr) \\[3pt]
  \includegraphics{Duv.pdf} \end{gather*}

\vspace{-10pt}

\noindent where, assuming it is $u(\theta_1) \mkern0.5mu v(\theta_2)$ being 
integrated, the shading shows the regions of integration in the 
$(\theta_1 ,\mkern2mu} \def\DWS{\mkern3.5mu \theta_2)$-plane. Another way to interpret $D(u ,\mkern2mu} \def\DWS{\mkern3.5mu 
\mkern-1mu v)$ is as the integration by parts formula. Consider the 
following bilinear form $M$:

\vspace{-12pt}

 \begin{gather*}
   \begin{picture}(240,37)(0,0)
     \qbezier(1.5,35)(13.5,3)(25.5,35)
     \put(20.5,18){\makebox(0,0){\color{white}\circle*{14}}}
     \put(13.6,18){\makebox(0,0){$M$}}
     \put(20.5,18){\makebox(0,0){\circle{14}}}
     \put(3,35){\makebox(0,0){\circle*{3}}}
     \put(27,35){\makebox(0,0){\circle*{3}}}
     \put(48,29){\makebox(0,0){\scalebox{1.7}{$=$}}}
     \put(79.55,0){\includegraphics[height=22pt]{root.pdf}}
     \put(76.72,29.72){\makebox(0,0){\rotatebox{135}{\rule{17pt}{1.44pt}}}}
     \put(76.72,29.72){\makebox(0,0){\rotatebox{135}{\color{white}\rule{17pt}{0.72pt}}}}
     \put(90.28,29.72){\makebox(0,0){\rotatebox{45}{\rule{17pt}{1.44pt}}}}
     \put(73,35){\makebox(0,0){\circle*{3}}}
     \put(97,35){\makebox(0,0){\circle*{3}}}
     \put(85,23){\makebox(0,0){\circle*{3}}}
     \put(119,29){\makebox(0,0){\scalebox{1.7}{$+$}}}
     \put(147.55,0){\includegraphics[height=22pt]{root.pdf}}
     \put(144.72,29.72){\makebox(0,0){\rotatebox{135}{\rule{17pt}{1.44pt}}}}
     \put(158.28,29.72){\makebox(0,0){\rotatebox{45}{\rule{17pt}{1.44pt}}}}
     \put(158.28,29.72){\makebox(0,0){\rotatebox{45}{\color{white}\rule{17pt}{0.72pt}}}}
     \put(141,35){\makebox(0,0){\circle*{3}}}
     \put(165,35){\makebox(0,0){\circle*{3}}}
     \put(153,23){\makebox(0,0){\circle*{3}}}
     \put(186,29){\makebox(0,0){\scalebox{1.7}{$-$}}}
     \put(207.55,12){\includegraphics[height=22pt]{root.pdf}}
     \put(213,35){\makebox(0,0){\circle*{3}}}
     \put(231.55,12){\includegraphics[height=22pt]{root.pdf}}
     \put(237,35){\makebox(0,0){\circle*{3}}}
   \end{picture} \end{gather*}

\vspace{-7pt}

\noindent whose matrix $\MAT{M} = \bigl[ m_{\mkern1mu i \mkern-1mu j} \bigr]$, 
with $m_{\mkern1mu i \mkern-1mu j} = b_{i \mkern1mu} {a_{\mkern0.85mu i \mkern-0.85mu j}} + 
b_{\mkern-1.2mu j \mkern1.2mu} a_{\mkern-1.2mu j \mkern1.2mu i} - b_{i 
\mkern1mu} b_{\mkern-1mu j \mkern1mu}$, is symmetric. The property 
$D(\VEC{u} ,\mkern2mu} \def\DWS{\mkern3.5mu \mkern-1mu \VEC{v})$ can be written as $M(\VEC{u} ,\mkern2mu} \def\DWS{\mkern3.5mu 
\mkern-1mu \VEC{v}) = \VEC{u}^{\textrm{T}} \MAT{M} \mkern1mu \VEC{v} = 
0$.

An order $p$ Runge--Kutta method is said to be pseudo-symplectic of 
order $(p ,\mkern2mu} \def\DWS{\mkern3.5mu q)$ if $D \bigl( \VEC{\Phi}(\textrm{t}_1) ,\mkern2mu} \def\DWS{\mkern3.5mu 
\VEC{\Phi}(\textrm{t}_2) \bigr)$ for all trees $\mathrm{t}_1$ and 
$\mathrm{t}_2$ such that $\vert \mathrm{t}_1 \vert + \vert \mathrm{t}_2 
\vert \le q$ \citep[cor.~2.2 and eq.~(2.7)]{AuCh98a}. $B$-series 
coefficients, $\mathrm{t} \mapsto {\VEC{b} \mkern2mu \VEC{\Phi}}(\textrm{t})$, for a method of 
order $(p ,\mkern2mu} \def\DWS{\mkern3.5mu q)$ are in the coset $\bigl( \mathrm{E} + \mathsf{O}_{p + 
1} \bigr)$ and in the subgroup $\bigl( \cap_{\mkern1mu r = 1}^{\mkern1mu 
q - 1} \mathsf{D}_{\mkern1.25mu r ,\mkern2mu} \def\DWS{\mkern3.5mu q - r} \bigr)$ of the Butcher 
group $\mathsf{B}$ \citep[pp.~166 and 173]{But21}.

The property $\bigl( D \bigl( \VEC{\Phi} ,\mkern2mu} \def\DWS{\mkern3.5mu \VEC{u} \bigr)$ for all 
$\VEC{u} \in \pmb{\mathbf{R}}^s \bigr)$ will be shortened as $\VEC{D} 
\bigl( \VEC{\Phi} \bigr)$ (see \citep[def.~3.1]{AuCh98a}). It is 
equivalent to the linear form $\VEC{u} \mapsto M \bigl( \VEC{\Phi} ,\mkern2mu} \def\DWS{\mkern3.5mu 
\VEC{u} \bigr)$ being a zero function, or $\MAT{M} \mkern0.5mu 
\VEC{\Phi} = \VEC{0}$. The \NTH{\smash{j{\mkern1.2mu}}} component of 
$\VEC{D} \bigl( \VEC{\Phi} \bigr)$, the property $D \bigl( \VEC{\Phi} 
,\mkern2mu} \def\DWS{\mkern3.5mu \VEC{e}_{\mkern-1.5mu j} \bigr)$, where $\VEC{e}_{\mkern-1.5mu j}$ 
is the \NTH{\smash{j{\mkern1.2mu}}} vector of the standard basis of 
$\pmb{\mathbf{R}}^s$, will be denoted by $D_{\mkern-1.5mu j} \bigl( 
\VEC{\Phi} \bigr)$. Assuming that $\VEC{b} \mkern0.5mu \VEC{1} = 1$ and 
$\VEC{b} \mkern1mu \VEC{c} = \smash{\frac{1}{2}}$, here are the diagrams 
for $\VEC{D}(\VEC{1} \mkern-0.75mu )$ and $\VEC{D}(\VEC{c})$:

 \vspace{-13pt}

 \begin{gather*}
   \begin{picture}(318,162)(0,0)
     \put(130.55,113.45){\includegraphics[height=49pt]{birch.pdf}}
     \put(137.53,75){\includegraphics[height=22pt]{root.pdf}}
     \put(143,104.5){\makebox(0,0){\rotatebox{90}{\rule[2.88pt]{17pt}{1.44pt}}}}
     \put(143,113){\makebox(0,0){\circle*{3}}}
     \put(143,96){\makebox(0,0){\circle*{3}}}
     \put(172,95.5){\makebox(0,0){\scalebox{2}{$+$}}}
     \put(187.7,106.8){\rotatebox{22.5}{\includegraphics[height=19pt]{leaf.pdf}}}
     \put(203.3,101.9){\rotatebox{-30}{\includegraphics[height=49pt]{birch.pdf}}}
     \put(208.55,73){\includegraphics[height=22pt]{root.pdf}}
     \put(205.72,102.72){\makebox(0,0){\rotatebox{-45}{\rule{17pt}{1.44pt}}}}
     \put(202,108){\makebox(0,0){\circle*{3}}}
     \put(214,96){\makebox(0,0){\circle*{3}}}
     \put(255,94.5){\makebox(0,0){\scalebox{2}{$=$}}}
     \put(273.55,96.45){\includegraphics[height=49pt]{birch.pdf}}
     \put(280.55,73){\includegraphics[height=22pt]{root.pdf}}
     \put(286,96){\makebox(0,0){\circle*{3}}}
     \put(32,113){\makebox(0,0)[lb]{$\DS {\VEC{b} \mkern1.5mu \MAT{A} 
\mkern0.5mu \VEC{u}} + {\VEC{b} \mkern0.5mu (\VEC{c} . \mkern1mu 
\VEC{u})} = {\VEC{b} \mkern1mu \VEC{u}}$}}
     \put(-3.3,33.8){\rotatebox{22.5}{\includegraphics[height=19pt]{leaf.pdf}}}
     \put(17.5,37.7){\includegraphics[height=49pt]{larch.pdf}}
     \put(17.55,0){\includegraphics[height=22pt]{root.pdf}} 
     \put(14.72,29.72){\makebox(0,0){\rotatebox{-45}{\rule{17pt}{1.44pt}}}} 
     \put(27.25,30.35){\makebox(0,0){\rotatebox{60}{\rule[2.6pt]{17pt}{1.44pt}}}} 
     \put(11,35){\makebox(0,0){\circle*{3}}}
     \put(31.5,37.7){\makebox(0,0){\circle*{3}}} 
     \put(23,23){\makebox(0,0){\circle*{3}}} 
     \put(55,22.5){\makebox(0,0){\scalebox{1.8}{$+$}}} 
     \put(69.5,52.4){\includegraphics[height=19pt]{leaf.pdf}}
     \put(80.7,29.3){\rotatebox{-30}{\includegraphics[height=49pt]{larch.pdf}}}
     \put(87.55,0){\includegraphics[height=22pt]{root.pdf}} 
     \put(78.54,44){\makebox(0,0){\rotatebox{105}{\rule[2.78pt]{17pt}{1.44pt}}}}
     \put(84.72,29.72){\makebox(0,0){\rotatebox{-45}{\rule{17pt}{1.44pt}}}} 
     \put(76.6,51.4){\makebox(0,0){\circle*{3}}}
     \put(81,35){\makebox(0,0){\circle*{3}}} 
     \put(93,23){\makebox(0,0){\circle*{3}}} 
     \put(131,21.5){\makebox(0,0){\scalebox{1.8}{$=$}}} 
     \put(147,23){\makebox(0,0)[l]{$\DS\frac{1}{2}$}}
     \put(155,23){\includegraphics[height=49pt]{larch.pdf}}
     \put(163.55,0){\includegraphics[height=22pt]{root.pdf}} 
     \put(169,23){\makebox(0,0){\circle*{3}}} 
     \put(318,40){\makebox(0,0)[rb]{$\DS {\VEC{b} \bigl( \VEC{c} . (\MAT{A} 
\mkern0.5mu \VEC{u}) \mkern-1mu \bigr)} + {\VEC{b} \bigl( \mkern-1mu 
(\MAT{A} \mkern0.5mu \VEC{c}) . \mkern1mu \VEC{u} \bigr)} = 
{{\textstyle} \def\d{{\rm d}} \def\PHI{\varphi\frac12} \mkern0.5mu \VEC{b} \mkern1mu \VEC{u}}$}} \end{picture} 
\end{gather*}

 \vspace{-7pt}

\noindent The property $\VEC{D}(\VEC{1} \mkern-0.75mu )$ coincides with $D(1)$ 
\citep[p.~52]{But64a}, \citep[fig.~321(ii)]{But16}, \citep[pp.~175 and 
208]{HNW93}, \citep[p.~193]{But21}. If $C(2)$ and $D(1)$ are satisfied, 
then $\VEC{D}(\VEC{c})$ is equivalent to $D(2)$. For an explicit 
Runge--Kutta method that satisfies the explicit version $\bigl( (\MAT{A} 
\mkern1mu \VEC{c})_i \ne c_i^2 / 2 \bigr) \Rightarrow \bigl( (i = 2) 
\land (b_2 = 0) \bigr)$ of the property $C(2)$, it is impossible to have 
both $\VEC{D}(\VEC{1} \mkern-0.75mu )$ and $\VEC{D}(\VEC{c})$ \citep[prop.~3.6]{AuCh98a}.

Let a function $u(\theta)$ be called $[ \mkern1mu 0, 1 \mkern0.5mu ]${}-even or $[ \mkern1mu 0, 1 \mkern0.5mu ]${}-odd, if 
$u(1 - \theta) = u(\theta)$ or $u(1 - \theta) = -u(\theta)$, 
respectively. The following parity properties are true:

 \vspace{-14pt}

 \begin{gather} \begin{array}{l}
   {\DS \mbox{if $u(\theta)$ is $[ \mkern1mu 0, 1 \mkern0.5mu ]${}-odd, then $U(\theta) = 
\smash{\int_0^\theta \d \theta' \, u(\theta')}$ is $[ \mkern1mu 0, 1 \mkern0.5mu ]${}-even}} 
\\[6pt]
   {\DS \mbox{if $u(\theta)$ is $[ \mkern1mu 0, 1 \mkern0.5mu ]${}-even and $\smash{\int_0^1 \d 
\theta \, u(\theta) = 0}$, then $U(\theta) = \smash{\int_0^\theta \d 
\theta' \, u(\theta')}$ is $[ \mkern1mu 0, 1 \mkern0.5mu ]${}-odd}}
  \end{array} \label{even_odd} \end{gather}

 \vspace{-9pt}

\section{Family of pseudo-symplectic methods} \label{family}

\hskip\parindent{}A pseudo-symplectic Runge--Kutta method of order $(4 
,\mkern2mu} \def\DWS{\mkern3.5mu 8)$ will be searched for within $8$-stage explicit methods with 
$c_4 = c_5 = \frac12$, $c_6 = 1 - c_3$, $c_7 = 1 - c_2$, $c_8 = 1$, $b_6 
= b_3$, $b_7 = b_2$, and $b_8 = b_1$. Vectors $\VEC{x}$ and $\VEC{y}$ of 
$8$ components are going to be called ``even'' and ``odd'' if $x_6 = 
x_3$, $x_7 = x_2$, $x_8 = x_1$, and $y_1 + y_8 = y_2 + y_7 = y_3 + y_6 = 
y_4 = y_5 = 0$, respectively. Let $\VEC{p}_1 = 2 \VEC{c} - \VEC{1}$, 
$\VEC{p}_2 = 6 \mkern1.5mu \VEC{c}^2 - 6 \mkern1.5mu \VEC{c} + \VEC{1}$, 
and $\VEC{p}_3 = 20 \mkern1.5mu \VEC{c}^3 - 30 \mkern1.5mu \VEC{c}^2 + 
12 \VEC{c} - \VEC{1}$ be the analogs of shifted Legendre polynonials, 
orthogonal on $[ \mkern1mu 0, 1 \mkern0.5mu ]${} (see, \EG, \citep[s.~342]{But16}, 
\citep[p.~203]{But21}). The vectors $\VEC{1}$, $\VEC{b}^{\textrm{T}}$, 
and $\VEC{p}_2$ are ``even''; while $\VEC{p}_1$ and $\VEC{p}_3$ are 
``odd''. With such a definition of ``even'' and ``odd'' vectors, the 
eq.~(\ref{even_odd}) could be viewed as the transfiguration of the 
statements

 \vspace{-15pt}

 \begin{gather*} \begin{array}{l}
   {\DS \mbox{if $\VEC{u}$ is ``odd'', then $\MAT{A} \mkern0.5mu 
\VEC{u}$ is ``even''}} \\[2pt]
   {\DS \mbox{if $\VEC{u}$ is ``even'' and $\VEC{b} \mkern1mu \VEC{u} = 
0$, then $\MAT{A} \mkern0.5mu \VEC{u}$ is ``odd''}}
  \end{array} \end{gather*}

 \vspace{-7pt}

\noindent To imitate these parity properties, it will be demanded that 
$\MAT{A} \VEC{p}_1$ is ``even'' (this implies that the vector $\VEC{q}_1 
= \MAT{A} \VEC{c} - \frac12 \VEC{c}^2 = \frac12 \MAT{A} \VEC{p}_1 + 
\frac12 \VEC{c} (\VEC{1} - \VEC{c})$ is ``even'' too), and that $\MAT{A} 
\VEC{p}_2$ and $\MAT{A} \mkern0.5mu \VEC{q}_1$ are ``odd''. These 
demands can be viewed as simplifying assumptions, as they increase the 
redundancy in the order conditions (for example, the statements $\VEC{b} 
\VEC{1} = 1$ and $\VEC{b} \VEC{p}_2 = \VEC{b} \mkern1mu \VEC{q}_1 = 0$ 
imply that a method is of order $4$). Also the properties $\VEC{D}(\VEC{1} \mkern-0.75mu )$ 
and $\VEC{D}(\VEC{c})$ (or an easier to exploit $\VEC{D}(\VEC{p}_1 \mkern-0.5mu )$) will be 
assumed.\footnote{~The derivation was done in interaction with computer 
algebra system Wolfram Mathematica~8.0, mainly using commands 
\scalebox{0.725}[0.825]{\scriptsize\tt\bfseries Solve} to symbolically 
solve linear equations, \scalebox{0.725}[0.825]{\scriptsize\tt\bfseries 
Simplify}, and \scalebox{0.725}[0.825]{\scriptsize\tt\bfseries Factor}. 
Step-by-step description of the process (of course, this could be done 
in many ways) allows for an easy reproduction.}

The vectors $\VEC{p}_3 - \VEC{p}_1$, $\VEC{p}_1 \mkern-0.5mu . 
\mkern1.75mu \VEC{q}_1$, $\MAT{A} \VEC{p}_2$, and $\MAT{A} \mkern0.5mu 
\VEC{q}_1$ are ``odd'' and have the \FIRST{} and the \NTH{8} components 
being equal to $0$. Thus, they lie in the span of two vectors: 
$\VEC{v}_2 = \bigl[ \; 0 ~~~ 1 ~~~ 0 ~~~ 0 ~~~ 0 ~~~ 0 ~~ -\!1 ~~~ 0 \; 
\bigr]{}^{\textrm{T}}$ and $\VEC{v}_3 = \bigl[ \; 0 ~~~ 0 ~~~ 1 ~~~ 0 
~~~ 0 ~~ -\!1 ~~~ 0 ~~~ 0 \; \bigr]{}^{\textrm{T}}$. To satisfy the 
order conditions $D \bigl( \VEC{\Phi}(\textrm{t}_1), 
\VEC{\Phi}(\textrm{t}_2) \bigr)$ with $\vert \textrm{t}_1 \vert = \vert 
\textrm{t}_2 \vert = 4$, it is necessary and sufficient that $\VEC{b} 
\bigl( \VEC{v}_2 . (\MAT{A} \VEC{v}_2) \bigr) = -a_{72} b_2 = 0$, 
$\VEC{b} \bigl( \VEC{v}_3 . (\MAT{A} \VEC{v}_3) \bigr) = -a_{63} b_3 = 
0$, and $\VEC{b} \bigl( \VEC{v}_2 . (\MAT{A} \VEC{v}_3) \bigr) + \VEC{b} 
\bigl( (\MAT{A} \VEC{v}_2) . \VEC{v}_3 \bigr) = (a_{76} - a_{73}) b_2 + 
(a_{32} - a_{62}) b_3 = 0$. To support the observance of the parity 
properties and not to lose the flexibility in choosing $b_2$ and $b_3$, 
let $a_{63} = a_{72} = 0$, $a_{62} = a_{32}$, and $a_{76} = a_{73}$.

The conditions $-c_2 = (\MAT{A} \VEC{p}_1 \mkern-0.5mu )_2 = (\MAT{A} \VEC{p}_1 \mkern-0.5mu )_7 = -a_{71}$ and $(\MAT{A} \VEC{p}_1 \mkern-0.5mu )_3 = 
(\MAT{A} \VEC{p}_1 \mkern-0.5mu )_6$ imply $a_{71} = c_2$ and $a_{61} = a_{31}$, respectively. The 
property $D_1{(\VEC{p}_1 \mkern-0.5mu )}$ becomes $a_{81} b_1 = 0$, thus $a_{81} = 0$.

The property $D_2{(\VEC{p}_1 \mkern-0.5mu )}$ implies $a_{82} = c_2 b_2 / b_1$. Both 
$D_7{(\VEC{1} \mkern-0.75mu )}$ and $D_7{(\VEC{p}_1 \mkern-0.5mu )}$ imply $a_{87} = c_2 b_2 / b_1$. Now all 
three $D_3{(\VEC{p}_1 \mkern-0.5mu )} \Leftrightarrow D_6{(\VEC{p}_1 \mkern-0.5mu )}$, $D_8{(\VEC{p}_1 \mkern-0.5mu )}$, and 
$(\MAT{A} \VEC{p}_1 \mkern-0.5mu )_8 = 0$ imply $a_{86} = a_{83}$.

If $b_2 = 0$, then $D_6{(\VEC{1} \mkern-0.75mu )}$ implies $a_{83} = c_3 b_3 / b_1$, and 
then $D_6{(\VEC{p}_1 \mkern-0.5mu )}$ means $b_3 a_{32} c_2 = 0$. Both cases $b_3 = 0$ and 
$a_{32} = 0$ lead to contradictions, thus $b_2 \ne 0$. With $a_{31} = 
c_3 - a_{32}$, the properties $D_6{(\VEC{1} \mkern-0.75mu )}$ and $D_6{(\VEC{p}_1 \mkern-0.5mu )}$ imply 
$a_{83} = a_{31} b_3 / b_1$ and $a_{73} = a_{32} b_3 / b_2$. The Butcher 
tableau \citep[p.~191]{But64b} now looks like

 \vspace{-11pt}

 \begin{gather*}
   \begin{array}{r|ccccccccc}
     0\;\,\vphantom{\vert_\vert} \\
     {c_2}\;\,\vphantom{\vert_\vert^\vert} & \phantom{-}c_2 \\
     {c_3}\;\,\vphantom{\vert_\vert^\vert} & \phantom{-}a_{31} & \phantom{-}a_{32} \\
     \frac{1}{2}\;\,\vphantom{\vert_\vert^\vert} & \phantom{-}a_{41} & \phantom{-}a_{42} & 
\phantom{-}a_{43} \\
     \frac{1}{2}\;\,\vphantom{\vert_\vert^\vert} & \phantom{-}a_{51} & \phantom{-}a_{52} & 
\phantom{-}a_{53} & \phantom{-}a_{54} \\
     {1 - c_3}\;\,\vphantom{\vert_\vert^\vert} & \,\phantom{-}a_{31}\, & \phantom{-}a_{32} & 
\phantom{-}0 & \,\phantom{-}a_{64}\, & \,\phantom{-}a_{65}\, \\
     {1 - c_2}\;\,\vphantom{\vert_\vert^\vert} & \phantom{-}c_2 & \phantom{-}0 & 
\phantom{-}a_{32} \frac{b_3}{b_2} & \phantom{-}a_{74} & 
\phantom{-}a_{75} & \phantom{-}a_{32} \frac{b_3}{b_2} \\
     1\;\,\vphantom{\vert_\vert^\vert}\vphantom{\vert_{\vert_\vert}} & \phantom{-}0 & 
\phantom{-}c_{2} \frac{b_2}{b_1} & \phantom{-}a_{31} \frac{b_3}{b_1} & 
\phantom{-}a_{84} & \phantom{-}a_{85} & \phantom{-}a_{31} 
\frac{b_3}{b_1} & \phantom{-}c_{2} \frac{b_2}{b_1} \\
     \hline
     \vphantom{\vert^{\vert^\vert}} & \phantom{-}b_1 & \phantom{-}b_2 & 
\phantom{-}b_3 & \phantom{-}b_4 & \phantom{-}b_5 & \phantom{-}b_3 & 
\phantom{-}b_2 & \,\phantom{-}b_1
   \end{array} \end{gather*} 

 \vspace{-2pt}

As ${(\MAT{A} \VEC{1} \mkern-0.75mu )}_6 = 1 - c_3$, let $a_{64} = 1 - 2 c_3 - a_{65}$. Then 
$(\MAT{A} \VEC{p}_2)_3 + (\MAT{A} \VEC{p}_2)_6 = 0$ implies $c_3 \ne 
\frac{1}{6}$ and $a_{32} = (6 c_3 - 1) / 24 c_2 (1 - c_2)$. As 
${(\MAT{A} \VEC{1} \mkern-0.75mu )}_4 = \frac12$, let $a_{41} = \frac12 - a_{42} - a_{43}$. Then 
$(\MAT{A} \VEC{p}_2)_4 = 0$ implies $a_{42} = \bigl( 1 - 12 a_{43} c_3 
(1 - c_3) \bigr) / 12 c_2 (1 - c_2)$. Now $(\MAT{A} \mkern1mu \VEC{q}_1 
\mkern-1.5mu )_4 = 0$ is equivalent to $a_{43} = c_2 / \bigl( 6 c_3 - 12 
c_3 (c_3 - c_2) - 1 \bigr)$, with the denominator necessarily being 
non-zero. As ${(\MAT{A} \VEC{1} \mkern-0.75mu )}_5 = \frac12$, let $a_{51} = \frac12 - a_{52} - 
a_{53} - a_{54}$. Then $(\MAT{A} \VEC{p}_2)_5 = 0$ implies $a_{54} = 
\frac13 - 4 a_{52} c_2 (1 - c_2) - 4 a_{53} c_3 (1 - c_3)$.

As $(\MAT{A} \VEC{1})_7 = 1 - c_2$, let $a_{74} = 1 - 2 c_2 - a_{75} - 
b_3 (6 c_3 - 1) / 12 b_2 c_2 (1 - c_2)$. Then $(\MAT{A} \VEC{p}_2)_2 + 
(\MAT{A} \VEC{p}_2)_7 = \bigl( b_3 (6 c_3 - 1) (1 - 2 c_3)^2 - 4 b_2 c_2 
(1 - c_2) (1 - 6 c_2) \bigr) / 8 b_2 c_2 (1 - c_2) = 0$, thus $b_3 = 4 
b_2 c_2 (1 - c_2) (1 - 6 c_2) / (1 - 2 c_3)^2 (6 c_3 - 1)$. The 
coefficient $a_{84}$ is determined from ${(\MAT{A} \VEC{1} \mkern-0.75mu )}_8 = 1$. Then the 
derivative of $b_1 (\MAT{A} \VEC{p}_2)_8$ with respest to $b_1$ is equal 
to $-\frac{1}{2} \ne 0$, so $(\MAT{A} \VEC{p}_2)_8 = 0$ determines $b_1 
= b_2 \bigl( 1 - 12 c_2 + 24 c_2^2 (1 - c_2) - 6 c_3 + 96 c_2 c_3 - 312 
c_2^2 c_3 + 288 c_2^3 c_3 \bigr) / (6 c_3 - 1)$. Now $\MAT{A} \VEC{p}_2 
= \bigl[ \; 0 ~~~ c_2 ~~~ \frac{1 - 2 c_3}{4} ~~~ 0 ~~~ 0 ~~~ 
{-{\mkern-1mu}\frac{1 - 2 c_3}{4}} ~~ {-{\mkern-1mu}c_2} ~~~ 0 \; 
\bigr]{}^{\textrm{T}}$ is ``odd'', also $\VEC{q}_1$ is ``even''. The 
coefficient $a_{85}$ and the weight $b_4$ are determined from 
$D_5(\VEC{1})$ and $\VEC{b} \mkern0.5mu \VEC{1} = 1$, respectively.

Now $\VEC{b} \mkern1mu \VEC{c}^2 = \frac{1}{4} + 2 b_2 (1 - 2 c_3)^2 
\bigl( 18 c_3 (1 - 2 c_2)^2 - 1 - 2 c_2 \bigr) / (6 c_3 - 1)$. Thus, $18 
c_3 (1 - 2 c_2)^2 - 1 - 2 c_2 \ne 0$ and $b_3 = (6 c_3 - 1) / 24 c_2 
\bigl( 18 c_3 (1 - 2 c_2)^2 - 1 - 2 c_2 \bigr)$.

The coefficient $a_{75}$ is determined from $D_5{(\VEC{p}_1 \mkern-0.5mu )}$.

The properties $D_2{(\VEC{1} \mkern-0.75mu )}$, $D_3{(\VEC{1} \mkern-0.75mu )}$, and $(\MAT{A} \mkern0.5mu \VEC{q}_1 \mkern-0.5mu )_5 = 0$ can be 
written as $b_5 a_{52} + \phi_{2b} b_5 = \phi_{20}$, $b_5 a_{53} + 
\phi_{3b} b_5 = \phi_{30}$, and $\phi_2 a_{52} + \phi_3 a_{53} = 
\phi_0$, respectively, where all the coefficients $\phi_{2b}$, 
$\phi_{20}$, $\phi_{3b}$, $\phi_{30}$, $\phi_2$, $\phi_3$, and $\phi_0$ 
are rational functions of $c_2$ and $c_3$. It can be checked that 
$\phi_2 \phi_{2b} + \phi_3 \phi_{3b} + \phi_0 = 0$. The following 
quantity should be equal to zero: $\phi_2 (b_5 a_{52} + \phi_{2b} b_5 - 
\phi_{20}) + \phi_3 (b_5 a_{53} + \phi_{3b} b_5 - \phi_{30}) = b_5 
(\phi_2 a_{52} + \phi_3 a_{53} + \phi_2 \phi_{2b} + \phi_3 \phi_{3b}) - 
\phi_2 \phi_{20} - \phi_3 \phi_{30} = b_5 \bigl( \phi_0 + (-\phi_0) 
\bigr) - \phi_2 \phi_{20} - \phi_3 \phi_{30} = -(\phi_2 \phi_{20} + 
\phi_3 \phi_{30})$. This can be rewritten as

 \vspace{-18pt}

 \begin{gather*}
   \zeta(c_2 ,\mkern2mu} \def\DWS{\mkern3.5mu c_3) = \sum_{m, n} \zeta_{\mkern1mu m \mkern0.5mu n} 
\mkern1.5mu c_2^m \mkern0.5mu (2 c_3)^n = 0 \qquad
   \begin{array}{r|rrrrr}
     \mbox{\footnotesize$5$}\; & \phantom{-}9 & -9 \\
     \mbox{\footnotesize$4$}\; & {\;}-21 & {\;}-15 & {\;}-108 & 
\phantom{-}216 \\
     \mbox{\footnotesize$3$}\; & \phantom{-}15 & \phantom{-}75 & 
\phantom{-}294 & {\;}-576 \\
     \mbox{\footnotesize$2$}\; & -1 & -93 & -198 & \phantom{-}396 & 
{\;}\phantom{-}72 \\
     \mbox{\footnotesize$n = 1$}\; & -3 & \phantom{-}53 & \phantom{-}18 
& -132 \\
     \mbox{\footnotesize$n = 0$}\; & \phantom{-}1 & -13 & \phantom{-}20 
\\
     \hline
     \zeta_{\mkern1mu m \mkern0.5mu n}\vphantom{\vert^\vert}\; & 
\mbox{\footnotesize$m = 0$} & \mbox{\footnotesize$m = 1$} & 
\mbox{\footnotesize$\phantom{-}2$}~ & 
\mbox{\footnotesize$\phantom{-}3$}~~ & 
\mbox{\footnotesize$\phantom{-}4$}~
   \end{array} \end{gather*}

 \vspace{-3pt}

 \begin{figure}
   \vspace{10pt}
   \centerline{\includegraphics{plane.pdf}}
   \vspace{3pt}
   \caption{The curve $\zeta(c_2 ,\mkern2mu} \def\DWS{\mkern3.5mu c_3) = 0$ (solid curve); the line 
$c_3 = \frac{1}{6}$ (dashed curve); and the hyperbola $6 c_3 - 12 c_3 
(c_3 - c_2) - 1 = 0$ (dotted curve, the coefficient $a_{43}$ is infinite 
on it) which intersects with the curve $\zeta = 0$ at the points 
$(\frac{1}{6} ,\mkern2mu} \def\DWS{\mkern3.5mu \frac{1}{6})$, $(\frac{1}{6} ,\mkern2mu} \def\DWS{\mkern3.5mu \frac{1}{2})$, 
$\mathsf{Q} = \bigl( (6 - \sqrt{21}) / 15 ,\mkern2mu} \def\DWS{\mkern3.5mu (9 - \sqrt{21}) / 12 
\bigr) = (0.09449..., 0.36811...)$, and $\mathsf{Q}' = \bigl( (6 + 
\sqrt{21}) / 15 ,\mkern2mu} \def\DWS{\mkern3.5mu (9 + \sqrt{21}) / 12 \bigr) = (0.70550..., 
1.13188...)$. Some branches of the curve $\zeta = 0$ are not shown as 
they correspond to too large values of $c_2$ or $c_3$ to be visible on 
the graph. The points $\mathsf{C}$ and $\mathsf{A}$ are indicated by 
open dots and correspond to the method in eq.~(\ref{method_C}) and the 
$1$-dimensional family of methods that includes the method in 
eq.~(\ref{the_method}), respectively.}
 \end{figure}

\noindent Both coefficients $\phi_{20}$ and $\phi_{30}$ are equal to zero at 
the point $(c_2 ,\mkern2mu} \def\DWS{\mkern3.5mu c_3) = \mathsf{C} = (2 \gamma ,\mkern2mu} \def\DWS{\mkern3.5mu 4 \gamma)$, 
where $\gamma = \bigl( 2 + 2^{1 / 3} + 2^{-1 / 3} \bigr) / 12 = 1 / 4 
\mkern0.5mu (2 - 2^{1/3}) = 0.33780...$ There $\phi_{2b} = 0$ and 
$\phi_{3b} = 2 \gamma - \frac{1}{2}$, thus $a_{53} = \frac{1}{2} (-5 + 
48 \gamma - 96 \gamma^2) \mkern1mu a_{52} + \frac{1}{2} - 2 \gamma$, 
also $b_5 a_{52} = 0$ (from $D_2{(\VEC{1} \mkern-0.75mu )}$) and $a_{65} a_{52} = 0$ (from 
$(\MAT{A} \mkern0.5mu \VEC{q}_1 \mkern-0.5mu )_3 + (\MAT{A} \mkern0.5mu \VEC{q}_1 \mkern-0.5mu )_6 = 0$). In both cases of $a_{52} = 0$ (then 
$\VEC{a}_{4*} = \VEC{a}_{5*}$, \IE, the \NTH{4} and the \NTH{5} stages 
duplicate each other) and $a_{52} \ne 0$ (then $b_5 = 0$ and 
$\VEC{a}_{*5} = \VEC{0}$, \IE, the \NTH{5} stage is not used), the 
method becomes equivalent to

 \vspace{-12pt}

 \begin{gather}
   \begin{array}{r|cccccccc}
     0\;\,\vphantom{\vert_\vert} \\
     {2 \gamma}\;\,\vphantom{\vert_\vert^\vert} & \phantom{-}2 \gamma \\
     {4 \gamma}\;\,\vphantom{\vert_\vert^\vert} & \phantom{-}0 & \;\,\phantom{-}4 \gamma \\
     \frac{1}{2}\;\,\vphantom{\vert_\vert^\vert} & \phantom{-}2 \gamma & \;\;\phantom{-}0 & 
\phantom{-}\frac{1}{2} - 2 \gamma \\
     {1 - 4 \gamma}\;\,\vphantom{\vert_\vert^\vert} & \phantom{-}0 & \;\;\phantom{-}4 \gamma & 
\phantom{-}0 & \phantom{-}1 - 8 \gamma \\
     {1 - 2 \gamma}\;\,\vphantom{\vert_\vert^\vert} & \phantom{-}2 \gamma & \;\;\phantom{-}0 & 
\phantom{-}\frac{1}{2} - 2 \gamma & \phantom{-}0 & 
\phantom{-}\frac{1}{2} - 2 \gamma \\
     {1}\;\,\vphantom{\vert_\vert^\vert} & \phantom{-}0 & \;\;\phantom{-}4 \gamma & 
\phantom{-}0 & \phantom{-}1 - 8 \gamma & \phantom{-}0 & \,\phantom{-}4 
\gamma \\
     \hline
     \vphantom{\vert^{\vert^{\vert^\vert}}} & \phantom{-}\gamma & 
\;\;\phantom{-}2 \gamma & \phantom{-}\frac{1}{4} - \gamma & 
\phantom{-}\frac{1}{2} - 4 \gamma & \phantom{-}\frac{1}{4} - \gamma & 
\;\phantom{-}2 \gamma & \;\;\,\phantom{-}\gamma
   \end{array} \label{method_C} \end{gather}

 \vspace{-3pt}

\noindent This $7$-stage pseudo-symplectic Runge--Kutta method of order $(4 
,\mkern2mu} \def\DWS{\mkern3.5mu 9)$ is strongly connected to a well-known explicit symplectic 
method for separated Hamiltonians \citep[eqs.~(4.9) and (4.12)]{FoRu90}, 
\citep[eq.~(2.11)]{Yos90}.

Besides the point $\mathsf{C}$ one must have $b_5 \ne 0$, and thus 
$a_{52} = \phi_{20} / b_5 - \phi_{2b}$ and $a_{53} = \phi_{30} / b_5 - 
\phi_{3b}$. The cubic equation $z \mkern1mu (z - \frac{1}{2}) (z - 1) = 
\frac{1}{24}$ has three roots: $z_1 = 1 / 2 - \sin(2 \pi / 9) / \sqrt{3} 
= 0.12888...$, $z_2 = 1 / 2 - \sin(2 \pi / 9 + 2 \pi / 3) / \sqrt{3} = 1 
/ 8 \sin^2 (2 \pi / 9) = 0.30253...$, and $z_3 = 1 / 2 - \sin(2 \pi / 9 
+ 4 \pi / 3) / \sqrt{3} = 1.06857...$ The quantity $b_5 \bigl( (\MAT{A} \mkern0.5mu \VEC{q}_1 \mkern-0.5mu )_3 
+ (\MAT{A} \mkern0.5mu \VEC{q}_1 \mkern-0.5mu )_6 \bigr)$, which should be equal to zero, is a linear function 
of $b_5$ and $a_{65}$, and its derivative with respect to $b_5$ is not 
equal to zero besides the points $(c_2 ,\mkern2mu} \def\DWS{\mkern3.5mu c_3) = \mathsf{A} = (z_1 
,\mkern2mu} \def\DWS{\mkern3.5mu z_2)$, $\mathsf{A}' = (z_2 ,\mkern2mu} \def\DWS{\mkern3.5mu z_3)$, and $\mathsf{A}'' = (z_3 
,\mkern2mu} \def\DWS{\mkern3.5mu z_1)$. There the weight $b_5$ is determined from $(\MAT{A} \mkern0.5mu \VEC{q}_1 \mkern-0.5mu )_3 + 
(\MAT{A} \mkern0.5mu \VEC{q}_1 \mkern-0.5mu )_6 = 0$, and then $a_{65}$ is found from $(\MAT{A} \mkern0.5mu \VEC{q}_1 \mkern-0.5mu )_8 = 0$. The 
vector $\MAT{A} \mkern0.5mu \VEC{q}_1 = \bigl[ \; 0 ~~~ 0 ~~~ 
{-{\mkern-1mu}\frac{c_2 (6 c_3 - 1)}{48 (1 - c_2)}} ~~~ 0 ~~~ 0 ~~~ 
{\frac{c_2 (6 c_3 - 1)}{48 (1 - c_2)}} ~~ 0 ~~~ 0 \; 
\bigr]{}^{\textrm{T}}$ is now ``odd''. The result is a $1$-dimensional 
(along the curve $\zeta(c_2 ,\mkern2mu} \def\DWS{\mkern3.5mu c_3) = 0$) family of pseudo-symplectic 
methods of order $(4 ,\mkern2mu} \def\DWS{\mkern3.5mu 6)$ that satisfy $\VEC{D}{(\VEC{1} \mkern-0.75mu )}$, 
$\VEC{D}(\VEC{c})$, and $\VEC{D}(\VEC{c}^2)$, but not $\VEC{D}(\MAT{A} 
\mkern0.5mu \VEC{c})$. For example, at the point $(c_2 ,\mkern2mu} \def\DWS{\mkern3.5mu c_3) = 
\mathsf{R} = (\frac{1}{6} ,\mkern2mu} \def\DWS{\mkern3.5mu \frac{11}{30})$ a method is equivalent 
(the \NTH{6} stage is not used) to

 \vspace{-12.5pt}

 \begin{gather*}
   \begin{array}{r|cccccccc}
     0\;\,\vphantom{\vert_\vert} \\
     {\frac{1}{6}}\;\,\vphantom{\vert_\vert^\vert} & \;\phantom{-}\frac{1}{6} \\
     {\frac{11}{30}}\;\,\vphantom{\vert_\vert^\vert} & \;\phantom{-}\frac{1}{150} & 
\;\phantom{-}\frac{9}{25} \\
     \frac{1}{2}\;\,\vphantom{\vert_\vert^\vert} & \;\phantom{-}\frac{1}{4} & \;-\frac{13}{48} 
& \;\phantom{-}\frac{25}{48} \\
     \frac{1}{2}\;\,\vphantom{\vert_\vert^\vert} & \;-\frac{37}{50} & 
\;\phantom{-}\frac{59}{25} & \;-2 & \;\phantom{-}\frac{22}{25} \\
     {\frac{5}{6}}\;\,\vphantom{\vert_\vert^\vert} & \;\phantom{-}\frac{1}{6} & \;\phantom{-}0 
& \;\phantom{-}0 & \;\phantom{-}\frac{4}{11} & 
\;\phantom{-}\frac{10}{33} \\
     1\;\,\vphantom{\vert_\vert^\vert}\vphantom{\vert_{\vert_\vert}} & \;\phantom{-}0 & 
\;\phantom{-}\frac{3}{8} & \;\phantom{-}0 & \;\phantom{-}\frac{4}{11} & 
\;-\frac{5}{44} & \;\phantom{-}\frac{3}{8} \\
     \hline
     \vphantom{\vert^{\vert^{\vert^\vert}}} & \;\phantom{-}\frac{1}{12} 
& \;\phantom{-}\frac{3}{16} & \;\phantom{-}0 & \;\phantom{-}\frac{4}{11} 
& \;\phantom{-}\frac{25}{264} & \;\phantom{-}\frac{3}{16} & 
\;\phantom{-}\frac{1}{12}
   \end{array} \end{gather*}

At all the three points $\mathsf{A}$, $\mathsf{A}'$, and $\mathsf{A}''$ 
the vector $\MAT{A} \mkern0.5mu \VEC{q}_1$ is ``odd'', plus 
$\VEC{D}{(\VEC{1} \mkern-0.75mu )}$ and $\VEC{D}(\VEC{c})$ are satisfied. The coefficient 
$a_{65}$ is determined from $D_5(\VEC{c}^2)$. Now $\VEC{D}(\VEC{c}^2)$ 
and $\VEC{D}(\MAT{A} \mkern0.5mu \VEC{c})$ are satisfied, and the result 
is a $1$-dimensional family of $8$-stage pseudo-symplectic Runge--Kutta 
methods of order $(4 ,\mkern2mu} \def\DWS{\mkern3.5mu 8)$ indexed by a parameter $\psi$:

 \vspace{-13.5pt}

 \begin{gather*}
   \begin{array}{r|cccccccc}
     0\;\,\vphantom{\vert_\vert^\vert} \\
     {c_2}\;\,\vphantom{\vert_\vert^\vert} & \;\,c_2 \\
     {c_3}\;\,\vphantom{\vert_\vert^\vert} & \;\,0 & \;\,c_3 \\
     \frac{1}{2}\;\,\vphantom{\vert_\vert^\vert} & \;\,\frac12 - c_2 & \;\,c_2 + c_3 - 1 & 
\;\,1 - c_3 \\
     \frac{1}{2}\;\,\vphantom{\vert_\vert^\vert} & \;\,a_{51} & \;\,a_{52} & \;\,a_{53} & 
\phantom{-}a_{54} \\
     {1 - c_3}\;\,\vphantom{\vert_\vert^\vert} & \;\,0 & \;\,c_3 & \;\,0 & \phantom{-}a_{64} & 
\phantom{-}a_{65} \\
     {1 - c_2}\;\,\vphantom{\vert_\vert^\vert} & \;\,c_2 & \;\,0 & \;\,\frac12 - 2 c_2 & 
\phantom{-}a_{74} & \phantom{-}a_{75} & \;\;\frac12 - 2 c_2 \\
     1\;\,\vphantom{\vert_\vert^\vert} & \;\,0 & \;\,c_3 & \;\,0 & \phantom{-}a_{64} & 
\phantom{-}a_{65} & \;\;0 & \phantom{-}c_3 \\
     \hline
     \vphantom{\vert^{\vert^{\vert^\vert}}} & \;\,\frac12 c_2 & 
\;\,\frac12 c_3 & \;\,\frac14 - c_2 & \phantom{-}b_4 & \phantom{-}b_5 & 
\;\;\frac14 - c_2 & \phantom{-}\frac12 c_3 & \phantom{-}\frac12 c_2
   \end{array} \\[6pt]
   \left[ \begin{array}{c} \VEC{a}_{4*}\! \\ \VEC{a}_{5*}\! \end{array} 
\right] = \left[ \begin{array}{cccccccc}
     \vphantom{\vert_\vert^\vert}{}\PHI c_2 & \phantom{a}(1 - \PHI) c_3 & \phantom{a}\PHI 
(\frac12 - 2 c_2) & \phantom{a}\PHI c_2 + (1 - \PHI) (\frac12 - c_3) & 
\phantom{m}0 & \phantom{-}0 & \phantom{-}0 & \phantom{-}0\, \\
     \vphantom{\vert_\vert^\vert}{}\psi c_2 & \phantom{a}(1 - \psi) c_3 & \phantom{a}\psi 
(\frac12 - 2 c_2) & \phantom{a}\psi c_2 + (1 - \psi) (\frac12 - c_3) & 
\phantom{m}0 & \phantom{-}0 & \phantom{-}0 & \phantom{-}0\,
   \end{array} \right] \\[6pt] \left[ \begin{array}{cc}
     a_{64} &  \;\,a_{65} \\
     a_{74} & \;\,a_{75} \\
     b_4 & \;\,b_5
   \end{array} \right] = \left[ \begin{array}{cc}
    \vphantom{\vert_\vert^\vert}{}2 (\frac12 - c_3) (1 - \chi) & \phantom{-}2 (\frac12 - c_3) 
\chi \\
    \vphantom{\vert_\vert^\vert}{}2 c_2 (1 + \chi) & -2 c_2 \chi \\
    \vphantom{\vert_\vert^\vert}{}\frac12 (a_{64} + a_{74}) & \phantom{-}\frac12 (a_{65} + 
a_{75})
  \end{array} \right] , \qquad
  \chi = \frac {4 (1 - 3 c_2)} {(1 - 6 c_2) (\PHI - \psi)}
 \end{gather*}

 \vspace{-4.5pt}

\noindent Here $c_2 (c_2 - \frac{1}{2}) (c_2 - 1) = \frac{1}{24}$, $c_3 = c_2 
(1 - c_2) / (\frac12 - c_2) = 1 / 6 (1 - 2 c_2)^2$, and $\PHI = (\frac12 
- c_3) / (\frac12 - c_2 - c_3) = 1 / 2 c_2 - 1$. Out of these three 
points, the methods at the point $(c_2 ,\mkern2mu} \def\DWS{\mkern3.5mu c_3) = \mathsf{A} = (z_1 
,\mkern2mu} \def\DWS{\mkern3.5mu z_2) = (0.12888... ,\mkern2mu} \def\DWS{\mkern3.5mu 0.30253...)$ are the most efficient.

The residuals in the properties $D \bigl( \VEC{\Phi}(\textrm{t}_1), 
\VEC{\Phi}(\textrm{t}_2) \bigr)$ with $\vert \textrm{t}_1 \vert = 4$ and 
$\vert \textrm{t}_2 \vert = 5$ are inversely proportional to $b_5$. In 
the limit $b_5 \to \infty$, or $\psi \to \PHI$, one gets a method that 
uses information about the derivative of the
 r.h.s.{\ }function $\VEC{f}$ \citep[s.~224]{But16}, 
\citep[s.~IV.7]{HW96}. As the sum $b_4 + b_5$ does not depend on $\psi$, 
maximizing $b_5$ but keeping $b_4 \ge 0$ leads to $b_4 = 0$, $\psi = 
(\frac{1}{2} - c_3) / (\frac{1}{2} + c_2 - c_3) = 2 c_3$:
 \begin{align} \begin{array}{r|cccccccc}
     0{\;\,}\vphantom{\vert_\vert^\vert} & & & & & & \multicolumn{3}{l}{\mbox{\small$c_2 = 
\frac{1}{2} - \frac{1}{\sqrt{3}} \sin \bigl( \mkern-1mu \frac{2 \pi}{9} 
\bigr)$}} \\
     c_2{\;\,}\vphantom{\vert_\vert^\vert} & \phantom{l}c_2 & & & & & 
\multicolumn{3}{l}{\mbox{\small$c_3 = \frac{1}{2} - \frac{1}{\sqrt{3}} 
\sin \bigl( \mkern-0.5mu \frac{\pi}{9} \bigr) $}} \\
     c_3{\;\,}\vphantom{\vert_\vert^\vert} & \phantom{l}0 & c_3 \\
     \frac{1}{2}{\;\,}\vphantom{\vert_\vert^\vert} & \phantom{l}\frac12 - c_2 & c_2 + c_3 - 1 
& 1 - c_3 \\
     \frac{1}{2}{\;\,}\vphantom{\vert_\vert^\vert} & \phantom{l}2 c_2 c_3 & (1 - 2 c_3) c_3 & 
(1 - 4 c_2) c_3 & 4 c_2 c_3 \\
     {1 - c_3}{\;\,}\vphantom{\vert_\vert^\vert} & \phantom{l}0 & c_3 & 0 & 4 c_2 - 2 & 
\frac{1}{2 c_2} - 2 \\
     {1 - c_2}{\;\,}\vphantom{\vert_\vert^\vert} & \phantom{l}c_2 & 0 & \frac12 - 2 c_2 & 2 - 
4 c_2 & 6 c_2 - 2 & \frac12 - 2 c_2 \\
     1{\;\,}\vphantom{\vert_\vert^\vert} & \phantom{l}0 & c_3 & 0 & 4 c_2 - 2 & \frac{1}{2 
c_2} - 2 & 0 & \phantom{n}c_3 \\
     \hline
     \vphantom{\vert^{\vert^{\vert^\vert}}} & \phantom{l}\frac12 c_2 & 
\frac12 c_3 & \frac14 - c_2 & 0 & \frac12 + c_2 - c_3 & \frac14 - c_2 & 
\phantom{n}\frac12 c_3 & \phantom{m}\frac12 c_2
   \end{array} \hskip-12pt \label{the_method} \end{align} All the nodes, 
weights, and coefficients are elements of the algebraic extension 
$\pmb{\mathbf{Q}}(c_2) = \bigl\{ \alpha_1 + \alpha_2 c_2 + \alpha_3 c_3 
\mkern2mu \big\vert \mkern2mu \alpha_1 , \alpha_2 , \alpha_3 \in 
\pmb{\mathbf{Q}} \bigr\}$ of the field of rational numbers 
$\pmb{\mathbf{Q}}$. For example, $c_2^2 = -\frac{1}{12} + \frac{7}{6} 
c_2 - \frac{1}{6} c_3$, $c_2 c_3 = -\frac{1}{12} + \frac{1}{6} c_2 + 
\frac{1}{3} c_3$, $c_3^2 = -\frac{1}{3} + \frac{1}{6} c_2 + \frac{4}{3} 
c_3$, and $a_{65} = a_{85} = 1 / 2 c_2 - 2 = \PHI - 1 = 3 - 4 c_2 - 2 
c_3$.

The stability function of the method in eq.~(\ref{the_method}) is equal 
to

\vspace{-14.5pt}

\begin{gather*}
   R(z) = 1 + z + {\textstyle} \def\d{{\rm d}} \def\PHI{\varphi\frac{1}{2}} z^2 + {\textstyle} \def\d{{\rm d}} \def\PHI{\varphi\frac{1}{6}} z^3 + 
{\textstyle} \def\d{{\rm d}} \def\PHI{\varphi\frac{1}{24}} z^4 + {\textstyle} \def\d{{\rm d}} \def\PHI{\varphi\frac{1}{120}} r_5 \mkern1mu z^5 + 
{\textstyle} \def\d{{\rm d}} \def\PHI{\varphi\frac{1}{720}} r_6 \mkern1mu z^6 + {\textstyle} \def\d{{\rm d}} \def\PHI{\varphi\frac{1}{5040}} r_7 \mkern1mu 
z^7 + {\textstyle} \def\d{{\rm d}} \def\PHI{\varphi\frac{1}{40320}} r_8 \mkern1mu z^8
 \end{gather*}

\vspace{-3.5pt}

\noindent where $r_5 = \frac{5}{6} + \frac{10}{3} c_2 - \frac{5}{6} c_3 = 
 1.0108...$, $r_6 = 20 c_2 - 5 c_3 = 6 r_5 - 5 = 1.0650...$, $r_7 = 
 -\frac{35}{8} + 105 c_2 - \frac{105}{4} c_3 = 1.2165...$, and $r_8 = 
 -\frac{70}{3} + \frac{1400}{3} c_2 - \frac{350}{3} c_3 =
1.5179...$

 \begin{table} \vspace{6pt} \begin{tabular}{r|rccclll}
     & $s$ & $p$ & $q$ & $10^4 \times T_4$ & $\phantom{-} 10^3 \times 
T_5$ & $\phantom{-} 10^3 \times T_6$ & $\phantom{-} R(z) R(-z) \,-\, 1$ 
\\
     \hline
     RK4$\vphantom{\vert_\vert^\vert}$ & $4$ & $4$ & $4$ & $0$ & $\phantom{0}14.504...$ & 
$\phantom{0}16.035...$ & $\phantom{-}0.01388888... \; z^6$ \\
     AC36 & $5$ & $3$ & $6$ & $7.5690...$ & $\phantom{00}2.3451...$ & 
$\phantom{00}3.9611...$ & $-0.00057870... \; z^8$ \\
     CLMR47$\vphantom{\vert_\vert^\vert}$ & $6$ & $4$ & $7$ & $0$ & $\phantom{00}0.29185...$ & 
$\phantom{00}0.39787...$ & $\phantom{-}0.00024760... \; z^8$ \\
     CCRL47 & $6$ & $4$ & $7$ & $0$ & $\phantom{00}1.4813...$ & 
$\phantom{00}1.6278...$ & $-0.00030365... \; z^8$ \\
     eq.~(\ref{method_C})$\vphantom{\vert_\vert^\vert}$ & $7$ & $4$ & $9$ & $0$ & 
$112.99...$ & $132.54...$ & $-0.00144678... \; z^{10}$ \\
     eq.~(\ref{the_method}) & $8$ & $4$ & $8$ & $0$ & 
$\phantom{00}0.64048...$ & $\phantom{00}0.91796...$ & $\phantom{-} 
0.00000950... \; z^{10}$ \\
     CV8$\vphantom{\vert_\vert^\vert}$ & $11\mkern-0.5mu$ & $8$ & $8$ & $0$ & 
$\phantom{0000}0$ & $\phantom{0000}0$ & $\phantom{-}0.00000627... \; 
z^{10}$ \\
     GL4$\vphantom{\vert_{\vert_\vert}}$ & $2$ & $4$ & $\infty$ & $0$ & 
$\phantom{00}4.3306...$ & $\phantom{00}5.6178...$ & $\phantom{-}0$
 \end{tabular} \begin{tabular}{r|cccccp{4pt}lp{0pt}l}
     & $C(2)\vphantom{\vert^{\big\vert}}$ & 
$\mkern-1mu{}\VEC{D}(\VEC{1})\mkern-1mu$ & $\VEC{D}(\VEC{c})$ & 
$\mkern-3mu{}\VEC{D}(\VEC{c}^2)\mkern-3mu$ & 
$\mkern-7mu{}\VEC{D}(\MAT{A} \mkern0.5mu \VEC{c})\mkern-7mu$ & & 
{$\mbox{max}_{ij} \vert \mkern1mu {a_{\mkern0.85mu i \mkern-0.85mu j}} \vert$} & &
{$\phantom{-}\mbox{min}_j \mkern1.5mu {b_{\mkern-1.25mu j}} \vphantom{\vert_\vert^\vert}$} \\
     \hline
     RK4$\vphantom{\vert_\vert^\vert}$ & \FALSE & $\textsf{T}$} \def\FALSE{$\tiny\textsf{F}$ & \FALSE & \FALSE & \FALSE & & 
$\phantom{0}1$ & & $\phantom{-}0.1666...$ \\
     AC36 & \FALSE & $\textsf{T}$} \def\FALSE{$\tiny\textsf{F}$ & $\textsf{T}$} \def\FALSE{$\tiny\textsf{F}$ & \FALSE & \FALSE & & 
$\phantom{0}2.1621...$ & & $-0.3054...$ \\
     CLMR47$\vphantom{\vert_\vert^\vert}$ & \FALSE & $\textsf{T}$} \def\FALSE{$\tiny\textsf{F}$ & $\textsf{T}$} \def\FALSE{$\tiny\textsf{F}$ & $\textsf{T}$} \def\FALSE{$\tiny\textsf{F}$ & $\textsf{T}$} \def\FALSE{$\tiny\textsf{F}$ & & 
$\phantom{0}4.4309...$ & & $\phantom{-}0.0277...$ \\
     CCRL47 & \FALSE & $\textsf{T}$} \def\FALSE{$\tiny\textsf{F}$ & $\textsf{T}$} \def\FALSE{$\tiny\textsf{F}$ & $\textsf{T}$} \def\FALSE{$\tiny\textsf{F}$ & $\textsf{T}$} \def\FALSE{$\tiny\textsf{F}$ & & 
$\phantom{0}2.9265...$ & & $-0.4488...$ \\
     eq.~(\ref{method_C})$\vphantom{\vert_\vert^\vert}$ & \FALSE & $\textsf{T}$} \def\FALSE{$\tiny\textsf{F}$ & $\textsf{T}$} \def\FALSE{$\tiny\textsf{F}$ & 
$\textsf{T}$} \def\FALSE{$\tiny\textsf{F}$ & $\textsf{T}$} \def\FALSE{$\tiny\textsf{F}$ & & $\phantom{0}1.7024...$ & & $-0.8512...$ \\
     eq.~(\ref{the_method}) & \FALSE & $\textsf{T}$} \def\FALSE{$\tiny\textsf{F}$ & $\textsf{T}$} \def\FALSE{$\tiny\textsf{F}$ & $\textsf{T}$} \def\FALSE{$\tiny\textsf{F}$ & $\textsf{T}$} \def\FALSE{$\tiny\textsf{F}$ & & 
$\phantom{0}1.8793...$ & & $\phantom{-}0.0644...$ \\
     CV8$\vphantom{\vert_\vert^\vert}$ & $\textsf{T}$} \def\FALSE{$\tiny\textsf{F}$ & $\textsf{T}$} \def\FALSE{$\tiny\textsf{F}$ & \FALSE & \FALSE & \FALSE & & 
$14.728...$ & & $\phantom{-}0.05$ \\
     GL4$\vphantom{\vert_{\vert_\vert}}$ & $\textsf{T}$} \def\FALSE{$\tiny\textsf{F}$ & $\textsf{T}$} \def\FALSE{$\tiny\textsf{F}$ & $\textsf{T}$} \def\FALSE{$\tiny\textsf{F}$ & $\textsf{T}$} \def\FALSE{$\tiny\textsf{F}$ 
& $\textsf{T}$} \def\FALSE{$\tiny\textsf{F}$ & & $\phantom{0}0.5386... $ & & $\phantom{-}0.5$ \end{tabular} 
\vspace{4pt} \caption{A comparison of eight $s$-stage Runge--Kutta 
methods of order $(p ,\mkern2mu} \def\DWS{\mkern3.5mu q)$. Error coefficients are defined as $T_p^2 
= \sum_{\;\textrm{t},\;\vert \mkern1mu \textrm{t} \mkern1mu \vert = p} 
\bigl( \VEC{b} \mkern2mu \VEC{\Phi}(\textrm{t}) - 1 / \textrm{t}! 
\bigr){}^2 \mkern-1mu / \mkern1mu \sigma^2(\textrm{t}) = (1 / p!)^2 
\sum_{\;\textrm{t},\;\vert \mkern1mu \textrm{t} \mkern1mu \vert = p} 
\alpha^2(t) \bigl( \textrm{t}! \mkern3mu \VEC{b} \mkern2mu 
\VEC{\Phi}(\textrm{t}) - 1 \bigr){}^2$, where $\sigma(\textrm{t})$ is 
the order of the symmetry group of the tree $\textrm{t}$, and 
$\alpha(\textrm{t})$ is the number of monotonic labelings of 
$\textrm{t}$ (see, \EG, \cite[ss.~304 and 318]{But16}, \citep[pp.~147 
and 158]{HNW93}, \citep[pp.~58 and 60]{But21}, \citep[pp.~57 and 
58]{HLW06}). The column about the function $R(z) R(-z) - 1$ shows the 
first non-zero term in its Taylor expansion about $z = 0$. Within the 
$1$-dimensional families of methods at the points $\mathsf{A}$, 
$\mathsf{A}'$, and $\mathsf{A}''$ the error coefficients $T_5$, $T_6$, 
and the stability function $R(z)$ do not depend on the parameter $\psi$. 
The symbols $\textsf{T}$} \def\FALSE{$\tiny\textsf{F}$$\mkern2mu$/$\mkern3mu$\raisebox{0.5pt}{\FALSE} stand 
for true$\mkern1mu$/$\mkern1mu$false. The explicit version of the 
property $C(2)$ is used. The $\mbox{min}_j \mkern1.5mu {b_{\mkern-1.25mu j}}$ column shows 
the minimal value of a non-zero weight.} \label{table_methods} 
\end{table}

 \begin{figure} \centerline{\includegraphics{stair.pdf} \vspace{6pt}} 
\caption{Quadrature scheme graphical depiction for the eight methods 
listed in Table~\ref{table_methods}. Here ${\Sigma_{\mkern1mu i} 
\mkern1.5mu \VEC{b}} = \sum_{j = 1}^{\mkern1.5mu i} {b_{\mkern-1.25mu j}}$, where $1 \le i 
\le s$, are the cumulative weights, with $\Sigma_{\mkern1mu 0} 
\mkern1.5mu \VEC{b} = 0$. Due to the order condition $\VEC{b} 
\mkern0.5mu \VEC{1} = 1$ for an $s$-stage method one has 
$\Sigma_{\mkern1mu s} \mkern1.75mu \VEC{b} = 1$. On the diagram the 
points $\bigl( {c_{\mkern0.5mu i}} ,\mkern2mu} \def\DWS{\mkern3.5mu {\Sigma_{\mkern1mu i - 1} \mkern1.5mu \VEC{b}} 
\bigr)$ and $\bigl( {c_{\mkern0.5mu i}} ,\mkern2mu} \def\DWS{\mkern3.5mu {\Sigma_{\mkern1mu i} \mkern1.5mu \VEC{b}} 
\bigr)$, where $1 \le i \le s$, are connected by a thick line (with 
white filling if ${\Sigma_{\mkern1mu i} \mkern1.5mu \VEC{b}} < 
{\Sigma_{\mkern1mu i - 1} \mkern1.5mu \VEC{b}}$). To easier follow the 
progression of the stages, the points $\bigl( {c_{\mkern0.5mu i}} ,\mkern2mu} \def\DWS{\mkern3.5mu 
{\Sigma_{\mkern1mu i} \mkern1.5mu \VEC{b}} \bigr)$ and $\bigl( 
c_{\mkern0.5mu i + 1} ,\mkern2mu} \def\DWS{\mkern3.5mu {\Sigma_{\mkern1mu i} \mkern1.5mu \VEC{b}} 
\bigr)$, where $1 \le i < s$, are connected by dotted lines.} 
\label{stair} \end{figure}

The basic properties of new methods eq.~(\ref{method_C}), 
eq.~(\ref{the_method}) and of some known Runge--Kutta methods are 
compared in Figure~\ref{stair} and in Table~\ref{table_methods}, where 
(and also in Figures~\ref{rigid} and \ref{pendulum}) the latter methods 
are named as follows: RK4 is the classical Runge--Kutta method 
\citep[p.~448]{Kut01}; AC36 is 
\citep[fig.~4.1]{AuCh98a};$\mkern1mu$\footnote{~In 
\citep[p.~453]{AuCh98a} it is stated that a certain optimization over 
methods satisfying $\hat{D}([\tau])$ (same as property 
$\VEC{D}(\VEC{c})$) results in $c_2 = 0.000050587...$ Nevertheless, 
optimization over methods not necessarily satisfying $\hat{D}([\tau])$ 
led to the method \citep[fig.~4.1]{AuCh98a} that satisfies 
$\hat{D}([\tau])$ and has $c_2 = 0.13502...$} CLMR47 is 
\citep[p.~262]{CLMR10};$\mkern1mu$\footnote{~In \citep[p.~262]{CLMR10} 
the denominators ``355,568'', ``831,328'', and ``9,955,904'' in 
$a_{42}$, $a_{51}$, and $a_{62}$ should be read as ``3,555,680'', 
``8,313,280'', and ``99,559,040'', respectively.} CCRL47 is 
\citep[p.~90]{CCRL17}; CV8 is \citep[tab.~1]{CoVe72}, 
\citep[p.~210]{But16}; and GL4 is the $2$-stage Gauss--Legendre method 
\citep[]{HaHo55}, \citep[p.~56]{But64a}.

\section{Numerical tests} \label{numerical_tests}

\hskip\parindent{}In this section the Runge--Kutta methods from 
Table~\ref{table_methods} are tested on two numerical examples: 
torque-free rotation of a rigid body \citep[s.~36]{LaLi76} and a 
pendulum with a non-separable Hamiltonian \citep[ex.~16.2, 
p.~313]{HNW93}.

Euler's equations describing free rotation of a rigid body with 
principal moments of inertia $I_1 = 1$, $I_2 = 2$, and $I_3 = 3$ are $\d 
\omega_{\mkern1mu 1} / \d \mkern1mu t = -\omega_{\mkern1.5mu 2} \omega_{\mkern1.5mu 3}$, $\d \omega_{\mkern1.5mu 2} / \d \mkern1mu t = \omega_{\mkern1mu 1} \omega_{\mkern1.5mu 3}$, and 
$\d \omega_{\mkern1.5mu 3} / \d \mkern1mu t = -\frac{1}{3} \omega_{\mkern1mu 1} \omega_{\mkern1.5mu 2}$. The solution with 
initial condition $(\omega_{\mkern1mu 1} ,\mkern2mu} \def\DWS{\mkern3.5mu \omega_{\mkern1.5mu 2} ,\mkern2mu} \def\DWS{\mkern3.5mu \omega_{\mkern1.5mu 3}) \big\vert_{\mkern1mu 
t = 0} = (12 ,\mkern2mu} \def\DWS{\mkern3.5mu 0 ,\mkern2mu} \def\DWS{\mkern3.5mu 7)$ is shown in Figure~\ref{csdn}. The pair of 
quadratic invariants, twice the kinetic energy $2 T = I_1 
\omega_{\mkern1mu 1}^{\mkern1mu 2} + I_2 \omega_{\mkern1.5mu 2}^{\mkern1mu 2} + I_3 \omega_{\mkern1.5mu 3}^{\mkern1mu 
2} = 291$ and squared magnitude of the angular momentum $L^2 = I_1^2 
\omega_{\mkern1mu 1}^{\mkern1mu 2} + I_2^2 \omega_{\mkern1.5mu 2}^{\mkern1mu 2} + I_3^2 
\omega_{\mkern1.5mu 3}^{\mkern1mu 2} = 585$, can be reduced to $Q_1(\VEC{\omega}) = 
\omega_{\mkern1mu 1}^{\mkern1mu 2} + \omega_{\mkern1.5mu 2}^{\mkern1mu 2} = 144$ and 
$Q_2(\VEC{\omega}) = \omega_{\mkern1.5mu 2}^{\mkern1mu 2} + 3 \omega_{\mkern1.5mu 3}^{\mkern1mu 2} = 
147$. When this system of differential equations is solved numerically, 
the values of the invariants $Q_1$ and $Q_2$ gradually drift from their 
initial values $144$ and $147$, see Figure~\ref{rigid}. Due to round-off 
errors, $Q_1 \bigl( \VEC{\tilde\omega}(n h) \bigr)$ and $Q_2 \bigl( 
\VEC{\tilde\omega}(n h) \bigr)$ as functions of $n$ can be approximately 
viewed as random walks. In the absence of a systematic drift the typical 
deviation of $(Q_1 ,\mkern2mu} \def\DWS{\mkern3.5mu Q_2)$ from $(144 ,\mkern2mu} \def\DWS{\mkern3.5mu 147)$ grows roughly 
proportional to $\sqrt{n}$, \EG, the slope of the curves near the bottom 
of Figure~\ref{rigid}, where the round-off errors dominate, is about 
$\frac{1}{2}$. Even a small bias makes the deviation to eventually, \IE, 
for large enough $n$, be proportional to $n$: in the log-log plots of 
Figure~\ref{rigid} the slopes of most of the curves are equal to $1$. 
The performances of AC36, CLMR47, and CCRL47 methods are close, with 
AC36 doing the best out of the three. Besides the implicit GL4, the 
methods in eq.~(\ref{method_C}) and eq.~(\ref{the_method}), and also 
CV8, due to their higher pseudo-symplecticity order have the fastest 
rate of decrease of the error in $Q_1$ and $Q_2$ with the step size $h$ 
(see the two right panels in Figure~\ref{rigid}). The method in 
eq.~(\ref{the_method}), out of the explicit methods tested, has the best 
performance.

 \begin{figure} \centerline{\includegraphics{csdn.pdf} \vspace{6pt}} 
\caption{A solution of Euler's equations with moments of inertia $I_1 = 
1$, $I_2 = 2$, and $I_3 = 3$: $(\omega_{\mkern1mu 1} ,\mkern2mu} \def\DWS{\mkern3.5mu \omega_{\mkern1.5mu 2} ,\mkern2mu} \def\DWS{\mkern3.5mu \omega_{\mkern1.5mu 3}) ( \mkern0.25mu t \mkern-0.5mu ) = 
(12 \mkern1mu \textrm{cn} ,\mkern2mu} \def\DWS{\mkern3.5mu 12 \mkern1mu \textrm{sn} ,\mkern2mu} \def\DWS{\mkern3.5mu 7 \mkern1mu 
\textrm{dn}) \bigl( 7 \mkern0.5mu t ,\mkern2mu} \def\DWS{\mkern3.5mu \frac{48}{49} \bigr)$, where 
$\textrm{cn}$, $\textrm{sn}$, and $\textrm{dn}$ are the Jacobi elliptic 
functions. The motion is periodic with the period $4 K \bigl( 
\frac{48}{49} \bigr) / 7 = 1.9109...$, where $K(m) = 
\smash{\int_{\mkern1mu 0}^{\mkern1mu \pi/2}} \d \theta \mkern1.5mu / 
\bigl( 1 - m \mkern1mu \sin^2 \theta \bigr){}^{1/2}$ is the complete 
elliptic integral of the \FIRST{} kind. The circular and elliptic 
cylinders correspond to the quadratic invariants $Q_1(\VEC{\omega}) = 
\omega_{\mkern1mu 1}^{\mkern1mu 2} + \omega_{\mkern1.5mu 2}^{\mkern1mu 2} = 144$ and 
$Q_2(\VEC{\omega}) = \omega_{\mkern1.5mu 2}^{\mkern1mu 2} + 3 \omega_{\mkern1.5mu 3}^{\mkern1mu 2} = 
147$, respectively. The thin lines on them divide the period into $72$ 
equal parts. Points of the trajectory that are connected by vertical 
lines with $\omega_{\mkern1mu 1}$ and $\omega_{\mkern1.5mu 2}$ axes are $(12 ,\mkern2mu} \def\DWS{\mkern3.5mu 0 ,\mkern2mu} \def\DWS{\mkern3.5mu 7)$ and $(0 
,\mkern2mu} \def\DWS{\mkern3.5mu 12 ,\mkern2mu} \def\DWS{\mkern3.5mu 1)$. The alternating behavior of $\omega_{\mkern1.5mu 2}( \mkern0.25mu t \mkern-0.5mu )$ is due to 
the rotation around the intermediate principal axis being unstable 
\cite[p.~48]{Poi34}, \citep[s.~37]{LaLi76}.} \label{csdn} \end{figure}

 \begin{figure}
   \vspace{7pt}
   \centerline{\includegraphics{rigid.pdf}}
   \vspace{3pt}
 \caption{Growth of the error in $Q_1(\VEC{\tilde\omega})$ and 
$Q_2(\VEC{\tilde\omega})$ with time $t$. Panels from left to right: (1) 
RK4 (curves from top to near bottom) and GL4 (bottom curves); (2) AC36 
(lower curves in groups of three), CLMR47 (middle curves), and CCRL47 
(upper curves); (3) CV8 and eq.~(\ref{method_C}) (lower and upper curves 
in pairs); and (4) eq.~(\ref{the_method}). The step size is $h = s 
\mkern1mu h_1$, where $s$ is the number of stages. Upper solid, dotted, 
and dashed curves correspond to $h_1 = 2^{-7}$, $h_1 = 2^{-8}$, and $h_1 
= 2^{-9}$, respectively. Next to upper solid or dotted curves correspond 
to $h_1 = 2^{-10}$ or $h_1 = 2^{-11}$, respectively, and so on, \IE, 
going to the next, with a different line style, curve down divides $h_1$ 
by $2$. Upper solid curves, upper dotted curves, and so on, on different 
panels correspond to the same value of $h_1$, and thus to the same 
number of the r.h.s. function $\VEC{f}(t ,\mkern2mu} \def\DWS{\mkern3.5mu \VEC{x})$ evaluations.} 
\label{rigid} \end{figure}

 \begin{figure}
   \vspace{10pt}
   \centerline{\includegraphics{pendulum.pdf}}
   \vspace{3pt}
 \caption{The speed of the systematic drift of the Hamiltonian 
$\mathcal{H}$ as the function of the step size $h = s \mkern1mu h_1$, 
where $s$ is the number of stages. The number of the r.h.s. function 
$\VEC{f}(t ,\mkern2mu} \def\DWS{\mkern3.5mu \VEC{x})$ evaluations is equal to the time duration of 
the simulation (which is $10020000$ here) divided by $h_1$, and with 
$h_1$ being fixed is independent of the method. The curves correspond to 
RK4, AC36, eq.~(\ref{the_method}), and GL4 (solid curves from top to 
bottom), CLMR47 and CV8 (upper and lower dotted curves), and CCRL47 and 
eq.~(\ref{method_C}) (upper and lower dashed curves). Three thin solid 
lines on this log-log plot have slopes $5$, $7$, and $9$, which 
corresponds to the speed of the drift being proportional to $h^5$, 
$h^7$, and $h^9$. The weighted moving average of the Hamiltonian that 
was used is $\langle \mathcal{H} \rangle ( \mkern0.25mu t \mkern-0.5mu ) = \sum_{\mkern1.5mu n 
\mkern1.5mu \in \mathcal{S} ( \mkern0.25mu t \mkern-0.5mu )} \sin^2 \bigl( \pi (n h - t) / 20000 
\bigr) \mkern1mu \mathcal{H} \bigl( \tilde{p}(n h) ,\mkern2mu} \def\DWS{\mkern3.5mu \tilde{x}(n h) 
\bigr) \mkern1mu \big/ \mkern1mu \sum_{\mkern1.5mu n \mkern1.5mu \in 
\mathcal{S} ( \mkern0.25mu t \mkern-0.5mu )} \sin^2 \bigl( \pi (n h - t) / 20000 \bigr)$, where 
$\mathcal{S} ( \mkern0.25mu t \mkern-0.5mu ) = \bigl\{ n \mkern2.5mu \big\vert \mkern3.5mu t < n h 
< t + 20000 \bigr\}$. The simulations were run at the University of 
Arizona High Performance Computing center.} \label{pendulum} 
\end{figure}

The \SECOND{} numerical example is a system $\d \mkern0.5mu x / \d \mkern1mu t = 
\partial \mathcal{H} / \partial p$, $\d p / \d \mkern1mu t = -\partial \mathcal{H} 
/ \partial x$, with the Hamiltonian function $\mathcal{H}(p ,\mkern2mu} \def\DWS{\mkern3.5mu x) = 
\frac{1}{2} p^2 - \bigl( 1 - \frac{1}{6} p \bigr) \cos(x)$. The solution 
with initial condition $(x ,\mkern2mu} \def\DWS{\mkern3.5mu p) \big\vert_{\mkern1mu t = 0} = \bigl( 
\arccos(-0.8) ,\mkern2mu} \def\DWS{\mkern3.5mu 0 \bigr)$ forms a periodic trajectory along the 
closed level curve $\mathcal{H}(p ,\mkern2mu} \def\DWS{\mkern3.5mu x) = 0.8$ around the origin. 
Because of notable variations of the Hamiltonian within a period in a 
numerical solution (see, \EG, Figure~\ref{linear}), here the speed of 
the systematic drift of the Hamiltonian $\mathcal{H}$ is estimated 
through the difference between weighted moving average values at two 
different moments of time. For the eight methods in 
Table~\ref{table_methods}, how the drift of the Hamiltonian depends on 
the step size $h$, is shown in Figure~\ref{pendulum}. Again, the 
performances of AC36, CLMR47, and CCRL47 methods are close, with AC36 
doing the best out of the three. In the triple of higher order explicit 
methods, the ones in eq.~(\ref{method_C}) and eq.~(\ref{the_method}) 
outperform CV8.\footnote{~Although with the same amount of work the 
drift of the Hamiltonian for CV8 is about $15$ times faster than for the 
method in eq.~(\ref{the_method}), the trajectory $\bigl( x(t) ,\mkern2mu} \def\DWS{\mkern3.5mu p(t) 
\bigr)$ is computed more accurately by the \NTH{8} order method CV8.}

 \section{Conclusions}

\hskip\parindent{}There are known $5$- and $6$-stage pseudo-symplectic 
Runge--Kutta methods of order $(3 ,\mkern2mu} \def\DWS{\mkern3.5mu 6)$ and $(4 ,\mkern2mu} \def\DWS{\mkern3.5mu 7)$, 
respectively, see, \EG, Table~\ref{table_methods}. With $7$ stages it is 
possible to come up with a method of order $(4 ,\mkern2mu} \def\DWS{\mkern3.5mu 9)$, see 
eq.~(\ref{method_C}). Utilising $8$ stages, one can construct robust, 
with non-negative weights and monotonically increasing nodes, methods of 
order $(4 ,\mkern2mu} \def\DWS{\mkern3.5mu 8)$, see eq.~(\ref{the_method}). The newly derived 
methods, largely due to their higher order, have better quadratic 
invariants and energy preservation properties than previously known 
pseudo-symplectic methods.

 \begin{figure}[H]
   \vspace{10pt}
   \centerline{\includegraphics{linear.pdf}}
   \vspace{3pt}
 \caption{The Hamiltonian $\mathcal{H} \bigl( \tilde{p}(t) ,\mkern2mu} \def\DWS{\mkern3.5mu 
\tilde{x}(t) \bigr)$ computed on the numerical solution obtained by the 
method in eq.~(\ref{the_method}). The vertical scale in the lower panel 
(time step $h = \frac{1}{16}$) is $16$ times larger than in the upper 
one ($h = \frac{1}{8}$). Similarity of heights of the depicted 
Hamiltonian deviations in the two panels is the consequence of the 
method being of order $4$. In order to show both the deviations small 
time scale structure and the systematic drift of the Hamiltonian, from 
the whole simulation (which lasts $1000500$) eleven fragments of 
duration $1000$ are shown. As the method conserves the symplectic 
structure up to the order $8$, the speed of the systematic drift 
decreases with $h$ much faster than the magnitude of the deviations 
within a period.} \label{linear} \end{figure}

 