Implicit A-Stable Peer Triplets for ODE Constrained Optimal Control Problems

Jens Lang
Technical University Darmstadt, Department of Mathematics
Dolivostraße 15, 64293 Darmstadt, Germany
lang@mathematik.tu-darmstadt.de

Bernhard A. Schmitt
Philipps-Universität Marburg, Department of Mathematics,
Hans-Meerwein-Straße 6, 35043 Marburg, Germany
schmitt@mathematik.uni-marburg.de

Abstract

This paper is concerned with the construction and convergence analysis of novel implicit Peer triplets of two-step nature with four stages for nonlinear ODE constrained optimal control problems. We combine the property of superconvergence of some standard Peer method for inner grid points with carefully designed starting and end methods to achieve order four for the state variables and order three for the adjoint variables in a first-discretize-then-optimize approach together with A-stability. The notion triplets emphasizes that these three different Peer methods have to satisfy additional matching conditions. Four such Peer triplets of practical interest are constructed. Also as a benchmark method, the well-known backward differentiation formula BDF4, which is only $A(73.35^{o})$ -stable, is extended to a special Peer triplet to supply an adjoint consistent method of higher order and BDF type with equidistant nodes. Within the class of Peer triplets, we found a diagonally implicit $A(84^{o})$ -stable method with nodes symmetric in $[0,1]$ to a common center that performs equally well. Numerical tests with three well established optimal control problems confirm the theoretical findings also concerning A-stability.

Key words. Implicit Peer two-step methods, BDF-methods, nonlinear optimal control, first-discretize-then-optimize, discrete adjoints

1 Introduction

The design of efficient time integrators for the numerical solution of optimal control problems constrained by systems of ordinary differential equations (ODEs) is still an active research field. Such systems typically arise from semi-discretized partial differential equations describing, e.g., the dynamics of heat and mass transfer or fluid flow in complex physical systems. Symplectic one-step Runge-Kutta methods [14, 21] exploit the Hamiltonian structure of the first-order optimality system - the necessary conditions to find an optimizer - and automatically yield a consistent approximation of the adjoint equations, which can be used to compute the gradient of the objective function. The first-order symplectic Euler, second-order Störmer-Verlet and higher-order Gauss methods are prominent representatives of this class, which are all implicit for general Hamiltonian systems, see the monograph [8]. Generalized partitioned Runge-Kutta methods which allow to compute exact gradients with respect to the initial condition are studied in [16]. To avoid the solution of large systems of nonlinear equations, stabilized explicit Runge-Kutta-Chebyshev methods have been recently proposed, too [2]. However, as all one-step methods, also symplectic Runge-Kutta schemes join the structural suffering of order reductions, which may lead to poor results in their application, e.g., to boundary control problems such as external cooling and heating in a manufacturing process; see [15, 18] for a detailed study of this behaviour.

In contrast, multistep methods including Peer two-step methods avoid order reductions and allow a simple implementation [5, 6]. However, the discrete adjoint schemes of linear multistep methods are in general not consistent or show a significant decrease of their approximation order [1, 20]. Recently, we have developed implicit Peer two-step methods [12] with three stages to solve ODE constrained optimal control problems of the form

$\displaystyle\mbox{minimize }C\big{(}y(T)\big{)}$		(1)
$\displaystyle\mbox{subject to }y^{\prime}(t)=$	$\displaystyle\,f\big{(}y(t),u(t)\big{)},\quad u(t)\in U_{ad},\;t\in(0,T],$	(2)
$\displaystyle y(0)=$	$\displaystyle\,y_{0},$	(3)

with the state $y(t)\in{\mathbb{R}}^{m}$ , the control $u(t)\in{\mathbb{R}}^{d}$ , $f:{\mathbb{R}}^{m}\times{\mathbb{R}}^{d}\mapsto{\mathbb{R}}^{m}$ , the objective function $C:{\mathbb{R}}^{m}\mapsto{\mathbb{R}}$ , where the set of admissible controls $U_{ad}\subset{\mathbb{R}}^{d}$ is closed and convex. Introducing for any $u\in U_{ad}$ the normal cone mapping

\displaystyle N_{U}(u)=

\displaystyle\,\{w\in{\mathbb{R}}^{d}:w^{T}(v-u)\leq 0\mbox{ for all }v\in U_{ad}\},

(4)

the first-order Karush-Kuhn-Tucker (KKT) optimality system [7, 25] reads

$\displaystyle y^{\prime}(t)=$	$\displaystyle\,f\big{(}y(t),u(t)\big{)},\quad t\in(0,T],\quad y(0)=y_{0},$	(5)
$\displaystyle p^{\prime}(t)=$	$\displaystyle\,-\nabla_{y}f\big{(}y(t),u(t)\big{)}^{\sf T}p(t),\quad t\in[0,T),\quad p(T)=\nabla_{y}C\big{(}y(T)\big{)}^{\sf T},$	(6)
	$\displaystyle\,-\nabla_{u}f\big{(}y(t),u(t)\big{)}^{\sf T}p(t)\in N_{U}\big{(}u(t)\big{)},\quad t\in[0,T].$	(7)

In this paper, we assume the existence of a unique local solution $(y^{\star},p^{\star},u^{\star})$ of the KKT system with sufficient regularity properties to justify the use of higher order Peer triplets, see, e.g., the smoothness assumption in [7, Section 2].

Following a first-discretize-and-then-optimize approach, we apply an $s$ -stage implicit Peer two-step method to (2)-(3) with approximations $Y_{ni}\approx y(t_{n}+c_{i}h)$ and $U_{ni}\approx u(t_{n}+c_{i}h)$ , $i=1,\ldots,s,$ on an equi-spaced time grid $\{t_{0},\ldots,t_{N+1}\}\subset[0,T]$ with step size $h=(T-t_{0})/(N+1)$ and nodes $c_{1},\ldots,c_{s}$ , which are fixed for all time steps, to get the discrete constraint nonlinear optimal control problem

$\displaystyle\mbox{minimize }C\big{(}y_{h}(T)\big{)}$		(8)
$\displaystyle\mbox{subject to }A_{0}Y_{0}=$	$\displaystyle\,a\otimes y_{0}+hb\otimes f(y_{0},u_{0})+hK_{0}F(Y_{0},U_{0}),$	(9)
$\displaystyle A_{n}Y_{n}=$	$\displaystyle\,B_{n}Y_{n-1}+hK_{n}F(Y_{n},U_{n}),\ n=1,\ldots,N,$	(10)
$\displaystyle y_{h}(T)=$	$\displaystyle\,(w^{T}\otimes I)Y_{N},$	(11)

with long vectors $Y_{n}=(Y_{ni})_{i=1}^{s}\in{\mathbb{R}}^{sm}$ , $U_{n}=(U_{ni})_{i=1}^{s}\in{\mathbb{R}}^{sd}$ , and $F(Y_{n},U_{n})=\big{(}f(Y_{ni},U_{ni})\big{)}_{i=1}^{s}$ . Further, $y_{h}(T)\approx y(T)$ , $u_{0}\approx u(0)$ , and $a,b,w\in{\mathbb{R}}^{s}$ are additional parameter vectors at both interval ends, $A_{n},B_{n},K_{n}\in{\mathbb{R}}^{s\times s}$ , and $I\in{\mathbb{R}}^{m\times m}$ being the identity matrix. We will use the same symbol for a coefficient matrix like $A$ and its Kronecker product $A\otimes I$ as a mapping from the space ${\mathbb{R}}^{sm}$ to itself. In contrast to one-step methods, Peer two-step methods compute $Y_{n}$ from the previous stage vector $Y_{n-1}$ . Hence, also a starting method, given in (9), for the first time interval $[t_{0},t_{1}]$ is required. On each subinterval, Peer methods may be defined by three coefficient matrices $A_{n},B_{n},K_{n}$ , where $A_{n}$ and $K_{n}$ are assumed to be nonsingular. For practical reasons, this general version will not be used, the coefficients in the inner grid points will belong to a fixed Peer method $(A_{n},B_{n},K_{n})\equiv(A,B,K)$ , $n=1,\ldots,N-1$ . The last forward step has the same form as the standard steps but uses exceptional coefficients $(A_{N},B_{N},K_{N})$ to allow for a better approximation of the end conditions.

The KKT conditions (5)-(7) for ODE-constrained optimal control problems on a time interval $[0,T]$ lead to a boundary value problem for a system of two differential equations, see Section 2. The first one corresponds to the original forward ODE for the state solution $y(t)$ and the second one is a linear, adjoint ODE for a Lagrange multiplier $p(t)$ . It is well known that numerical methods for such problems may have to satisfy additional order conditions for the adjoint equation [3, 7, 9, 13, 19, 21]. While these additional conditions are rather mild for one-step methods they may lead to severe restrictions for other types of methods like multistep and Peer methods, especially at the right-hand boundary at the end point $T$ . Here, the order for the approximation of the adjoint equation may be limited to one.

For Peer methods, this question was discussed first in [24] and the adjoint boundary condition at $T$ was identified as the most critical point. In a more recent article [12], these bottlenecks could be circumvented by two measures. First, equivalent formulations of the forward method are not equivalent for the adjoint formulation and using a redundant formulation of Peer methods with three coefficient matrices $(A,B,K)$ adds additional free parameters. The second measure is to consider different methods for the first and last time interval. Hence, instead of one single Peer method (which will be called the standard method) we discuss triplets of Peer methods consisting of a common standard method $(A,B,K)$ for all subintervals of the grid from the interior of $[0,T]$ , plus a starting method $(A_{0},K_{0})$ for the first subinterval and an end method $(A_{N},B_{N},K_{N})$ for the last one. These two boundary methods may have lower order than the standard method since they are used only once.

The present work extends the results from [12] which considered methods with $s=3$ stages only, in two ways. We will now concentrate on methods with four stages and better stability properties like A-stability. The purpose of an accurate solution of the adjoint equation increases the number of conditions on the parameters of the method. Requiring high order of convergence $s$ for the state variable $y(t)$ and order $s-1$ for the adjoint variable $p(t)$ – which we combine to the pair $(s,s-1)$ from now on – a variant of the method BDF3 was identified in [12] as the most attractive standard method. However, this method is not A-stable, with a stability angle of $\alpha=86^{o}$ . In order to obtain A-stability, we will reduce the required orders by one. Still, we will show that convergence with the orders $(s,s-1)$ can be regained by a superconvergence property.

The paper is organized as follows: In Section 2, the boundary value problem arising from the KKT system by eliminating the control, and its discretization by means of discrete adjoint Peer two-step triplets are formulated. An extensive error analysis concentrating on the superconvergence effect is presented in Section 3. The restrictions imposed by the starting and end method on the standard Peer two-step method is studied in Section 4. The following Section 5 describes the actual construction principles of Peer triplets. Numerical tests are done in Section 6. The paper concludes with a summary in Section 7.

2 The Boundary Value Problem

Following the usual Lagrangian approach applied in [12], the first order discrete optimality conditions now consist of the forward equations (9)-(11), the discrete adjoint equations, acting backwards in time,

	$\displaystyle A_{N}^{\sf T}P_{N}=$	$\displaystyle\,w\otimes p_{h}(T)+h\nabla_{Y}F(Y_{N},U_{N})^{\sf T}K_{N}^{\sf T}P_{N},$		(12)
	$\displaystyle A_{n}^{\sf T}P_{n}=$	$\displaystyle\,B_{n+1}^{\sf T}P_{n+1}+h\nabla_{Y}F(Y_{n},U_{n})^{\sf T}K_{n}^{\sf T}P_{n},\ N-1\geq n\geq 0,$		(13)

and the control conditions

	$\displaystyle\,-h\nabla_{U}F(Y_{n},U_{n})^{\sf T}K_{n}^{\sf T}P_{n}\in N_{U^{s}}\big{(}U_{n}\big{)},\ 0\leq n\leq N,$		(14)
	$\displaystyle\,-h\nabla_{u_{0}}f(y_{0},u_{0})^{\sf T}(b^{\sf T}\otimes I)P_{0}\in N_{U}\big{(}u_{0}\big{)}.$		(15)

Here, $p_{h}(T)=\nabla_{y}C\big{(}y_{h}(T)\big{)}^{\sf T}$ and the Jacobians of $F$ are block diagonal matrices $\nabla_{Y}F(Y_{n},U_{n})=\,\mbox{diag}_{i}\big{(}\nabla_{Y_{ni}}f(Y_{ni},U_{ni})\big{)}$ and $\nabla_{U}F(Y_{n},U_{n})=\,\mbox{diag}_{i}\big{(}\nabla_{U_{ni}}f(Y_{ni},U_{ni})\big{)}$ . The generalized normal cone mapping $N_{U^{s}}\big{(}U_{n}\big{)}$ is defined by

\displaystyle N_{U^{s}}(u)=

\displaystyle\,\left\{w\in{\mathbb{R}}^{sd}:w^{T}(v-u)\leq 0\mbox{ for all }v\in U_{ad}^{s}\subset{\mathbb{R}}^{sd}\right\}.

(16)

The discrete KKT conditions (9)-(15) should be good approximations to the continuous ones (5)-(7). In what follows, we assume sufficient smoothness of the optimal control problem such that a local solution $(y^{\star},u^{\star},p^{\star})$ of the KKT system (1)-(3) exists. Furthermore, let the Hamiltonian $H(y,u,p):=p^{\sf T}f(y,u)$ satisfies a coercivity assumption, which is a strong form of a second-order condition. Then the first-order optimality conditions are also sufficient [7]. If $(y,p)$ is sufficiently close to $(y^{\star},p^{\star})$ , the control uniqueness property introduced in [7] yields the existence of a locally unique minimizer $u=u(y,p)$ of the Hamiltonian over all $u\in U_{ad}$ . Substituting $u$ in terms of $(y,p)$ in (5)-(6), gives then the two-point boundary value problem

	$\displaystyle y^{\prime}(t)=$	$\displaystyle\,g\big{(}y(t),p(t)\big{)},\quad y(0)=y_{0},$		(17)
	$\displaystyle p^{\prime}(t)=$	$\displaystyle\,\phi\big{(}y(t),p(t)\big{)},\quad p(T)=\nabla_{y}C\big{(}y(T)\big{)}^{\sf T},$		(18)

with the source functions defined by

\displaystyle g(y,p):=

\displaystyle f\big{(}y,u(y,p)\big{)},\quad\phi(y,p):=-\nabla_{y}f\big{(}y,u(y,p)\big{)}^{\sf T}p.

(19)

The same arguments apply to the discrete first-order optimality system (9)-(15). Substituting the discrete controls $U_{n}=U_{n}(Y_{n},P_{n})$ in terms of $(Y_{n},P_{n})$ and defining

\displaystyle\Phi(Y_{n},K_{n}^{\sf T}P_{n}):=

\displaystyle\,\left(\phi(Y_{ni},(K_{n}^{\sf T}P_{n})_{i})\right)_{i=1}^{s},\quad G(Y_{n},P_{n}):=\left(g(Y_{ni},P_{ni})\right)_{i=1}^{s},

(20)

the approximations for the forward and adjoint differential equations read in a compact form

$\displaystyle A_{0}Y_{0}=$	$\displaystyle\,a\otimes y_{0}+b\otimes g\big{(}y_{0},p_{h}(0)\big{)}+hK_{0}G(Y_{0},P_{0}),$	(21)
$\displaystyle A_{n}Y_{n}=$	$\displaystyle\,B_{n}Y_{n-1}+hK_{n}G(Y_{n},P_{n}),\quad 1\leq n\leq N,$	(22)
$\displaystyle y_{h}(T)=$	$\displaystyle\,(w^{\sf T}\otimes I)Y_{N},$	(23)
$\displaystyle p_{h}(0)=$	$\displaystyle\,(v^{\sf T}\otimes I)P_{0},$	(24)
$\displaystyle A_{n}^{\sf T}P_{n}=$	$\displaystyle\,B_{n+1}^{\sf T}P_{n+1}-h\Phi(Y_{n},K_{n}^{\sf T}P_{n}),\ 0\leq n\leq N-1,$	(25)
$\displaystyle A_{N}^{\sf T}P_{N}=$	$\displaystyle\,w\otimes p_{h}(T)-h\Phi(Y_{N},K_{N}^{\sf T}P_{N}),\ n=N.$	(26)

Here, the value of $p_{h}(0)$ is determined by an interpolant $p_{h}(0)=(v^{\sf T}\otimes I)P_{0}\approx p(0)$ with $v\in{\mathbb{R}}^{s}$ of appropriate order. In a next step, these discrete equations are now treated as a discretization of the two-point boundary value problem (17)-(18). We will derive order conditions and give bounds for the global error.

3 Error Analysis

3.1 Order Conditions

The local error of the standard Peer method and the starting method is easily analyzed by Taylor expansion of the stage residuals, if the exact ODE solutions are used as stages. Hence, defining ${\bf y}_{n}^{(k)}(h{\bf c}):=\big{(}y^{(k)}(t_{n}+hc_{i})\big{)}_{i=1}^{s},\,k=0,1,$ for the forward Peer method, where ${\bf c}=(c_{1},\ldots,c_{s})^{\sf T}$ , local order $q_{1}$ means that

\displaystyle A_{n}{\bf y}_{n}(h{\bf c})-B_{n}{\bf y}_{n-1}(h{\bf c})-hK_{n}{\bf y}^{\prime}(h{\bf c})=O(h^{q_{1}}).

(27)

In all steps of the Peer triplet, requiring local order $q_{1}$ for the state variable and order $q_{2}$ for the adjoint solution leads to the following algebraic conditions from [12]. These conditions depend on the Vandermonde matrix $V_{q}=\big{(}\mbox{1\hskip-2.40005ptl},{\bf c},{\bf c}^{2},\ldots,{\bf c}^{q-1})\in{\mathbb{R}}^{s\times q}$ , the Pascal matrix ${\cal P}_{q}=\Big{(}{j-1\choose i-1}\Big{)}_{i,j=1}^{q}$ and the nilpotent matrix $\tilde{E}_{q}=\big{(}i\delta_{i+1,j}\big{)}_{i,j=1}^{q}$ which commutes with ${\cal P}_{q}=\exp(\tilde{E}_{q})$ . For the different steps (21)-(23) and (24)-(26), and in the same succession we write down the order conditions from [12] when $K_{n}$ is diagonal. The forward conditions are

$\displaystyle A_{0}V_{q_{1}}=$	$\displaystyle ae_{1}^{\sf T}+be_{2}^{\sf T}+K_{0}V_{q_{1}}\tilde{E}_{q_{1}},\ n=0,$	(28)
$\displaystyle A_{n}V_{q_{1}}=$	$\displaystyle B_{n}V_{q_{1}}{\cal P}_{q_{1}}^{-1}+K_{n}V_{q_{1}}\tilde{E}_{q_{1}},\ 1\leq n\leq N,$	(29)
$\displaystyle w^{\sf T}V_{q_{1}}=$	$\displaystyle\mbox{1\hskip-2.40005ptl}^{\sf T},$	(30)

with the cardinal basis vectors $e_{i}\in{\mathbb{R}}^{s}$ , $i=1,\ldots,s$ . The conditions for the adjoint methods are given by

$\displaystyle v^{\sf T}V_{q_{2}}=$	$\displaystyle e_{1}^{\sf T},$	(31)
$\displaystyle A_{n}^{\sf T}V_{q_{2}}=$	$\displaystyle B_{n+1}^{\sf T}V_{q_{2}}{\cal P}_{q_{2}}-K_{n}V_{q_{2}}\tilde{E}_{q_{2}},\ 0\leq n\leq N-1,$	(32)
$\displaystyle A_{N}^{\sf T}V_{q_{2}}=$	$\displaystyle w\mbox{1\hskip-2.40005ptl}^{\sf T}-K_{N}V_{q_{2}}\tilde{E}_{q_{2}},\ n=N.$	(33)

3.2 Bounds for the Global Error

In this section, the errors $\check{Y}_{nj}:=y(t_{nj})-Y_{nj}$ , $\check{P}_{nj}:=p(t_{nj})-P_{nj},$ $n=0,\ldots,N$ , $j=1,\ldots,s$ , are analyzed. According to [12], the equation for the errors $\check{Y}^{\sf T}=(\check{Y}_{0}^{\sf T},\ldots,\check{Y}_{N}^{\sf T})$ and $\check{P}^{\sf T}=(\check{P}_{0}^{\sf T},\ldots,\check{P}_{N}^{\sf T})$ is a linear system of the form

\displaystyle{\mathbb{M}}_{h}\check{Z}=\tau,\ \check{Z}=\begin{pmatrix}\check{Y}\\ \check{P}\end{pmatrix},\ \tau=\begin{pmatrix}\tau^{Y}\\ \tau^{P}\end{pmatrix},

(34)

where the matrix ${\mathbb{M}}_{h}$ has a $2\times 2$ -block structure and $(\tau^{Y},\tau^{P})$ denote the corresponding truncation errors. Deleting all $h$ -depending terms from ${\mathbb{M}}_{h}$ , the block structure of the remaining matrix ${\mathbb{M}}_{0}$ is given by

\displaystyle{\mathbb{M}}_{0}=\begin{pmatrix}M_{11}\otimes I_{m}&0\\ M_{21}\otimes\nabla_{yy}C_{N}&M_{22}\otimes I_{m}\end{pmatrix}

(35)

with $M_{11},M_{21},M_{22}\in{\mathbb{R}}^{s(N+1)\times s(N+1)}$ and a mean value $\nabla_{yy}C_{N}\in{\mathbb{R}}^{m\times m}$ of the symmetric Hessian matrix of $C$ . The index ranges of all three matrices are copied from the numbering of the grid, $0,\ldots,N$ , for convenience. In fact, $M_{21}=(e_{N}\otimes\mbox{1\hskip-2.40005ptl})(e_{N}\otimes w)^{\sf T}$ has rank one only with $e_{N}=(\delta_{Nj})_{j=0}^{N}$ . The diagonal blocks of ${\mathbb{M}}_{0}$ are nonsingular and its inverse has the form

\displaystyle{\mathbb{M}}_{0}^{-1}=\begin{pmatrix}M_{11}^{-1}\otimes I_{m}&0\\ -(M_{22}^{-1}M_{21}M_{11}^{-1})\otimes\nabla_{yy}C_{N}&M_{22}^{-1}\otimes I_{m}\end{pmatrix}.

(36)

The diagonal blocks $M_{11},M_{22}$ have a bi-diagonal block structure with identity matrices $I_{s}$ in the diagonal. The individual $s\times s$ -blocks of their inverses are easily computed with $M_{11}^{-1}$ having lower triangular block form and $M_{22}^{-1}$ upper triangular block form, with blocks

\displaystyle(M_{11}^{-1})_{nk}=\bar{B}_{n}\cdots\bar{B}_{k+1},\,k\leq n,\quad(M_{22}^{-1})_{nk}=\bar{B}_{n+1}^{\sf T}\cdots\bar{B}_{k}^{\sf T},\,k\geq n,

(37)

with the abbreviations $\bar{B}_{n}:=A_{n}^{-1}B_{n}$ , $1\leq n\leq N$ and $\tilde{B}_{n+1}^{\sf T}:=(B_{n+1}A_{n}^{-1})^{\sf T}$ , $0\leq n<N$ . Empty products for $k=n$ mean the identity $I_{s}$ .

Defining ${\mathbb{U}}:=h^{-1}({\mathbb{M}}_{h}-{\mathbb{M}}_{0})$ and rewriting (34) in fixed-point form

\displaystyle\check{Z}=

\displaystyle\,h{\mathbb{M}}_{0}^{-1}{\mathbb{U}}\check{Z}+{\mathbb{M}}_{0}^{-1}\tau\,,

(38)

it has been shown in the proof of Theorem 4.1 of [12] for smooth right hand sides $f$ and $h\leq h_{0}$ that

\displaystyle\|\check{Z}\|\leq 2\max\{\|M_{11}^{-1}\tau^{Y}\|,\|M_{22}^{-1}\tau^{P}\|\}

(39)

in suitable norms, where these norms are discussed in more detail in Section 3.3 below. Moreover, due to the lower triangular block structure of ${\mathbb{M}}_{0}$ the estimate for the error in the state variable may be refined (Lemma 4.2 in [12]) to

\displaystyle\|\check{Y}\|\leq\|M_{11}^{-1}\tau^{Y}\|+hL\|\check{Z}\|

(40)

with some constant $L$ . Without additional conditions, estimates of the terms on the right-hand side of (39) in the form $\|\check{Z}\|=O(h^{-1}\|\tau\|)$ lead to the loss of one order in the global error. However, this loss may be avoided by one additional superconvergence condition on the forward and the adjoint method each, which will be considered next.

3.3 Superconvergence of the Standard Method

For $s$ -stage Peer methods, global order $s$ may be attained in many cases if other properties of the method have lower priority. For optimal stiff stability properties like A-stability, however, it may be necessary to sacrifice one order of consistency as in [17, 23]. Accordingly, in this paper the order conditions for the standard method are lowered by one compared to the requirements in the recent paper [12] to local orders $(s,s-1)$ , see Table 1. Still, the higher global orders may be preserved to some extent by the concept of superconvergence which prevents the order reduction in the global error by cancellation of the leading error term.

Superconvergence is essentially based on the observation that the powers of the forward stability matrix $\bar{B}:=A^{-1}B$ may converge to a rank-one matrix which maps the leading error term of $\tau^{Y}$ to zero. This is the case if the eigenvalue 1 of $\bar{B}$ is isolated. Indeed, if the eigenvalues $\lambda_{j},\,j=1,\ldots,s,$ of the stability matrix $\bar{B}$ satisfy

\displaystyle 1=\lambda_{1}>|\lambda_{2}|\geq\ldots|\lambda_{s}|,

(41)

then its powers $\bar{B}^{n}$ converge to the rank-one matrix $\mbox{1\hskip-2.40005ptl}\mbox{1\hskip-2.40005ptl}^{\sf T}A$ since 1l and $\mbox{1\hskip-2.40005ptl}^{\sf T}A$ are the right and left eigenvectors having unit inner product $\mbox{1\hskip-2.40005ptl}^{\sf T}A\mbox{1\hskip-2.40005ptl}=1$ due to the preconsistency conditions $A^{-1}B\mbox{1\hskip-2.40005ptl}=\mbox{1\hskip-2.40005ptl}$ and $A^{\sf-T}B^{\sf T}\mbox{1\hskip-2.40005ptl}=\mbox{1\hskip-2.40005ptl}$ of the forward and backward standard Peer method, see (29), (32). The convergence is geometric, i.e.

\displaystyle\|\bar{B}^{n}-\mbox{1\hskip-2.40005ptl}\mbox{1\hskip-2.40005ptl}^{\sf T}A\|=\|\big{(}\bar{B}-\mbox{1\hskip-2.40005ptl}\mbox{1\hskip-2.40005ptl}^{\sf T}A\big{)}^{n}\|=O(\gamma^{n})\to 0,\,n\to\infty,

(42)

for any $\gamma\in(|\lambda_{2}|,1)$ . Some care has to be taken here since the error estimate (39) depends on the existence of special norms satisfying $\|\bar{B}\|_{X_{1}}:=\|X_{1}^{-1}\bar{B}X_{1}\|_{\infty}\!=\!1$ , resp. $\|\tilde{B}^{\sf T}\|_{X_{2}}\!=\!1$ . Concentrating on the forward error $\|M_{11}^{-1}\tau^{Y}\|$ , a first transformation of $\bar{B}$ is considered with the matrix

\displaystyle X=\begin{pmatrix}1&-\beta^{\sf T}\\ \mbox{1\hskip-2.40005ptl}_{s-1}&I_{s-1}-\mbox{1\hskip-2.40005ptl}_{s-1}\beta^{\sf T}\end{pmatrix},\ X^{-1}=\begin{pmatrix}\beta_{1}&\beta^{\sf T}\\ -\mbox{1\hskip-2.40005ptl}_{s-1}&I_{s-1}\end{pmatrix},

(43)

where $(\beta_{1},\beta^{\sf T})=\mbox{1\hskip-2.40005ptl}^{\sf T}A$ . Since $Xe_{1}=\mbox{1\hskip-2.40005ptl}$ and $e_{1}^{\sf T}X^{-1}=\mbox{1\hskip-2.40005ptl}^{\sf T}A$ , the matrix $X^{-1}\bar{B}X$ is block-diagonal with the dominating eigenvalue 1 in the first diagonal entry. Due to (41) there exists an additional nonsingular matrix $\Xi\in{\mathbb{R}}^{(s-1)\times(s-1)}$ such that the lower diagonal block of $X^{-1}BX$ has norm smaller one, i.e.

\displaystyle\|\Xi^{-1}\big{(}-\mbox{1\hskip-2.40005ptl}_{s-1},I_{s-1}\big{)}\bar{B}\begin{pmatrix}-\beta^{\sf T}\\ I_{s-1}-\mbox{1\hskip-2.40005ptl}_{s-1}\beta^{\sf T}\end{pmatrix}\Xi\|_{\infty}=\gamma<1.

(44)

Hence, with the matrix

\displaystyle X_{1}:=

\displaystyle\,X\begin{pmatrix}1&\\ &\Xi\end{pmatrix}

(45)

the required norm is found, satisfying $\|\bar{B}\|_{X_{1}}:=\|X_{1}^{-1}\bar{B}X_{1}\|_{\infty}=\max\{1,\gamma\}=1$ and $\|\bar{B}-\mbox{1\hskip-2.40005ptl}\mbox{1\hskip-2.40005ptl}^{\sf T}A\|_{X_{1}}=\gamma<1$ in (42). Using this norm in (39) and (40) for $\epsilon^{Y}:=M_{11}^{-1}\tau^{Y}$ , it is seen with (37) that

\displaystyle\epsilon^{Y}_{n}=

\displaystyle\sum_{k=0}^{n}\mbox{1\hskip-2.40005ptl}(\mbox{1\hskip-2.40005ptl}^{\sf T}A\tau_{n-k}^{Y})+\underbrace{\sum_{k=0}^{n}\big{(}\bar{B}-\mbox{1\hskip-2.40005ptl}\mbox{1\hskip-2.40005ptl}^{\sf T}A\big{)}^{k}\tau_{n-k}^{Y}}_{=O(h^{s})},\ 0\leq n<N.

(46)

Only for $\epsilon_{N}^{Y}$ a slight modification is required and the factors in the second sum have to be replaced by $(\bar{B}_{N}-\mbox{1\hskip-2.40005ptl}\mbox{1\hskip-2.40005ptl}^{\sf T}A)(\bar{B}-\mbox{1\hskip-2.40005ptl}\mbox{1\hskip-2.40005ptl}^{\sf T}A)^{k-1}$ for $k>1$ with norms still of size $O(\gamma^{k})$ . Now, in all cases the loss of one order in the first sum in (46) may be avoided if the leading $O(h^{s})$ -term of $\tau_{n-k}^{Y}$ is canceled in the product with the left eigenvector, i.e. if $\mbox{1\hskip-2.40005ptl}^{\sf T}A\tau_{n-k}^{Y}=O(h^{s+1})$ . An analogous argument may be applied to the second term $\|M_{22}^{-1}\tau^{P}\|$ in (39). The adjoint stability matrix $\tilde{B}^{\sf T}=(BA^{-1})^{\sf T}$ possesses the same eigenvalues as $\bar{B}$ and its leading eigenvectors are also known: $\tilde{B}^{\sf T}\mbox{1\hskip-2.40005ptl}=\mbox{1\hskip-2.40005ptl}$ and $(A\mbox{1\hskip-2.40005ptl})^{\sf T}\tilde{B}^{\sf T}=\mbox{1\hskip-2.40005ptl}^{\sf T}B^{\sf T}=(A\mbox{1\hskip-2.40005ptl})^{\sf T}$ by preconsistency and an analogous construction applies.

Under the conditions corresponding to local orders $(s,s-1)$ the leading error terms in $\tau_{n}^{Y}=\frac{1}{s!}h^{s}\eta_{s}\otimes y^{(s)}(t_{n})+O(h^{s+1})$ and $\tau_{n}^{P}=\frac{1}{(s-1)!}\eta_{s-1}^{\ast}\otimes p^{(s-1)}(t_{n})+O(h^{s})$ are given by

	$\displaystyle\eta_{s}=$	$\displaystyle\,{\bf c}^{s}-A^{-1}B({\bf c}-\mbox{1\hskip-2.40005ptl})^{s}-sA^{-1}K{\bf c}^{s-1},$		(47)
	$\displaystyle\eta_{s-1}^{\ast}=$	$\displaystyle\,{\bf c}^{s-1}-A^{-T}B^{\sf T}({\bf c}+\mbox{1\hskip-2.40005ptl})^{s-1}+(s-1)A^{-T}K{\bf c}^{s-2}.$		(48)

Steps	forward: $q_{1}=s$	adjoint: $q_{2}=s-1$
Start, $n=0$	(28)	(32),(58), $\beta=0$
Standard, $1\leq n<N$	(29)	(32)
Superconvergence	(49)	(50)
Compatibility	(73)	(73)
Last step	(29), $n=N$	(32), $n=N-1$
End point	(30)	(33),(58), $\beta=N$

Table 1: Combined order conditions for the Peer triplet, including the compatibility condition (73) and the condition (58) for full matrices

K_{N}

Considering now $(\mbox{1\hskip-2.40005ptl}^{\sf T}A)\eta_{s}$ in (46) and similarly $(A\mbox{1\hskip-2.40005ptl})^{\sf T}\eta_{s-1}^{\ast}$ the following result is obtained.

Theorem 3.1

Let the Peer triplet with $s>1$ stages satisfy the order conditions collected in Table 1 and let the solutions satisfy $y\in C^{s+1}[0,T]$ , $p\in C^{s}[0,T]$ . Let the coefficients of the standard Peer method satisfy the conditions

	$\displaystyle\mbox{1\hskip-2.40005ptl}^{\sf T}\big{(}A{\bf c}^{s}-B({\bf c}-\mbox{1\hskip-2.40005ptl})^{s}-sK{\bf c}^{s-1}\big{)}=$	$\displaystyle\,0,$		(49)
	$\displaystyle\mbox{1\hskip-2.40005ptl}^{\sf T}\big{(}A^{\sf T}{\bf c}^{s-1}-B^{\sf T}({\bf c}+\mbox{1\hskip-2.40005ptl})^{s-1}+(s-1)K{\bf c}^{s-2}\big{)}=$	$\displaystyle\,0,$		(50)

and let (41) be satisfied. Assume, that a Peer solution $(Y^{\sf T},P^{\sf T})^{\sf T}$ exists and that $f$ and $C$ have bounded second derivatives. Then, for stepsizes $h\leq h_{0}$ the error of these solutions is bounded by

\displaystyle\|Y_{nj}-y(t_{nj})\|_{\infty}+h\|P_{nj}-p(t_{nj})\|_{\infty}=O(h^{s}),

(51)

$n=0,\ldots,N,\,j=1,\ldots,s$ .

Proof: Under condition (49) the representation (46) shows that $\|M_{11}^{-1}\tau^{Y}\|=O(h^{s})$ . In the same way follows $\|M_{22}^{-1}\tau^{P}\|=O(h^{s-1})$ from condition (50). In (39) this leads to a common error $\|\check{Z}\|=O(h^{s-1})$ which may be refined for the state variable with (40) to $\|\check{Y}\|=O(h^{s})$ . $\square$

Remark 3.1

The estimate (46) shows that superconvergence may be a fragile property and may be impaired if $|\lambda_{2}|$ is too close to one, leading to very large values in the bound

\sum_{k=0}^{n}\|\big{(}\bar{B}-\mbox{1\hskip-2.40005ptl}\mbox{1\hskip-2.40005ptl}^{\sf T}A\big{)}^{k}\|_{X_{1}}\leq\frac{1}{1-\gamma}

for the second term in (46). In fact, numerical tests showed that the value $\gamma\cong|\lambda_{2}|$ plays a crucial role. While $|\lambda_{2}|\doteq 0.9$ was appropriate for the 3-stage method AP3o32f which shows order 3 in the tests in Section 6, for two 4-stage methods with $|\lambda_{2}|\doteq 0.9$ superconvergence was not seen in any of our 3 test problems below. Reducing $|\lambda_{2}|$ farther with additional searches we found that a value below $\gamma=0.8$ may be safe to achieve order 4 in practice. Hence, the value $|\lambda_{2}|$ will be one of the data listed in the properties of the Peer methods developed below.

Remark 3.2

By Theorem 3.1, weakening the order conditions from the local order pair $(s+1,s)$ to the present pair $(s,s-1)$ combined with fewer conditions for superconvergence preserves global order $h^{s}$ for the state variable. However, it also leads to a more complicated structure of the leading error. Extending the Taylor expansion for $\tau^{Y},\tau^{P}$ and applying the different bounds, a more detailed representation of the state error may be derived,

	$\displaystyle\\|\check{Y}\\|\leq$	$\displaystyle h^{s}\big{(}\frac{\\|\eta_{s}\\|}{(1-\gamma)s!}\\|y^{(s)}\\|+\frac{\|\mbox{1\hskip-2.40005ptl}^{\sf T}A\eta_{s+1}\|}{(s+1)!}\\|y^{(s+1)}\\|\big{)}$
		$\displaystyle+h^{s}\big{(}\hat{L}\\|p^{(s)}\\|+\hat{L}\\|p^{(s-1)}\\|\big{)}+O(h^{s+1})$		(52)

with some modified constant $\hat{L}$ . Obviously, the leading error depends on four different derivatives of the solutions. Since the dependence on $p$ is rather indirect, we may concentrate on the first line in (3.2). Both derivatives there may not be correlated and it may be difficult to choose a reasonable combination of both error constants as objective function in the construction of efficient methods. Still, in the local error the leading term is obviously $\tau^{Y}\doteq h^{s}\frac{1}{s!}\eta_{s}\otimes y^{(s)}$ , with $\eta_{s}$ defined in (47), and it may be propagated through non-linearity and rounding errors. Hence, in the search for methods,

\displaystyle err_{s}:=\frac{1}{s!}\|\eta_{s}\|_{\infty}=\frac{1}{s!}\|{\bf c}^{s}-A^{-1}B({\bf c}-\mbox{1\hskip-2.40005ptl})^{s}-sA^{-1}K{\bf c}^{s-1}\|_{\infty}

(53)

is used as the leading error constant.

3.4 Adjoint Order Conditions for General Matrices $K_{n}$

The number of order conditions for the boundary methods is so large that they may not be fulfilled with the restriction to diagonal coefficient matrices $K_{0},K_{N}$ for $s\geq 3$ . Hence, it is convenient to make the step to full matrices $K_{0},K_{N}$ in the boundary methods and the order conditions for the adjoint schemes have to be derived for this case. Unfortunately, for such matrices the adjoint schemes (26) and (25) for $n=0$ have a rather unfamiliar form. Luckily, the adjoint differential equation $p^{\prime}=-\nabla_{y}f(y,u)^{\sf T}p$ is linear. We abbreviate the initial value problem for this equation as

\displaystyle p^{\prime}(t)=-J(t)p(t),\quad p(T)=p_{T},

(54)

with $J(t)=\nabla_{y}f\big{(}y(t),u(t)\big{)}^{\sf T}$ and for some boundary index $\beta\in\{0,N\}$ , we consider the matrices $A_{\beta}=(a_{ij}^{(\beta)})$ , $K_{\beta}=(\kappa_{ij}^{(\beta)})$ . Starting the analysis with the simpler end step (26), we have

\displaystyle\sum_{j=1}^{s}a_{ji}^{(N)}P_{Nj}=

\displaystyle\,w_{i}p_{h}(T)+hJ(t_{Ni})\sum_{j=1}^{s}\kappa_{ji}^{(N)}P_{Nj},\ i=1,\ldots,s,

(55)

which is some kind of half-one-leg form since it evaluates the Jacobian $J$ and the solution $p$ at different time points. This step must be analyzed for the linear equation (54) only. Expressions for the higher derivatives of the solution $p$ follow easily:

\displaystyle p^{\prime\prime}=\,

\displaystyle(J^{2}-J^{\prime})p,\quad p^{\prime\prime\prime}=(-J^{3}+2J^{\prime}J+JJ^{\prime}-J^{\prime\prime})p.

(56)

Lemma 3.1

For a smooth coefficient matrix $J(t)$ , the scheme (55) for the linear differential equation (54) has local order 3 under the conditions

\displaystyle A_{N}^{\sf T}V_{3}+K_{N}^{\sf T}V_{3}\tilde{E}_{3}=w\mbox{1\hskip-2.40005ptl}^{\sf T}

(57)

and with $\beta=N$

\displaystyle\sum_{i=1\atop i\not=j}^{s}(c_{i}-c_{j})\kappa_{ij}^{(\beta)}=0,\quad j=1,\ldots,s.

(58)

Proof: Considering the residual of the scheme with the exact solution $p(t)$ , Taylor expansion at $t_{N}$ and the Leibniz rule yield

	$\displaystyle\Delta_{i}=$	$\displaystyle\,\sum_{j=1}^{s}a_{ji}^{(N)}p(t_{Nj})-w_{i}p(T)-hJ(t_{Ni})\sum_{j=1}^{s}\kappa_{ji}^{(N)}p(t_{Nj})$
	$\displaystyle=$	$\displaystyle\,\sum_{k=0}^{q-1}\frac{h^{k}}{k!}\big{(}\sum_{j=1}^{s}a_{ji}^{(N)}c_{j}^{k}-w_{i}\big{)}p^{(k)}-\underbrace{h\sum_{k=0}^{q-2}\frac{h^{k}}{k!}c_{i}^{k}J^{(k)}\sum_{j=1}^{s}\kappa_{ji}^{(N)}\sum_{\ell=0}^{q-2}\frac{h^{\ell}}{\ell!}c_{j}^{\ell}p^{(\ell)}}_{:=\delta}+O(h^{q})$

where all derivatives are evaluated at $t_{N}$ . The second term can be further reformulated as

	$\displaystyle\delta=$	$\displaystyle\,\sum_{k=1}^{q-1}h^{k}\sum_{\ell=0}^{k-1}\sum_{j=1}^{s}\kappa_{ji}^{(N)}\frac{c_{i}^{\ell}c_{j}^{k-\ell-1}}{\ell!(k-\ell-1)!}J^{(\ell)}p^{(k-\ell-1)}$
	$\displaystyle=$	$\displaystyle\,\sum_{k=1}^{q-1}h^{k}\sum_{\ell=1}^{k}\sum_{j=1}^{s}\kappa_{ji}^{(N)}\frac{c_{i}^{\ell-1}c_{j}^{k-\ell}}{(\ell-1)!(k-\ell)!}J^{(\ell-1)}p^{(k-\ell)}.$

Looking at the factors of $h^{0},h^{1},h^{2}$ separately leads to the order conditions

	$\displaystyle h^{0}:\ 0=$	$\displaystyle\,A_{N}^{\sf T}\mbox{1\hskip-2.40005ptl}-w,$
	$\displaystyle h^{1}:\ 0=$	$\displaystyle\,(\sum_{j=1}^{s}a_{ji}^{(N)}c_{j}-w_{i})p^{\prime}-\sum_{j=1}^{s}\kappa_{ji}^{(N)}Jp=(-\sum_{j=1}^{s}a_{ji}^{(N)}c_{j}+w_{i}-\sum_{j=1}^{s}\kappa_{ji}^{(N)})Jp,\mbox{ i.e.}$
	$\displaystyle 0=$	$\displaystyle\,A_{N}^{\sf T}{\bf c}+K_{N}^{\sf T}\mbox{1\hskip-2.40005ptl}-w,$
	$\displaystyle h^{2}:\ 0=$	$\displaystyle\,\frac{1}{2}(\sum_{j}a_{ji}^{(N)}c_{j}^{2}-w_{i})p^{\prime\prime}-\sum_{j=1}^{s}\kappa_{ji}^{(N)}\big{(}c_{j}Jp^{\prime}+c_{i}J^{\prime}p\big{)}$
	$\displaystyle=$	$\displaystyle\,\frac{1}{2}(\sum_{j=1}^{s}c_{j}^{2}a_{ji}^{(N)}-w_{i}+2\sum_{j=1}^{s}c_{j}\kappa_{ji}^{(N)})J^{2}p-\frac{1}{2}(\sum_{j=1}^{s}c_{j}^{2}a_{ji}^{(N)}-w_{i}+2\sum_{j=1}^{s}\kappa_{ji}^{(N)}c_{i})J^{\prime}p.$

Cancellation of the factor of $J^{2}$ requires the condition $0=A_{N}^{\sf T}{\bf c}^{2}+2K_{N}^{\sf T}{\bf c}-w$ , which combines with the $h^{0},h^{1}$ -conditions to (57). The factor of $J^{\prime}$ , however, requires $0=A_{N}^{\sf T}{\bf c}^{2}+2D_{c}K_{N}^{\sf T}\mbox{1\hskip-2.40005ptl}-w$ with $D_{c}=\,\mbox{diag}(c_{i})$ . Subtracting this expression from the factor of $J^{2}$ leaves $0=K_{N}^{\sf T}{\bf c}-D_{c}K_{N}^{\sf T}\mbox{1\hskip-2.40005ptl}$ , which corresponds to (58). $\square$

Remark 3.3

Condition (57) is the standard version for diagonal $K_{N}$ from [12]. Hence, the half-one-leg form of (55) introduces $s$ additional conditions (58), only, for order 3 while the boundary coefficients $K_{\beta}$ , $\beta=0,N$ , may now contain $s(s-1)$ additional elements. In fact, a similar analysis for the step (25) with $n=\beta=0$ reveals again (58) as the only condition in addition to (32).

4 Existence of Boundary Methods Imposes Restrictions on the Standard Method

In the previous paper [12], the combination of forward and adjoint order conditions for the standard method $(A,B,K)$ into one set of equations relating only $A$ and $K$ already gave insight on some background of these methods like the advantages of using symmetric nodes. It also simplifies the actual construction of methods leading to shorter expressions during the elimination process with algebraic software tools. For ease of reference, this singular Sylvester type equation is reproduced here,

\displaystyle(V_{q_{2}}{\cal P}_{q_{2}})^{\sf T}A(V_{q_{1}}{\cal P}_{q_{1}})-V_{q_{2}}^{\sf T}AV_{q_{1}}=(V_{q_{2}}{\cal P}_{q_{2}})^{\sf T}KV_{q_{1}}{\cal P}_{q_{1}}\tilde{E}_{q_{1}}+(V_{q_{2}}\tilde{E}_{q_{2}})^{\sf T}KV_{q_{1}}.

(59)

A similar combination of the order conditions for the boundary methods, however, reveals crucial restrictions: the triplet of methods $(A_{0},K_{0})$ , $(A,B,K)$ , $(A_{N},B_{N},K_{N})$ has to be discussed together since boundary methods of appropriate order may not exist for any standard method $(A,B,K)$ , only for those satisfying certain compatibility conditions required by the boundary methods. Knowing these conditions allows to design the standard method alone without the ballast of two more methods with many additional unknowns. This decoupled construction also greatly reduces the dimension of the search space if methods are optimized by automated search routines. We start with the discussion of the end method.

4.1 Combined Conditions for the End Method

We remind that we now are looking for methods having local order $q_{1}=s$ and $q_{2}=s-1$ everywhere, which we abbreviate from now on as $(q_{1},q_{2})=(s,q)$ . In particular this means that $A^{\sf T}V_{q}+K^{\sf T}V_{q}\tilde{E}_{q}=B^{\sf T}V_{q}{\cal P}_{q}$ for the standard method. Looking for bottlenecks in the design of these methods, we try to identify crucial necessary conditions and consider the three order conditions for the end method $(A_{N},B_{N},K_{N})$ in combination

\displaystyle\begin{array}[]{rl}A_{N}V_{s}-B_{N}V_{s}{\cal P}_{s}^{-1}-K_{N}V_{s}\tilde{E}_{s}&=0,\\[2.84526pt] B_{N}^{\sf T}V_{q}{\cal P}_{q}&=A^{\sf T}V_{q}+K^{\sf T}V_{q}\tilde{E}_{q}=B^{\sf T}V_{q}{\cal P}_{q},\\[2.84526pt] A_{N}^{\sf T}V_{q}+K_{N}^{\sf T}V_{q}\tilde{E}_{q}&=w\mbox{1\hskip-2.40005ptl}_{q}^{\sf T}.\end{array}

(63)

From these conditions the matrices $A_{N},B_{N}$ may be eliminated, revealing the first restrictions on $B$ . Here, the singular matrix map

\displaystyle{\cal L}_{q,s}:\,{\mathbb{R}}^{q\times s}\to{\mathbb{R}}^{q\times s},\quad X\mapsto\tilde{E}_{q}^{\sf T}X+X\tilde{E}_{s}

(64)

plays a crucial role.

Lemma 4.1

A necessary condition for a boundary method $(A_{N},B_{N},K_{N})$ to satisfy (63) is

\displaystyle{\cal L}_{q,s}\big{(}V_{q}^{\sf T}K_{N}V_{s}\big{)}=

\displaystyle\mbox{1\hskip-2.40005ptl}_{q}\mbox{1\hskip-2.40005ptl}_{s}^{\sf T}-V_{q}^{\sf T}BV_{s}{\cal P}_{s}^{-1}.

(65)

Proof: The second condition, $V_{q}^{\sf T}B_{N}=V_{q}^{\sf T}B$ , in (63) leads to the necessary equation

V_{q}^{\sf T}A_{N}V_{s}-V_{q}^{\sf T}K_{N}V_{s}\tilde{E}_{s}=V_{q}^{\sf T}BV_{s}{\cal P}_{s}^{-1}

due to the first condition. The transposed third condition

V_{q}^{\sf T}A_{N}V_{s}+\tilde{E}_{q}^{\sf T}V_{q}^{\sf T}K_{N}V_{s}=\mbox{1\hskip-2.40005ptl}_{q}w^{\sf T}V_{s}=\mbox{1\hskip-2.40005ptl}_{q}\mbox{1\hskip-2.40005ptl}_{s}^{\sf T}

multiplied by the nonsingular matrix $V_{s}$ may be used to eliminate $A_{N}$ and leads to

\tilde{E}_{q}^{\sf T}V_{q}^{\sf T}K_{N}V_{s}+V_{q}^{\sf T}K_{N}V_{s}\tilde{E}_{s}=\mbox{1\hskip-2.40005ptl}_{q}\mbox{1\hskip-2.40005ptl}_{s}^{\sf T}-V_{q}^{\sf T}BV_{s}{\cal P}_{s}^{-1}

which is the equation (65) from the statement. $\square$

Unfortunately, this lemma leads to several restrictions on the design of the methods due to the properties of the map ${\cal L}_{q,s}$ . Firstly, for diagonal matrices $K_{N}$ the image ${\cal L}_{q,s}\big{(}V_{q}^{\sf T}K_{N}V_{s}\big{)}$ has a very restricted shape.

Lemma 4.2

If $K\in{\mathbb{R}}^{s\times s}$ is a diagonal matrix, then $V_{q}^{\sf T}KV_{s}$ and ${\cal L}_{q,s}\big{(}V_{q}^{\sf T}KV_{s}\big{)}$ are Hankel matrices with constant entries along anti-diagonals.

Proof: With $K=\,\mbox{diag}(\kappa_{i})$ , we have $x_{ij}:=e_{i}^{\sf T}(V_{q}^{\sf T}KV_{s})e_{j}=\sum_{k=1}^{s}\kappa_{k}c_{k}^{i+j-2}$ , showing Hankel form of $X=(x_{ij})=:(\xi_{i+j-1})$ for $i=1,\ldots,s-1$ , $j=1,\ldots,s$ . Now, $e_{i}^{\sf T}(\tilde{E}_{q}^{\sf T}X+X\tilde{E}_{s})e_{j}=(i-1)x_{i-1,j}+x_{i,j-1}(j-1)=(i+j-2)\xi_{i+j-2}$ , which shows again Hankel form. $\square$

This lemma means that an end method with diagonal $K_{N}$ only exists if also $V_{q}^{\sf T}BV_{s}{\cal P}_{s}^{-1}$ on the right-hand side of (65) has Hankel structure. Unfortunately it was observed that for standard methods with definite $K$ this is the case for $q_{2}\leq 2$ only (there exist methods with an explicit stage $\kappa_{33}=0$ ).

Remark 4.1

Trying to overcome this bottleneck with diagonal matrices $K_{N}$ , one might consider adding additional stages of the end method. However, using general end nodes $(\hat{c}_{1},\ldots,\hat{c}_{\hat{s}})$ with $\hat{s}\geq s$ does not remove this obstacle. The corresponding matrix ${\cal L}_{q,s}(\hat{V}_{q}^{\sf T}K_{N}\hat{V}_{s})$ with appropriate Vandermonde matrices $\hat{V}$ still has Hankel form.

But even with a full end matrix $K_{N}$ , Lemma 4.1 and the Fredholm alternative enforce restrictions on the standard method $(A,B,K)$ due to the singularity of ${\cal L}_{q,s}$ . This is discussed for the present situation with $s=4,q=3$ , only. The matrix belonging to the map ${\cal L}_{q,s}$ is $I_{s}\otimes\tilde{E}_{q}^{\sf T}+\tilde{E}_{s}^{\sf T}\otimes I_{q}$ and its transpose is $I_{s}\otimes\tilde{E}_{q}+\tilde{E}_{s}\otimes I_{q}$ . Hence, the adjoint of the map ${\cal L}_{q,s}$ is given by

{\cal L}_{q,s}^{\sf T}:\,{\mathbb{R}}^{q\times s}\to{\mathbb{R}}^{q\times s},\quad X\mapsto\tilde{E}_{q}X+X\tilde{E}_{s}^{\sf T}.

Component-wise the map acts as

{\cal L}_{q,s}^{\sf T}:\ (x_{ij})\mapsto\big{(}ix_{i+1,j}+jx_{i,j+1}\big{)}

with elements having indices $i>q$ or $j>s$ being zero. For $q=3,\,s=4,$ one gets

{\cal L}_{3,4}(X)^{\sf T}=\begin{pmatrix}x_{12}+x_{21}&2x_{13}+x_{22}&3x_{14}+x_{23}&x_{24}\\ x_{22}+2x_{31}&2x_{23}+2x_{32}&3x_{24}+2x_{33}&2x_{34}\\ x_{32}&2x_{33}&3x_{34}&0\end{pmatrix}.

It is seen that the kernel of ${\cal L}_{3,4}^{\sf T}$ has dimension 3 and is given by

\displaystyle X=\begin{pmatrix}\xi_{1}&\xi_{2}&\xi_{3}&0\\ -\xi_{2}&-2\xi_{3}&0&0\\ \xi_{3}&0&0&0\end{pmatrix}.

(66)

In (65), the Fredholm condition leads to restrictions on the matrix $B$ from the standard scheme. However, since the matrix $A$ should have triangular form, it is the more natural variable in the search for good methods and an equivalent reformulation of these conditions for $A$ is of practical interest.

Lemma 4.3

Assume that the standard method $(A,B,K)$ has local order $(s,q)=(4,3)$ . Then, end methods $(A_{N},B_{N},K_{N})$ of order $(s,q)=(4,3)$ only exist if the standard method $(A,B,K)$ satisfies the following set of 3 conditions, either for $B$ or for $A$ ,

\displaystyle\left.\begin{array}[]{r}\mbox{1\hskip-2.40005ptl}^{\sf T}B\mbox{1\hskip-2.40005ptl}=1\\[2.84526pt] \mbox{1\hskip-2.40005ptl}^{\sf T}B{\bf c}-{\bf c}^{\sf T}B\mbox{1\hskip-2.40005ptl}=1\\[2.84526pt] \mbox{1\hskip-2.40005ptl}^{\sf T}B{\bf c}^{2}-2{\bf c}^{\sf T}B{\bf c}+({\bf c}^{2})^{\sf T}B\mbox{1\hskip-2.40005ptl}=1\end{array}\right\}\left\{\begin{array}[]{l}1=\mbox{1\hskip-2.40005ptl}^{\sf T}A\mbox{1\hskip-2.40005ptl}\\[2.84526pt] 1=\mbox{1\hskip-2.40005ptl}^{\sf T}A{\bf c}-{\bf c}^{\sf T}A\mbox{1\hskip-2.40005ptl}\\[2.84526pt] 0=\mbox{1\hskip-2.40005ptl}^{\sf T}A{\bf c}^{2}-2{\bf c}^{\sf T}A{\bf c}+({\bf c}^{2})^{\sf T}A\mbox{1\hskip-2.40005ptl}.\end{array}\right.

(73)

Proof: Multiplying equation (65) by ${\cal P}_{s}$ from the right and using $\tilde{E}_{s}{\cal P}_{s}={\cal P}_{s}\tilde{E}_{s}$ , an equivalent form is

{\cal L}_{q,s}(V_{q}^{\sf T}K_{N}V_{s}{\cal P}_{s})=\mbox{1\hskip-2.40005ptl}_{q}\mbox{1\hskip-2.40005ptl}_{s}^{\sf T}{\cal P}_{s}-V_{q}^{\sf T}BV_{s}=:R,

with $\mbox{1\hskip-2.40005ptl}_{s}^{\sf T}{\cal P}_{s}=(1,2,4,8,\ldots)$ by the binomial formula. The Fredholm alternative requires that $tr(X^{\sf T}R)=0$ for all $X$ from (66). We now frequently use the identities $tr\big{(}(vu^{\sf T})R\big{)}=u^{\sf T}Rv$ and $Ve_{i}={\bf c}^{i-1}$ with $e_{i}\in{\mathbb{R}}^{s}$ . The kernel in (66) is spanned by 3 basis elements. The first, $X_{1}=\bar{e}_{1}e_{1}^{\sf T}$ (with the convention $\bar{e}_{i}\in{\mathbb{R}}^{q}$ ) leads to

0\stackrel{{\scriptstyle!}}{{=}}tr(X_{1}^{\sf T}R)=\bar{e}_{1}^{\sf T}\mbox{1\hskip-2.40005ptl}_{q}\mbox{1\hskip-2.40005ptl}_{s}^{\sf T}e_{1}-\bar{e}_{1}^{\sf T}V_{q}^{\sf T}BV_{s}e_{1}=1-\mbox{1\hskip-2.40005ptl}_{s}^{\sf T}B\mbox{1\hskip-2.40005ptl}_{s}.

The second basis element is $X_{2}=\bar{e}_{1}e_{2}^{\sf T}-\bar{e}_{2}e_{1}^{\sf T}$ . Here $tr(X_{2}^{\sf T}\mbox{1\hskip-2.40005ptl}_{q}\mbox{1\hskip-2.40005ptl}_{s}^{\sf T}{\cal P}_{s})=\mbox{1\hskip-2.40005ptl}_{s}^{\sf T}{\cal P}_{s}e_{2}-\mbox{1\hskip-2.40005ptl}_{s}^{\sf T}{\cal P}_{s}e_{1}=1$ and

tr(X_{2}^{\sf T}V_{q}^{\sf T}BV_{s})=\bar{e}_{1}^{\sf T}V_{q}BV_{s}e_{2}-\bar{e}_{2}^{\sf T}V_{q}BV_{s}e_{1}=\mbox{1\hskip-2.40005ptl}_{s}^{\sf T}B{\bf c}-{\bf c}^{\sf T}B\mbox{1\hskip-2.40005ptl}_{s}.

For the third element, $X_{3}=\bar{e}_{1}e_{3}^{\sf T}-2\bar{e}_{2}e_{2}^{\sf T}+\bar{e}_{3}e_{1}^{\sf T}$ , one gets $tr(X_{3}^{\sf T}\mbox{1\hskip-2.40005ptl}_{q}\mbox{1\hskip-2.40005ptl}_{s}^{\sf T}{\cal P}_{s})=\mbox{1\hskip-2.40005ptl}_{s}^{\sf T}{\cal P}_{s}(e_{3}-2e_{2}+e_{1})=4-4+1=1$ . The third condition on $B$ is

tr(X_{3}^{\sf T}V_{q}^{\sf T}BV_{s})=\mbox{1\hskip-2.40005ptl}_{s}^{\sf T}B{\bf c}^{2}-2{\bf c}^{\sf T}B{\bf c}+({\bf c}^{2})^{\sf T}B\mbox{1\hskip-2.40005ptl}_{s}.

The versions for $A$ follow from the order conditions. Let again $\mbox{1\hskip-2.40005ptl}:=\mbox{1\hskip-2.40005ptl}_{s}$ . The first columns of (29) and (32) show $B\mbox{1\hskip-2.40005ptl}=A\mbox{1\hskip-2.40005ptl}$ and $\mbox{1\hskip-2.40005ptl}^{\sf T}B=\mbox{1\hskip-2.40005ptl}^{\sf T}A$ , which gives $\mbox{1\hskip-2.40005ptl}^{\sf T}B\mbox{1\hskip-2.40005ptl}=\mbox{1\hskip-2.40005ptl}^{\sf T}A\mbox{1\hskip-2.40005ptl}=1$ . The second column of (29) reads $B{\bf c}=A{\bf c}+A\mbox{1\hskip-2.40005ptl}-K\mbox{1\hskip-2.40005ptl}$ and leads to $\mbox{1\hskip-2.40005ptl}^{\sf T}A{\bf c}=\mbox{1\hskip-2.40005ptl}^{\sf T}B{\bf c}+\mbox{1\hskip-2.40005ptl}^{\sf T}A\mbox{1\hskip-2.40005ptl}-\mbox{1\hskip-2.40005ptl}^{\sf T}K\mbox{1\hskip-2.40005ptl}$ showing also $\mbox{1\hskip-2.40005ptl}^{\sf T}K\mbox{1\hskip-2.40005ptl}=\mbox{1\hskip-2.40005ptl}^{\sf T}A\mbox{1\hskip-2.40005ptl}=1$ . Hence, the second condition in (73) is equivalent with

1\stackrel{{\scriptstyle!}}{{=}}\mbox{1\hskip-2.40005ptl}^{\sf T}(B{\bf c})-{\bf c}^{\sf T}(B\mbox{1\hskip-2.40005ptl})=\mbox{1\hskip-2.40005ptl}^{\sf T}(A{\bf c}+A\mbox{1\hskip-2.40005ptl}-K\mbox{1\hskip-2.40005ptl})-{\bf c}^{\sf T}A\mbox{1\hskip-2.40005ptl}=\mbox{1\hskip-2.40005ptl}^{\sf T}A{\bf c}-{\bf c}^{\sf T}A\mbox{1\hskip-2.40005ptl}.

In order to show the last equivalence in (73), we have to look ahead at the forward condition (29) for order 3, which is $B{\bf c}^{2}=A({\bf c}^{2}+2{\bf c}+\mbox{1\hskip-2.40005ptl})-2K({\bf c}+\mbox{1\hskip-2.40005ptl})$ . This leads to $\mbox{1\hskip-2.40005ptl}^{\sf T}A{\bf c}^{2}=\mbox{1\hskip-2.40005ptl}^{\sf T}B{\bf c}^{2}=\mbox{1\hskip-2.40005ptl}^{\sf T}A({\bf c}^{2}+2{\bf c}+\mbox{1\hskip-2.40005ptl})-2\mbox{1\hskip-2.40005ptl}^{\sf T}K({\bf c}+\mbox{1\hskip-2.40005ptl})$ , which is equivalent to $2(\mbox{1\hskip-2.40005ptl}^{\sf T}A{\bf c}-\mbox{1\hskip-2.40005ptl}^{\sf T}K{\bf c})=1$ . Now, this expression is required in the last reformulation which also uses the second adjoint order condition $B^{\sf T}{\bf c}=A^{\sf T}{\bf c}-A^{\sf T}\mbox{1\hskip-2.40005ptl}+K^{\sf T}\mbox{1\hskip-2.40005ptl}$ , yielding

	$\displaystyle 1\stackrel{{\scriptstyle!}}{{=}}$	$\displaystyle\,(\mbox{1\hskip-2.40005ptl}^{\sf T}B){\bf c}^{2}+({\bf c}^{2})^{\sf T}(B\mbox{1\hskip-2.40005ptl})-2{\bf c}^{\sf T}(B^{\sf T}{\bf c})$
	$\displaystyle=$	$\displaystyle\,\mbox{1\hskip-2.40005ptl}^{\sf T}A{\bf c}^{2}+({\bf c}^{2})^{\sf T}A\mbox{1\hskip-2.40005ptl}-2{\bf c}^{\sf T}A^{\sf T}{\bf c}+\underbrace{2({\bf c}^{\sf T}A^{\sf T}\mbox{1\hskip-2.40005ptl}-{\bf c}^{\sf T}K^{\sf T}\mbox{1\hskip-2.40005ptl})}_{=1}.\qquad\mbox{$\square$}$

Remark 4.2

It can be shown that only the first condition $\mbox{1\hskip-2.40005ptl}^{\sf T}A\mbox{1\hskip-2.40005ptl}=1$ in (73) is required if the matrix $K$ is diagonal. This first condition is merely a normalization fixing the free common factor in the class $\{\alpha\cdot(A,B,K):\,\alpha\not=0\}$ of equivalent methods. The other two conditions are consequences of the order conditions on the standard method $(A,B,K)$ with diagonal $K$ . However, the proof is rather lengthy and very technical and is omitted.

The restrictions (73) on the standard method also seem to be sufficient with (58) posing no further restrictions. In fact, with (73) the construction of boundary methods was always possible in Section 5.

4.2 Combined Conditions for the Starting Method

The starting method has to satisfy only two conditions

\displaystyle\begin{array}[]{rl}A_{0}V_{s}=&\,ae_{1}^{\sf T}+be_{2}^{\sf T}+K_{0}V_{s}\tilde{E}_{s},\\[2.84526pt] A_{0}^{\sf T}V_{q}=&\,B^{\sf T}V_{q}{\cal P}_{q}-K_{0}^{\sf T}V_{q}\tilde{E}_{q}.\end{array}

(76)

The first two columns of the first equation may be solved for $a=A_{0}V_{s}e_{1}$ and $b=(A_{0}V_{s}-K_{0}V_{s}\tilde{E}_{s})e_{2}$ leading to the reduced conditions

(A_{0}V_{s}-K_{0}V_{s}\tilde{E}_{s})Q_{3}=0,\quad Q_{3}:=I_{s}-e_{1}e_{1}^{\sf T}-e_{2}e_{2}^{\sf T}.

The presence of the projection $Q_{3}$ leads to changes in the condition compared to the end method.

Lemma 4.4

A necessary condition for the starting method $(A_{0},K_{0})$ to satisfy (76) is

\displaystyle\Big{(}{\cal L}_{q,s}\big{(}{\cal P}_{q}^{-T}V_{q}^{\sf T}K_{0}V_{s}\big{)}-V_{q}^{\sf T}BV_{s}\Big{)}Q_{3}=0.

(77)

Proof: Transposing the second condition from (76) and multiplying with $V_{s}Q_{3}$ gives

(V_{q}^{\sf T}A_{0}V_{s}+\tilde{E}_{q}^{\sf T}V_{q}^{\sf T}K_{0}V_{s}-{\cal P}_{q}^{\sf T}V_{q}^{\sf T}BV_{s})Q_{3}=0,

and $V_{q}^{\sf T}A_{0}V_{s}Q_{3}$ may now be eliminated from both equations, yielding

(\tilde{E}_{q}^{\sf T}V_{q}^{\sf T}K_{0}V_{s}+V_{q}^{\sf T}K_{0}V_{s}\tilde{E}_{s})Q_{3}={\cal P}_{q}^{\sf T}V_{q}^{\sf T}BV_{s}Q_{3}.

Again, $P_{q}^{\sf T}$ may be moved to the left side and leads to (77) since it commutes with $\tilde{E}_{q}^{\sf T}$ . $\square$

The situation is now similar to the one for the end method (also concerning a diagonal form of $K_{0}$ ) and we consider again the Fredholm condition. The matrix belonging to the matrix product ${\cal L}_{q,s}()\cdot Q_{3}$ is $Q_{3}\otimes\tilde{E}_{q}^{\sf T}+(\tilde{E}_{s}Q_{3})^{\sf T}\otimes I_{q}$ and it has the transpose $Q_{3}\otimes\tilde{E}_{q}+(\tilde{E}_{s}Q_{3})\otimes I_{q}$ . This matrix belongs to the map

\displaystyle X\mapsto\tilde{E}_{q}(XQ_{3})+(XQ_{3})\tilde{E}_{s}^{\sf T}={\cal L}_{q,s}^{\sf T}(XQ_{3}).

(78)

For $q=3,s=4$ , images of this map are given by

{\cal L}_{3,4}^{\sf T}(XQ_{3})=\begin{pmatrix}0&2x_{13}&x_{23}+3x_{14}&x_{24}\\ 0&2x_{23}&2x_{33}+3x_{24}&2x_{34}\\ 0&2x_{33}&3x_{34}&0\end{pmatrix}.

Here, the map ${\cal L}_{3,4}^{\sf T}$ alone introduces no new kernel elements, the kernel of (78) coincides with that of the map $X\mapsto XQ_{3}$ given by matrices of the form $X=\tilde{X}(I-Q_{3}),\,\tilde{X}\in{\mathbb{R}}^{q\times s}$ . Since the right-hand side of (77) is $V_{q}^{\sf T}BV_{s}Q_{3}$ , the condition for solvability

tr\big{(}X^{\sf T}V_{q}^{\sf T}BV_{s}Q_{3}\big{)}=tr\big{(}\tilde{X}^{\sf T}V_{q}^{\sf T}BV_{s}Q_{3}(I-Q_{3})\big{)}=0

is always satisfied since $Q_{3}(I-Q_{3})=0$ . Hence, no additional restrictions on the standard method are introduced by requiring the existence of starting methods.

5 Construction of Peer Triplets

The construction of Peer triplets requires the solution of the collected order conditions from Table 1 and additional optimization of stability and error properties. However, it has been observed that some of these conditions may be related in non-obvious ways, see e.g. Remark 4.2. This means that the accuracy of numerical solutions may be quite poor due to large and unknown rank deficiencies. Instead, all order conditions were solved here exactly by algebraic manipulation with rational coefficients as far as possible.

The construction of the triplets was simplified by the compatibility conditions (73) allowing the isolated construction of the standard method $(A,B,K)$ without the many additional parameters of the boundary methods. Furthermore, an elimination of the matrix $B$ from forward and adjoint conditions derived in [12], see (59), reduces the number of parameters in $A,K$ to $s(s+3)/2$ elements with $s-1$ additional parameters from the nodes. This is so since $(A,B,K)$ is invariant under a common shift of nodes and we chose the increments $d_{1}=c_{2}-c_{1}$ , $d_{j}=c_{j}-c_{2},\,j=3,\ldots,s$ as parameters. Still, due to the mentioned dependencies between conditions, for $s=4$ a 6-parameter family of methods exists which has been derived explicitly (with quite bulky expressions).

However, optimization of stability properties like $A(\alpha)$ -stability or error constants $err_{s}$ from (53) was not possible in Maple with 6 free parameters. Instead, the algebraic expressions were copied to Matlab scripts for some Monte-Carlo-type search routines. The resulting coefficients of the standard method were finally approximated by rational expressions and brought back to Maple for the construction of the two boundary methods.

At first glance, having the full 6-parameter family of standard methods at hand may seem to be a good work base. However, the large dimension of the search space may prevent optimal results with reasonable effort. This can be seen below, where the restriction to symmetric nodes or singly-implicit methods yielded methods with smaller $err_{4}$ than automated global searches.

$A(\alpha)$ -stability of the method may be checked [12] by considering the eigenvalue problem for the stability matrix $(A-zK)^{-1}B$ , i.e.

\displaystyle Bx=\lambda(A-zK)x\ \iff\ K^{-1}(A-\lambda^{-1}B)x=zx.

(79)

as an eigenvalue problem for $z\in{\mathbb{C}}$ where $|\lambda^{-1}|=1$ runs along the unit circle. Since we focus on A-stable methods, exact verification of this property would have been preferable, of course, but an algebraic proof of A-stability seemed to be out of reach. It turned out that the algebraic criterion of the second author [22] is rather restrictive (often corresponding to norm estimates, Lemma 2.8 ibd). On the other hand, application of the exact Schur criterion [10, 17] is not straight-forward and hardly feasible, since the (rational) coefficients of the stability polynomial are prohibitively large for optimized parameters (dozens of decimal places of the numerators).

5.1 Requirements for the Boundary Methods

Since Lemma 4.3 guarantees the existence of the two boundary methods $(A_{0},K_{0})$ and $(A_{N},B_{N},K_{N})$ , their construction can follow after that of the standard method. Requirements for these two members of the triplet may also be weakened since they are applied once only. This relaxation applies to the order conditions as shown in Table 1, but also to stability requirements. Still, the number of conditions at the boundaries is so large that the diagonal respectively triangular forms of $K_{0},K_{N}$ and $A_{0},A_{N}$ have to be sacrificed and replaced by some triangular block structure. Compared to the computational effort of the complete boundary value problem, the additional complexity of the two boundary steps should not be an issue. However, for non-diagonal matrices $K_{0},K_{N}$ and $s=4$ , the 4 additional one-leg-conditions (58) have to be obeyed.

Weakened stability requirements mean that the last forward Peer step (22) for $n=N$ and the two adjoint Peer steps (25) for $n=N-1$ and $n=0$ need not be A-stable and only nearly zero stable if the corresponding stability matrices have moderate norms. This argument also applies to the two Runge-Kutta steps without solution output (21) and (26). However, the implementation of these steps should be numerically safe for stiff problems and arbitrary step sizes. These steps require the solution of two linear systems with the matrices $A_{0}-hK_{0}J_{0},\ A_{N}^{\sf T}-hJ_{N}^{\sf T}K_{N}^{\sf T}$ or, rather

\displaystyle K_{0}^{-1}A_{0}-hJ_{0},\quad(K_{N}^{-1}A_{N})^{\sf T}-hJ_{N}^{\sf T},

(80)

where $J_{0},J_{N}$ are block diagonal matrices of Jacobians. These Jacobians are expected to have absolutely large eigenvalues in the left complex half-plane. For such eigenvalues, non-singularity of the matrices (80) is assured under the following eigenvalue condition:

\displaystyle\mu_{0}:=\min_{j}\mbox{Re}\lambda_{j}(K_{0}^{-1}A_{0})>0,\quad\mu_{N}:=\min_{j}\mbox{Re}\lambda_{j}(K_{N}^{-1}A_{N})>0.

(81)

In Table 3 the constants $\mu_{0},\mu_{N}$ are displayed for all designed Peer triplets as well as the spectral radii $\varrho(A_{N}^{-1}B_{N})$ and $\varrho(BA_{0}^{-1}),\,\varrho(B_{N}A^{-1})$ for the boundary steps of the Peer triplet. In the search for the boundary methods with their exact algebraic parameterizations, these spectral radii were minimized, and if they were close to one the values $\mu_{0},\mu_{N}>0$ were maximized.

5.2 $A(\alpha)$ -Stable 4-Stage Methods

Although our focus is on A-stable methods, we also shortly consider $A(\alpha)$ -stable methods. We like to consider BDF4 (backward differentiation formulas) as a benchmark, since triplets based on BDF3 were the most efficient ones in [12]. In order to distinguish the different methods, we denote them in the form AP $s$ o $q_{1}q_{2}$ aaa, where AP stands for Adjoint Peer method followed by the stage number and the smallest forward and adjoint orders $q_{1}$ and $q_{2}$ in the triplet. The trailing letters are related to properties of the method like diagonal or singly diagonal implicitness.

The Peer triplet AP4o43bdf based on BDF4 has equi-spaced nodes $c_{i}=i/4,\,i=1,\ldots,4$ , yielding $w=e_{4}$ . The coefficients of the full triplet are given in Appendix A.1. Obviously the method is singly-implicit and its well-known stability angle is $\alpha=73.35^{o}$ . We also monitor the norm of the zero-stability matrix $\|A^{-1}B\|_{\infty}$ , which may be a measure for the propagation of rounding errors. Its value is $\|A^{-1}B\|_{\infty}\leq 5.80$ . Since BDF4 has full global order 4, the error constant from (53) is $err_{4}=0$ . Still, the end methods were constructed with the local orders $(4,3)$ according to Table 1. The matrices of the corresponding starting method have a leading $3\times 3$ block and a separated last stage. We abbreviate this as block structure $(3+1)$ . All characteristics of the boundary method are given in Table 3.

Motivated by the beneficial properties of Peer methods with symmetric nodes seen in [12, 24], the nodes of the next triplet with diagonally-implicit standard method were chosen symmetric to a common center, i.e. $c_{1}+c_{4}=c_{2}+c_{3}$ . Unfortunately, searches for large stability angles with such nodes in the interval $[0,1]$ did not find A-stable methods, but the following method AP4o43dif with flip symmetric nodes and $\alpha=84.00^{o}$ , which is an improvement of 10 degrees over BDF4. Its coefficients are given in Appendix A.2. Although there exist A-stable methods with other nodes in $[0,1]$ , this method is of its own interest since its error constant $err_{4}\doteq 2.5$ e-3 is surprisingly small. The node vector of AP4o43dif includes $c_{4}=1$ , leading again to $w=e_{4}$ . Further properties of the standard method $(A,B,K)$ are $\|A^{-1}B\|_{\infty}\leq 2.01$ and the damping factor $|\lambda_{2}|=0.26$ . See Table 3 for the boundary methods.

5.3 A-stable Methods

By extensive searches with the full 6-parameter family of diagonally-implicit 4-stage methods many A-stable methods were found even with nodes in $[0,1]$ . In fact, regions with A-stable methods exist in at least 3 of the 8 octants in $(d_{1},d_{3},d_{4})$ -space. Surprisingly, however, for none of these methods the last node $c_{4}$ was the rightmost one. Also, it may be unexpected that some of the diagonal elements of $A$ and $K$ have negative signs. This does not cause stability problems if $a_{ii}\kappa_{ii}>0,\,i=1,\ldots,s$ , see also (80). A-stability assured, the search procedure tried to minimize a linear combination $err_{s}+\delta\|A^{-1}B\|_{\infty}$ of the error constant and the norm of the stability matrix with small $\delta<10^{-3}$ to account for the different magnitudes of these data. As mentioned in Remark 3.1, it was also necessary to include the damping factor $|\lambda_{2}|$ in the minimization process. One of the best A-stable standard methods with general nodes was named AP4o43dig. Its coefficients are given in Appendix A.3. It has an error constant $err_{4}\leq 0.0260$ and damping factor $|\lambda_{2}|\doteq 0.798$ , see Table 2. No block structure was chosen in the boundary methods in order to avoid large norms of the mixed stability matrices $BA_{0}^{-1},B_{N}A^{-1},A_{N}^{-1}B_{N}$ , see Table 3.

It was observed that properties of the methods may improve, if the nodes have a wider spread than the standard interval $[0,1]$ . In our setting, the general vector $w$ allows for an end evaluation $y_{h}(T)$ at any place between the nodes. Since an evaluation roughly in the middle of all nodes may have good properties, in a further search the nodes were restricted to the interval $[0,2]$ . Indeed, all characteristic data of the method AP4o43die with extended nodes presented in Appendix A.4 have improved. As mentioned before, the standard method is invariant under a common node-shift and a nearly minimal error constant was obtained with the node increments $d_{1}=\frac{10}{11}$ , $d_{3}=-1$ and $d_{4}=\frac{2}{3}$ . Then, the additional freedom in the choice of $c_{2}$ was needed for the boundary methods, since the conditions (81) could only be satisfied in a small interval around $c_{2}=\frac{5}{4}$ . The full node vector with this choice has alternating node increments since $c_{3}<c_{1}<c_{2}<c_{4}$ . The method is A-stable, its error constant $err_{4}\leq 0.0136$ is almost half as large as for the method AP4o43dig, and $\|A^{-1}B\|_{\infty}\leq 6.1$ and $|\lambda_{2}|\leq\frac{2}{3}$ are smaller, too. The data of the boundary methods can be found in Table 3.

$s$	triplet	nodes	$\alpha$	$\\|A^{-1}B\\|_{\infty}$	$\|\lambda_{2}\|$	$err_{s}$	Remarks
4	AP4o43bdf	BDF4	$73.35^{o}$	5.79	0.099	0	singly-implicit
	AP4o43dif	$[0,1]$	$84.0^{o}$	$2.01$	0.26	0.0025	diag.-implicit
	AP4o43dig	$[0,1]$	$90^{o}$	24.5	0.798	0.0260	$c_{3}$ =1
	AP4o43sil	$[0,1]$	$90^{o}$	32.2	0.60	0.0230	$c_{3}$ =1, sing.impl.
	AP4o43die	$[0,2]$	$90^{o}$	6.08	0.66	0.0135	nodes alternate
3	AP3o32f	$[0,1]$	$90^{o}$	15.3	0.91	0.0170	nodes alternate

Table 2: Properties of the standard methods of Peer triplets.

	Starting method			End method
triplet	blocks	$\mu_{0}$	$\varrho(BA_{0}^{-1})$	blocks	$\mu_{N}$	$\varrho(A_{N}^{-1}B_{N})$	$\varrho(B_{N}A^{-1})$
AP4o43bdf	3+1	5.47	1	1+3	3.81	1	1.15
AP4o43dif	3+1	6.27	1	1+3	4.40	1	1.03
AP4o43dig	4	0.99	1	4	0.89	1.001	1
AP4o43sil	3+1	1.88	1	4	0.72	1	1.03
AP4o43die	3+1	3.80	1	1+3	0.66	2.6	1.98
AP3o32f	1+1+1	1.50	1.02	1+2	0.94	1	1

Table 3: Properties of the boundary methods of Peer triplets.

For medium-sized ODE problems, where direct solvers for the stage equations may be used, an additional helpful feature is diagonal singly-implicitness of the standard method. In our context this means that the triangular matrices $K^{-1}A$ and $AK^{-1}$ have a constant value $\theta$ in the main diagonal. Using the ansatz

a_{ii}=\theta\kappa_{ii},\ i=1,\ldots,s,

for $A=(a_{ij})$ and $K=(\kappa_{ij})$ , the order conditions from Table 1 lead to a cubic equation for $\theta$ with no rational solutions, in general. In order to avoid pollution of the algebraic elimination through superfluous terms caused by rounding errors, numerical solutions for this cubic equations were not used until the very end. This means that also the order conditions for the boundary methods were solved with $\theta$ as a free parameter, only the final evaluation of the coefficients in Appendix A.5 was done with its numerical value. Also in the Matlab-search for A-stable methods, the cubic equation was solved numerically and it turned out that the largest positive solution gave the best properties. Hence, this Peer triplet was named AP4o43sil. For a first candidate with nearly minimal $err_{s}+\delta\|A^{-1}B\|_{\infty}$ , the damping factor $\gamma\doteq 0.89$ was again too close to one to ensure superconvergence in numerical tests, see Remark 3.1. However, further searches nearby minimizing the damping factor found a better standard method with $|\lambda_{2}|=0.60$ . Its nodes are ${\bf c}^{\sf T}=\left(\frac{1}{50},\frac{3}{5},1,\frac{41}{85}\right)$ , the diagonal parameter $\theta=3.34552931287687520$ is the largest zero of the cubic equation

112673616\,\theta^{3}+106686908\,\theta^{2}-2102637319\,\theta+1621264295=0.

Further properties of the standard method are A-stability, nodes in $(0,1]$ with $c_{3}=1$ , norm $\|A^{-1}B\|_{\infty}=32.2$ , and an error constant $err_{4}=0.0230$ . The end method $(A_{N},B_{N},K_{N})$ has full matrices $A_{N},K_{N}$ , see Table 3.

For the sake of completeness, we also present an A-stable diagonally-implicit 3-stage method, since in the previous paper [12] we could not find such methods with reasonable parameters. After relaxing the order conditions by using superconvergence, such methods exist. Applying all conditions for forward order $s=3$ and adjoint order $s-1=2$ , there remains a 5-parameter family depending on the node differences $d_{1}=c_{2}-c_{1}$ , $d_{3}=c_{3}-c_{2}$ and 3 elements of $A$ or $K$ . A-stable methods exist in all 4 corners of the square $[-\frac{1}{2},\frac{1}{2}]^{2}$ in the $(d_{1},d_{3})$ -plane, the smallest errors $err_{3}$ were observed in the second quadrant. The method AP3o32f with $(d_{1},d_{3})=(-\frac{5}{27},\frac{2}{5})$ has a node vector with $c_{3}=1$ . The coefficients can be found in Appendix A.6. The characteristic data are $err_{3}\leq 0.017$ , $\|A^{-1}B\|_{\infty}\leq 15.3$ , $|\lambda_{2}|=0.91$ . The starting method is of standard form with lower triangular $A_{0}$ and diagonal $K_{0}$ .

The main properties of the newly developed Peer triplets are summarized in Table 2 for the standard methods and Table 3 for the boundary methods.

6 Numerical Results

We present numerical results for all Peer triplets listed in Table 2 and compare them with those obtained for the third-order four-stage one-step W-method Ros3wo proposed in [13] which is also A-stable. All calculations have been done with Matlab-Version R2019a, using the nonlinear solver fsolve to approximate the overall coupled scheme (21)–(26) with a tolerance $1e\!-\!14$ . To illustrate the rates of convergence, we consider two nonlinear optimal control problems, the Rayleigh and a controlled motion problem. A linear wave problem is studied to demonstrate the practical importance of A-stability for larger time steps.

6.1 The Rayleigh Problem

The first problem is taken from [11] and describes the behaviour of a tunnel-diode oscillator. With the electric current $y_{1}(t)$ and the transformed voltage at the generator $u(t)$ , the Rayleigh problem reads

$\displaystyle\mbox{Minimize }\int_{0}^{2.5}\left(u(t)^{2}+y_{1}(t)^{2}\right)\,dt$		(82)
$\displaystyle\mbox{subject to }y^{\prime\prime}_{1}(t)-y^{\prime}_{1}\left(1.4-0.14y^{\prime}_{1}(t)^{2}\right)+y_{1}(t)=$	$\displaystyle\,4u(t),\quad t\in(0,2.5],$	(83)
$\displaystyle y_{1}(0)=y^{\prime}_{1}(0)=$	$\displaystyle\,-\!5.$	(84)

Introducing $y_{2}(t)=y_{1}^{\prime}(t)$ and eliminating the control $u(t)$ yields the following nonlinear boundary value problem (see [13] for more details):

$\displaystyle y^{\prime}_{1}(t)=$	$\displaystyle\,y_{2}(t),$	(85)
$\displaystyle y^{\prime}_{2}(t)=$	$\displaystyle\,-y_{1}(t)+y_{2}\left(1.4-0.14y_{2}(t)^{2}\right)-8p_{2}(t),$	(86)
	$\displaystyle\,y_{1}(0)=-5,\,y_{2}(0)=-5,$	(87)
$\displaystyle p^{\prime}_{1}(t)=$	$\displaystyle\,p_{2}(t)-2y_{1}(t),$	(88)
$\displaystyle p^{\prime}_{2}(t)=$	$\displaystyle\,-p_{1}(t)-(1.4-0.42y_{2}(t)^{2})p_{2}(t),$	(89)
	$\displaystyle\,p_{1}(2.5)=0,\,p_{2}(2.5)=0.$	(90)

Refer to caption — Figure 1: Rayleigh Problem: Convergence of the maximal state errors $\|(w^{\sf T}\otimes I)Y_{n}-y(t_{n+1})\|_{\infty}$ (left) and adjoint errors $\|(v^{\sf T}\otimes I)P_{n}-p(t_{n})\|_{\infty}$ (right), $n=0,\ldots,N$ .

To study convergence orders of our new methods, we compute a reference solution at the discrete time points $t=t_{n}$ by applying the classical fourth-order RK4 with $N\!=\!1280$ steps. Numerical results for the maximum state and adjoint errors are presented in Figure 1 for $N\!+\!1\!=\!20,40,80,160,320$ . AP3o32f and Ros3wo show their expected orders $(3,2)$ and $(3,3)$ for state and adjoint solutions, respectively. Order three for the adjoint solutions is achieved by all new four-stage Peer methods. The smaller error constants of AP4o43bdf, AP4o43dif and Ros3wo are clearly visible. The additional superconvergence order four for the state solutions shows up for AP4o43die and AP4o43sil and nearly for AP4o43dif and AP4o43dig. AP4o43bdf does not reach its full order four here, too.

6.2 A Controlled Motion Problem

This problem was studied in [14]. The motion of a damped oscillator is controlled in a double-well potential, where the control $u(t)$ acts only on the velocity $y_{2}(t)$ . The optimal control problem reads

$\displaystyle\mbox{Minimize }\frac{\alpha}{2}\\|y(6)-y_{f}\\|^{2}+\frac{1}{2}\int_{0}^{6}u(t)^{2}\,dt$		(91)
$\displaystyle\mbox{subject to }y^{\prime}_{1}(t)-y_{2}(t)=$	$\displaystyle\,0,$	(92)
$\displaystyle y^{\prime}_{2}(t)-y_{1}(t)+y_{1}(t)^{3}+\nu y_{2}(t)=$	$\displaystyle\,u(t),\quad t\in(0,6],$	(93)
$\displaystyle y_{1}(0)=-1,\,y_{2}(0)=$	$\displaystyle\,0,$	(94)

where $\nu>0$ is the damping parameter and $y_{f}$ the target final position. As in [14], we set $\nu=1$ , $y_{f}=(1,0)^{\sf T}$ , and $\alpha=10$ .

Eliminating the scalar control $u(t)$ yields the following nonlinear boundary value problem:

$\displaystyle y^{\prime}_{1}(t)=$	$\displaystyle\,y_{2}(t),$	(95)
$\displaystyle y^{\prime}_{2}(t)=$	$\displaystyle\,y_{1}(t)-y_{1}(t)^{3}-\nu y_{2}(t)-p_{2}(t),$	(96)
	$\displaystyle\,y_{1}(0)=-1,\,y_{2}(0)=0,$	(97)
$\displaystyle p^{\prime}_{1}(t)=$	$\displaystyle\,(3y_{1}(t)^{2}-1)p_{2}(t),$	(98)
$\displaystyle p^{\prime}_{2}(t)=$	$\displaystyle\,-p_{1}(t)+\nu p_{2}(t),$	(99)
	$\displaystyle\,p_{1}(6)=\alpha(y_{1}(6)-1),\,p_{2}(6)=\alpha y_{2}(6).$	(100)

The optimal control $u^{\star}=-p_{2}^{\star}$ must accelerate the motion of the particle to follow an optimal path $(y_{1}^{\star},y_{2}^{\star})$ through the total energy field $E=\frac{1}{2}y_{2}^{2}+\frac{1}{4}y_{1}^{4}-\frac{1}{2}y_{1}^{2}$ , shown in Figure 2 on the top, in order to reach the final target $y_{f}$ behind the saddle point. The cost obtained from a reference solution with $N\!=\!1279$ is $C(y(6))=0.77674$ , which is in good agreement with the lower order approximation in [14]. Numerical results for the maximum state and adjoint errors are presented in Figure 2 for $N\!+\!1\!=\!10,20,40,80,160$ . Worthy of mentioning is the repeated excellent performance of AP4o43bdf and AP4o43dif, but also the convincing results achieved by the third-order method Ros3wo. All theoretical orders are well observable, except for AP4o43dig, which tends to order three for the state solutions. A closer inspection reveals that this is caused by the second state $y_{2}$ , while the first one asymptotically converges with fourth order. However, the three methods AP4o43die, AP4o43dig, and AP4o43sil perform quite similar. Observe that AP3o32f has convergence problems for $N\!=\!9$ .

The stagnation of the state errors for the finest step sizes is due to the limited accuracy of Matlab’s fsolve – a fact which was already reported in [9].

6.3 A Wave Problem

The third problem is taken from [4] and demonstrates the practical importance of A-stability. We consider the optimal control problem

$\displaystyle\mbox{Minimize }y_{1}(1)+\frac{1}{2}\,\int_{0}^{1}u(t)^{2}\,dt$		(101)
$\displaystyle\mbox{subject to }y^{\prime\prime}_{1}(t)+(2\pi\kappa)^{2}y_{1}(t)=$	$\displaystyle\,u(t),\quad t\in(0,1],$	(102)
$\displaystyle y_{1}(0)=y^{\prime}_{1}(0)=$	$\displaystyle\,0,$	(103)

where $\kappa=16$ is used. Introducing $y_{2}(t)=y_{1}^{\prime}(t)$ and eliminating the control $u(t)$ yields the following linear boundary value problem:

$\displaystyle y^{\prime}_{1}(t)=$	$\displaystyle\,y_{2}(t),$	(104)
$\displaystyle y^{\prime}_{2}(t)=$	$\displaystyle\,-(2\pi\kappa)^{2}y_{1}(t)-p_{2}(t),$	(105)
	$\displaystyle\,y_{1}(0)=0,\,y_{2}(0)=0,$	(106)
$\displaystyle p^{\prime}_{1}(t)=$	$\displaystyle\,(2\pi\kappa)^{2}p_{2}(t),$	(107)
$\displaystyle p^{\prime}_{2}(t)=$	$\displaystyle\,-p_{1}(t),$	(108)
	$\displaystyle\,p_{1}(1)=1,\,p_{2}(1)=0.$	(109)

The exact solutions are given by

$\displaystyle y_{1}^{*}(t)=$	$\displaystyle\,\frac{1}{2(2\pi\kappa)^{3}}\sin(2\pi\kappa t)-\frac{t}{2(2\pi\kappa)^{2}}\cos(2\pi\kappa t),$	(110)
$\displaystyle y_{2}^{*}(t)=$	$\displaystyle\,\frac{t}{2(2\pi\kappa)}\sin(2\pi\kappa t),$	(111)
$\displaystyle p_{1}^{*}(t)=$	$\displaystyle\,\cos(2\pi\kappa t),\;p_{2}^{*}(t)=\frac{1}{2\pi\kappa}\sin(2\pi\kappa t),$	(112)

and the optimal control is $u^{*}(t)=-p_{2}^{*}(t)$ . The key observation here is that the eigenvalues of the Jacobian of the right-hand side in (104)-(109) are $\lambda_{1/2}=2\pi\kappa i$ and $\lambda_{3/4}=-2\pi\kappa i$ , which requests appropriate step sizes for the only $A(\alpha$ )-stable methods AP4o43bdf and AP4o34diw due to their stability restrictions along the imaginary axis. Indeed, a closer inspection of the stability region of the (multistep) BDF4 method near the origin reveals that $|\lambda h_{bdf4}|\leq 0.3$ is a minimum requirement to achieve acceptable approximations for problems with imaginary eigenvalues and moderate time horizon. For the four-stage AP4o43bdf with step size $h=4h_{bdf4}$ , this yields $|\lambda h|\leq 1.2$ and hence $h\leq 1/(32\pi)\approx 0.02$ for the wave problem considered. A similar argument applies to AP4o34diw, too. Numerical results for the maximum state and adjoint errors are plotted in Figure 3 for $N\!+\!1\!=\!20,40,80,160,320$ . They clearly show that both methods deliver first feasible results for $h=1/80$ and below only, but then again outperform the other Peer methods. Once again, Ros3wo performs remarkably well. The orders of convergence for the adjoint solutions are one better than the theoretical values, possibly due to the overall linear structure of the boundary value problem.

7 Summary

We have extended our three-stage adjoint Peer two-step methods constructed in [12] to four stages to not only improve the convergence order of the methods but also their stability. Combining superconvergence of a standard Peer method with a careful design of starting and end Peer methods with appropriately enhanced structure, discrete adjoint A-stable Peer methods of order $(4,3)$ could be found. Still, new requirements had to be dealt with for the higher order pair (4,3). The property of A-stability comes with larger error constants and a few other, minor structural disadvantages. As long as A-stability is not an issue to solve the boundary value problem arising from eliminating the control from the system of KKT conditions, a Peer variant AP4o43bdf of the $A(73.35^{o})$ -stable BDF4 is the most attractive method, closely followed by the $A(84^{o})$ -stable AP4o43dif, which performs equally well in nearly all of our test cases. The A-stable methods AP4o43dig and AP4o43die with diagonally implicit standard Peer methods are very good alternatives if eigenvalues close or on the imaginary axis are existent. We have also constructed the A-stable method AP4o43sil with singly-diagonal main part as an additional option if large linear systems can be still solved by a direct solver and hence the property of requesting one LU-decomposition only is highly valuable. In future work, we plan to train our novel methods in a projected gradient approach to also tackle large-scale PDE constrained optimal control problems with semi-discretizations in space. In these applications, Peer triplets may have to satisfy even more severe requirements.

Acknowledgements. The first author is supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) within the collaborative research center TRR154 “Mathematical modeling, simulation and optimisation using the example of gas networks” (Project-ID 239904186, TRR154/2-2018, TP B01).

References

[1] G. Albi, M. Herty, and L. Pareschi. Linear multistep methods for optimal control problems and applications to hyperbolic relaxation systems. Applied Mathematics and Computation, 354:460–477, 2019.
[2] I. Almuslimani and G. Vilmart. Explicit stabilized integrators for stiff optimal control problems. SIAM J. Sci. Comput., 43:A721–A743, 2021.
[3] F.J. Bonnans and J. Laurent-Varin. Computation of order conditions for symplectic partitioned Runge-Kutta schemes with application to optimal control. Numer. Math., 103:1–10, 2006.
[4] A.L. Dontchev, W.W. Hager, and V.M. Veliov. Second-order Runge-Kutta approximations in control constrained optimal control. SIAM J. Numer. Anal., 38:202–226, 2000.
[5] A. Gerisch, J. Lang, H. Podhaisky, and R. Weiner. High-order linearly implicit two-step peer – finite element methods for time-dependent PDEs. Appl. Numer. Math., 59:624–638, 2009.
[6] B. Gottermeier and J. Lang. Adaptive two-step peer methods for incompressible Navier-Stokes equations. In G. Kreiss, P. Lötstedt, A. Malqvist, and M. Neytcheva, editors, Proceedings of ENUMATH 2009, the 8th European Conference on Numerical Mathematics and Advanced Applications, Uppsala, July 2009, Numerical Mathematics and Advanced Applications, pages 387–395. Springer, 2009.
[7] W.W. Hager. Runge-Kutta methods in optimal control and the transformed adjoint system. Numer. Math., 87:247–282, 2000.
[8] E. Hairer, G. Wanner, and Ch. Lubich. Geometric Numerical Integration, Structure-Preserving Algorithms for Ordinary Differential Equations, volume 31 of Springer Series in Computational Mathematic. Springer, Heidelberg, Berlin, 1970.
[9] M. Herty, L. Pareschi, and S. Steffensen. Implicit-explicit Runge-Kutta schemes for numerical discretization of optimal control problems. SIAM J. Numer. Anal., 51:1875–1899, 2013.
[10] Z. Jackiewicz. General Linear Methods for Ordinary Differential Equations. John Wiley&Sons, Hoboken, New Jersey, 2009.
[11] D.H. Jacobson and D.Q. Mayne. Differential Dynamic Programming. American Elsevier Publishing, New York, 1970.
[12] J. Lang and B.A. Schmitt. Discrete adjoint implicit peer methods in optimal control. Technical Report arXiv:2002.12081, 2020.
[13] J. Lang and J.G. Verwer. W-methods in optimal control. Numer. Math., 124:337–360, 2013.
[14] X. Liu and J. Frank. Symplectic Runge-Kutta discretization of a regularized forward-backward sweep iteration for optimal control problems. J. Comput. Appl. Math., 383:113133, 2021.
[15] Ch. Lubich and A. Ostermann. Runge-Kutta approximation of quasi-linear parabolic equations. Math. Comp., 64:601–627, 1995.
[16] T. Matsuda and Y. Miyatake. Generalization of partitioned Runge–Kutta methods for adjoint systems. J. Comput. Appl. Math., 388:113308, 2021.
[17] J.I. Montijano, H. Podhaisky, L. Randez, and M. Calvo. A family of L-stable singly implicit Peer methods for solving stiff IVPs. BIT, 59:483–502, 2019.
[18] A. Ostermann and M. Roche. Runge-Kutta methods for partial differential equations and fractional orders of convergence. Math. Comp., 59:403–420, 1992.
[19] A. Sandu. On the properties of Runge-Kutta discrete adjoints. Lecture Notes in Computer Science, 3394:550–557, 2006.
[20] A. Sandu. Reverse automatic differentiation of linear multistep methods. In C. Bischof, H. Bücker, P. Hovland, U. Naumann, and J. Utke, editors, Advances in Automatic Differentiation, volume 64 of Lecture Notes in Computational Science and Engineering, pages 1–12. Springer, Berlin, 2008.
[21] J.M. Sanz-Serna. Symplectic Runge–Kutta schemes for adjoint equations, automatic differentiation, optimal control, and more. SIAM Review, 58:3–33, 2016.
[22] B.A. Schmitt. Algebraic criteria for A-stability of peer two-step methods. Technical Report arXiv:1506.05738, 2015.
[23] B.A. Schmitt and R. Weiner. Efficient A-stable peer two-step methods. J. Comput. Appl. Math., pages 319–32, 2017.
[24] D. Schröder, J. Lang, and R. Weiner. Stability and consistency of discrete adjoint implicit peer methods. J. Comput. Appl. Math., 262:73–86, 2014.
[25] J.L. Troutman. Variational Calculus and Optimal Control. Springer, New York, 1996.

Appendix

In what follows, we will give the coefficient matrices which define the Peer triplets discussed above. We have used the symbolic option in Maple as long as possible to avoid any roundoff errors which would pollute the symbolic manipulations by a great number of superfluous terms. If possible, we provide exact rational numbers for the coefficients and give numbers with $16$ digits otherwise. It is sufficient to only show pairs $(A_{n},K_{n})$ and the node vector ${\bf c}$ , since all other parameters can be easily computed from the following relations:

	$\displaystyle B_{n}=$	$\displaystyle\,(A_{n}V_{s}-K_{n}V_{s}\tilde{E}_{s})P_{s}V_{s}^{-1},$
	$\displaystyle a=$	$\displaystyle\,A_{0}\mbox{1\hskip-2.40005ptl},\;b=A_{0}c-K_{0}\mbox{1\hskip-2.40005ptl},\;w=V_{s}^{\sf-T}\mbox{1\hskip-2.40005ptl},\;v=V_{s}^{\sf-T}e_{1},\quad s=3,4,$

with $e_{1}=(1,0,\ldots,0)^{\sf T}\in{\mathbb{R}}^{s}$ and the special matrices

\displaystyle V_{s}=\big{(}\mbox{1\hskip-2.40005ptl},{\bf c},{\bf c}^{2},\ldots,{\bf c}^{s-1}),\quad{\cal P}_{q}=\Big{(}{j-1\choose i-1}\Big{)}_{i,j=1}^{s},\quad\tilde{E}_{s}=\big{(}i\delta_{i+1,j}\big{)}_{i,j=1}^{s}\,.

A.1. Coefficients of AP4o54bdf

{\bf c}^{\sf T}=\left(\frac{1}{4},\frac{1}{2},\frac{3}{4},1\right)

\displaystyle A_{0}=\begin{pmatrix}2&\frac{1}{2}&&\\[5.69054pt] -\frac{265}{96}&\frac{17}{96}&\frac{11}{288}&\\[5.69054pt] \frac{7}{6}&-\frac{47}{24}&\frac{25}{12}&\\[5.69054pt] -\frac{21}{32}&\frac{227}{96}&-\frac{1163}{288}&\frac{25}{12}\end{pmatrix},\,K_{0}=\begin{pmatrix}\frac{1}{2}&&&\\[5.69054pt] -\frac{77}{192}&\frac{3}{32}&&\\[5.69054pt] \frac{67}{192}&\frac{17}{96}&\frac{155}{576}&\\[5.69054pt] -\frac{19}{192}&-\frac{17}{192}&0&\frac{1}{4}\end{pmatrix}

\displaystyle A=\begin{pmatrix}\frac{25}{12}\\[2.84526pt] -4&\frac{25}{12}&\\[2.84526pt] 3&-4&\frac{25}{12}&\\[2.84526pt] -\frac{4}{3}&3&-4&\frac{25}{12}\end{pmatrix},\,K=\begin{pmatrix}\frac{1}{4}&&&\\[5.69054pt] &\frac{1}{4}&&\\[5.69054pt] &&\frac{1}{4}&\\[5.69054pt] &&&\frac{1}{4}\end{pmatrix}

\displaystyle A_{N}=\begin{pmatrix}\frac{635}{96}&&&\\[5.69054pt] -\frac{1235}{72}&\frac{35}{32}&\frac{67}{96}&-\frac{43}{288}\\[5.69054pt] \frac{4475}{288}&-\frac{35}{24}&0&\frac{43}{72}\\[5.69054pt] -5&\frac{35}{96}&-\frac{67}{96}&\frac{53}{96}\end{pmatrix},\,K_{N}=\begin{pmatrix}\frac{25}{32}&&&\\[5.69054pt] -\frac{5}{3}&\frac{61}{192}&-\frac{1}{192}&\\[5.69054pt] \frac{115}{64}&-\frac{13}{48}&\frac{23}{64}&\\[5.69054pt] -\frac{185}{288}&\frac{13}{96}&-\frac{1}{192}&\frac{43}{576}\end{pmatrix}

A.2. Coefficients of AP4o43dif

{\bf c}^{\sf T}=\left(\frac{3}{22},\frac{53}{132},\frac{97}{132},1\right)

\displaystyle\small A_{0}=\begin{pmatrix}\hfill 1.1582197171362010&\hfill 0.04624378638947835&&\\[2.84526pt] -0.7020381998219871&\hfill 1.55936795391810600&\hfill 0.02624560436867219&\\[2.84526pt] -0.5084723270832399&-3.71639374160988500&\hfill 2.21790664582949000&\\[2.84526pt] &\hfill 2.30412222276248700&-2.66055187655540200&1.275350214666080\end{pmatrix}

\displaystyle\small K_{0}=\begin{pmatrix}\hfill 0.10630513138050800&&&\\[2.84526pt] \hfill 0.39188176473763570&\hfill 0.18135555683856680&&\\[2.84526pt] -0.05671611789968874&\hfill 0.06216151000641002&0.3773797141792473&\\[2.84526pt] -0.08101130220826174&-0.03462160050989925&&0.1300000000000000\end{pmatrix}

\displaystyle\small A=\begin{pmatrix}\hfill 2.713996187194519&&&\\[2.84526pt] -5.753019558612675&\hfill 2.063116456071261&&\\[2.84526pt] \hfill 5.300000000000000&-4.392801539381829&\hfill 2.202673482804081&\\[2.84526pt] -2.313267438350870&\hfill 2.523025304770754&-2.619073109161321&1.275350214666080\end{pmatrix}

\displaystyle\small K=\,\mbox{diag}\left(0.2212740342685062,0.2910929443629617,0.3576330213685321,0.13\right)

\displaystyle\small A_{N}=\begin{pmatrix}\hfill 3.321208926131899&&&\\[2.84526pt] -6.825690771130220&\hfill 1.096545539465253&\hfill 0.5537291752348234&-0.08772005153993665\\[2.84526pt] \hfill 5.504481844998321&-1.589672639210835&\hfill 0.3213889335444567&\hfill 0.44690680951897520\\[2.84526pt] -2.000000000000000&\hfill 0.493127099745582&-0.8751181087792801&\hfill 0.64081324202096150\end{pmatrix}

\displaystyle\small K_{N}=\begin{pmatrix}\hfill 0.4780945554021703&&&\\[2.84526pt] -0.7500000000000000&\hfill 0.3443968810148793&0.03064999258349816&\\[2.84526pt] \hfill 0.9263584006629529&-0.2474620802680682&0.34743387239297130&\\[2.84526pt] -0.4116869618629235&\hfill 0.1378269814151266&0.03853141924782626&0.06599889592052376\end{pmatrix}

A.3. Coefficients of AP4o43dig

{\bf c}^{\sf T}=\left(\frac{139}{1159},\frac{11}{19},1,\frac{1375}{2014}\right)

\displaystyle\small A_{0}=\begin{pmatrix}-482.1874750642102&4.750000000000000&-5.916666666666667&-6.500000000000000\\[5.69054pt] 5295.612100386801&78.73229010791468&60.16394904407432&7.222222222222222\\[5.69054pt] 893.8003010294580&-4.061724320422766&9.736228361340601&19.67879886925837\\[5.69054pt] -5707.694317901957&-68.66254446706468&-58.05689396607474&-35.61626763467349\end{pmatrix}

\displaystyle\small K_{0}=\begin{pmatrix}-49.91295086094522&0.5250000000000000&-3.439024390243902&2.894736842105263\\[5.69054pt] 405.5730073881453&49.31516975831240&8.193548387096774&-3.428571428571429\\[5.69054pt] 53.62032171809015&12.67084977396168&-1.304859285573478&4.013292871986014\\[5.69054pt] -414.6382351541371&-49.08870633833948&-1.334271117642095&-15.23643896150272\end{pmatrix}

\displaystyle\small A=\begin{pmatrix}-2.604429828805958&&&\\[5.69054pt] 6.603320924494022&11.44234275562775&&\\[5.69054pt] 0.5317173544040980&-2.710438820206414&3.550000000000000&\\[5.69054pt] -5.000000000000000&2.026117385005894&2.376616772673509&-15.21524654319290\end{pmatrix}

\displaystyle\small K=\,\mbox{diag}\left(-0.8973222553064913,3.337407156628221,1.164566261468968,-2.604651162790698\right)

\displaystyle\small A_{N}=\begin{pmatrix}-3.754385964912281&0.01222493887530562&1.014925373134328&-0.1403508771929825\\[5.69054pt] 11.35280296428295e&45.64990373363066&-15.17493010383148&-20.79910559459065\\[5.69054pt] 0.03205698176794144&2.937595714981687&-3.756123160591242&2.666624105937791\\[5.69054pt] -7.630473981138614&-48.59972438748765&18.91612789128839&18.27283236584584\end{pmatrix}

\displaystyle\small K_{N}=\begin{pmatrix}-0.9578456075353955&-0.3387096774193548&0.1194029850746269&0.8045112781954887\\[5.69054pt] 11.57142857142857&-96.60667975350700&-20.12500000000000&102.4846028390543\\[5.69054pt] 3.761888534906390&-33.44512959236514&-5.865106921633463&34.94704625261853\\[5.69054pt] -15.32045224098695&134.2026156894527&26.37615509001594&-141.3196101415549\end{pmatrix}

A.4. Coefficients of AP4o43die

{\bf c}^{\sf T}=\left(\frac{15}{44},\frac{5}{4},\frac{1}{4},\frac{23}{12}\right)

\displaystyle A_{0}=\begin{pmatrix}\frac{1573}{27}&&&\\[5.69054pt] -\frac{29140496694667}{1728480384000}&\frac{37071572404007}{69715375488000}&\frac{37246788257}{8450348544}&\\[5.69054pt] -\frac{98934973036237}{2160600480000}&-\frac{59311823623513}{87144219360000}&-\frac{51770824817}{13203669600}&\\[5.69054pt] \frac{7653678714559}{1387052160000}&-\frac{7872573544487}{38730764160000}&-\frac{32557703329}{23473190400}&\frac{18014543}{144484600}\end{pmatrix}

\displaystyle K_{0}=\begin{pmatrix}\frac{110}{9}&&&\\[5.69054pt] \frac{4}{5}&\frac{168095353644187}{920242956441600}&-\frac{136571614975979}{36809718257664}&\\[5.69054pt] -\frac{9857504559041}{1080300240000}&\frac{398472962076949}{11503036955520000}&-\frac{18235836500357}{18404859128832}&\\[5.69054pt] -\frac{1423069729157}{1440400320000}&\frac{398472962076949}{7668691303680000}&\frac{136571614975979}{61349530429440}&\frac{2}{61}\end{pmatrix}

\displaystyle A=\begin{pmatrix}\frac{45808744223}{19505421000}&&&\\[5.69054pt] -\frac{279428522}{187552125}&\frac{3}{4}&&\\[5.69054pt] \frac{285647}{15004170}&\frac{22704013}{125034750}&-\frac{11}{10}&\\[5.69054pt] \frac{1}{4}&-\frac{11824391}{41678250}&\frac{832579}{4167825}&\frac{18014543}{144484600}\end{pmatrix}

\displaystyle K=\,\mbox{diag}\left(\frac{35085281}{25006950},\frac{2300653}{8335650},-\frac{1780019}{2500695},\frac{2}{61}\right)

\displaystyle A_{N}=\begin{pmatrix}\frac{46268635184890481}{6231747944448000}&&&\\[5.69054pt] -\frac{843924159892681}{239682613248000}&\frac{2095498352743}{2535104563200}&\frac{240331128931}{253510456320}&-\frac{10279031063}{122060590080}\\[5.69054pt] -4&\frac{3733202770769}{25351045632000}&-\frac{8883867896017}{6337761408000}&-\frac{192302368031}{24412118016000}\\[5.69054pt] \frac{39222471881369}{27696657530880}&-\frac{637146509711}{2816782848000}&-\frac{191242568381}{352097856000}&\frac{175009899277}{8137372672000}\end{pmatrix}

\displaystyle K_{N}=\begin{pmatrix}\frac{95205971609617}{35952391987200}&&&\\[5.69054pt] &\frac{28298016708823}{167316901171200}&-\frac{1296366536717}{4182922529280}&\\[5.69054pt] -\frac{958020525197}{864240192000}&-\frac{2425439590003}{104573063232000}&-\frac{27877023310129}{41829225292800}&\\[5.69054pt] -\frac{958020525197}{14980163328000}&-\frac{2425439590003}{69715375488000}&\frac{1296366536717}{6971537548800}&-\frac{22310489177}{4882423603200}\end{pmatrix}

A.5. Coefficients of AP4o43sil

{\bf c}^{\sf T}=\left(\frac{1}{50},\frac{3}{5},1,\frac{41}{85}\right)

\displaystyle\small A_{0}=\begin{pmatrix}-18.6770976012982273&-1.15212718448036531&-0.684527356670693701&\\ 30.2098963703001422&-19.0677876392318276&-7.55433120044842482&\\ -9.81986015015644262&-2.15227598175855777&4.86425591259034856&\\[5.69054pt] -\frac{57}{25}&\frac{695}{72}&8.28572795643498617&9.37534909694128499\end{pmatrix}

\displaystyle\small K_{0}=\begin{pmatrix}-11.4061014637853601&-0.0776818914116313719&-0.278650826939386227&\\[5.69054pt] \frac{59}{28}&-13.9738118565057040&1.13881886074868390&\\[5.69054pt] -\frac{2498819}{583100}&2.75133184568842863&0.477663277652266390&\\[5.69054pt] \frac{161}{25}&\frac{779}{80}&-0.352459713213735910&2.80235150260251923\end{pmatrix}

\displaystyle\small A=\begin{pmatrix}-3.40824065799546119&&&\\ -10.5240959029253065&-6.21116392867196304&&\\ 1.24215119761892880&-3.47742608697035235&3.58958480260299081&\\ 12.1231239821473188&-3.03082301205062472&1.32154050930321039&9.37534909694123776\end{pmatrix}

\displaystyle\small K=\,\mbox{diag}\left(-1.01874482010281778,-1.85655642136068493,1.07294973886097804,2.80235150260250919\right)

\displaystyle\small A_{N}=\begin{pmatrix}-3.93487127199009570&\frac{75}{26}&-\frac{7}{8}&\\[5.69054pt] -4.71932036763822580&42.7928650670737006&-3.60582429888170880&-36.9240613682294166\\[5.69054pt] -\frac{71}{202}&-2.31371287972058488&0.191639567731266976&1.90723457480523846\\[5.69054pt] 9.00567678814317298&-43.3637675719685003&5.28918473115044182&35.0168267934241781\end{pmatrix}

\displaystyle\small K_{N}=\begin{pmatrix}-0.687420439097535868&&&\\[5.69054pt] \frac{247}{72}&24.5972883813771674&-2.02081177860650990&-24.7095335501766993\\[5.69054pt] -0.427115478996059306&-\frac{1780}{289}&0.897376604330415086&5.61580307958561348\\[5.69054pt] -3.39815952896990200&-\frac{356}{17}&1.56153637437775765&22.4504633222483055\end{pmatrix}

A.6. Coefficients of AP3o32f

{\bf c}^{\sf T}=\left(\frac{106}{135},\frac{3}{5},1\right)

\displaystyle A_{0}=\begin{pmatrix}-\frac{13474483}{2809000}&&\\[5.69054pt] \frac{2765681}{1404500}&\frac{753641}{273375}&\\[5.69054pt] -\frac{48583191}{81461000}&-\frac{1538339}{1093500}&\frac{1783}{580}\end{pmatrix}

\displaystyle K_{0}=\,\mbox{diag}\left(-\frac{13474483}{7155000},\frac{2513302}{1366875},\frac{11}{10}\right)

\displaystyle A=\begin{pmatrix}-\frac{11}{2}&&\\[5.69054pt] \frac{6493}{2700}&\frac{64}{25}&\\[5.69054pt] -\frac{25757}{78300}&-\frac{121}{100}&\frac{1783}{580}\end{pmatrix}

\displaystyle K=\,\mbox{diag}\left(-\frac{93}{50},\frac{44}{25},\frac{11}{10}\right)

\displaystyle A_{N}=\begin{pmatrix}-3&&\\[5.69054pt] -\frac{559409}{391500}&\frac{5418793}{1458000}&\frac{2257039}{1691280}\\[5.69054pt] \frac{1733909}{391500}&-\frac{5418793}{1458000}&-\frac{565759}{1691280}\end{pmatrix}

\displaystyle K_{N}=\,\mbox{diag}\left(-\frac{1190159}{978750},\frac{5418793}{3645000},\frac{2257039}{4228200}\right)

	$\displaystyle\\|\check{Y}\\|\leq$	$\displaystyle h^{s}\big{(}\frac{\\|\eta_{s}\\|}{(1-\gamma)s!}\\|y^{(s)}\\|+\frac{\|\mbox{1\hskip-2.40005ptl}^{\sf T}A\eta_{s+1}\|}{(s+1)!}\\|y^{(s+1)}\\|\big{)}$
		$\displaystyle+h^{s}\big{(}\hat{L}\\|p^{(s)}\\|+\hat{L}\\|p^{(s-1)}\\|\big{)}+O(h^{s+1})$		(52)