\section{Introduction}
Strong stability preserving Runge--Kutta (RK) methods were first introduced by
Shu and Osher \cite{shu1988} in the context of time integration for first-order hyperbolic 
conservation laws:
\begin{align} \label{pde}
	{\mathcal U}_t + {\mathcal F}({\mathcal U})_x & = 0, &
    {\mathcal U}(x,t=0) = {\mathcal U}_0.
\end{align}
After semi-discretization, \eqref{pde} takes the form of an initial-value 
ordinary differential equation system:
\begin{align}
u'(t) & = f(u) & u(0) = u_0\, , \label{ode}
\end{align}
where $f$ is a discrete approximation to $-\mathcal F$.
In the scalar case $\mathcal U$ is dissipative, and it is natural to seek a
semi-discretization that is dissipative: 
\begin{align} \label{dissipative}
\frac{d}{dt}\|u\| \le 0,
\end{align}
where $\|\cdot\|$ denotes a convex functional (e.g., a norm, a semi--norm, \ldots). 
This is achieved by biasing the discretization $f$ in the upwind direction.
A necessary condition for \eqref{dissipative} is monotonicity under an
explicit Euler step \cite[p. 501]{Kraaijevanger1991}:
\begin{align} \label{FEcond1}
	\| v + h f(v)\| & \le \|v\|, & \mbox{for all $v$, and for $h$ satisfying } 0 \le h \le h_0,
\end{align}
where $h_0>0$ (in general $h_0$ may depend on $v$).  
Let $u_{n}, u_{n+1}$ denote approximations, computed by some numerical integrator,
to the solution at successive time steps $t_{n}$ and $t_{n+1}=t_{n}+h$.
Under the forward Euler monotonicity condition \eqref{FEcond1}, 
it is possible to prove 
that many Runge--Kutta and linear multistep methods also give monotone
solutions; i.e., solutions that satisfy
\begin{align}
	\|u_{n+1}\| & \le \|u_{n}\| & \mbox{for $h$ satisfying } 0 \le h \le R h_0.
\end{align}
Such methods are known as  strong stability preserving (SSP) methods, and
the factor $R$ is known as the radius of absolute monotonicity or SSP coefficient
of the method.
SSP methods necessarily have non-negative coefficients, since the monotonicity
property is proved using \eqref{FEcond1} and convexity.  

Monotonicity cannot be ensured using only assumption \eqref{FEcond1} in general for methods with negative coefficients \cite[Thm. 4.2]{Kraaijevanger1991}, or even for some methods 
(such as the classical fourth-order RK method) with non-negative coefficients \cite[Thm. 9.6]{Kraaijevanger1991}.
In order to accommodate such methods, a second discrete approximation to $-\mathcal F$
is introduced and referred to as $\tilde{f}$.  This discretization must be monotone
under an explicit Euler step with negative step size:
\begin{align} \label{FEcond2}
	\| v - h f(v)\| & \le \|v\|, & \mbox{for all $v$, and for $h$ satisfying } 
    										0 \le h \le \tilde{h}_0, 
\end{align}
where $\tilde{h}_0>0$. In the context of hyperbolic problems, $\tilde{f}$ must be biased in the downwind direction,
and typically $\tilde{h}_0 = h_0$. 
The downwind spatial discretization $\tilde{f}$ is to be used in place of $f$ wherever 
a negative coefficient  appears in the time integration method, in order to ensure 
monotonicity of the overall method.  Introduction of a downwind discretization
makes it possible to ensure monotonicity for a broader class of methods,
including the classical RK method of order four.  It also makes it possible
to ensure monotonicity for many methods under larger step sizes.

Methods that use both upwind and downwind operators can naturally be viewed as
{\em perturbed Runge--Kutta methods}. Although they are also connected to
additive RK methods (see \cite{higueras2005a,higueras2006}),
in the present work we will employ the perturbation viewpoint,
and refer to methods that use downwind discretization as perturbed RK methods.

During the last quarter century, a number of additional authors have studied 
monotonicity for methods that use downwind discretization.  
The main motivation for this work has been to
break the  ``order barrier'' that restricts explicit RK methods to order
four, or to find new methods with larger SSP coefficients.
In this context, numerical optimization of the SSP coefficient for RK methods with negative
coefficients was conducted for explicit methods in \cite{ruuth2004,ruuth2006,gottlieb2006} 
and for implicit methods in \cite{ketcheson2011a}.  
In each case, optimization was
carried out over methods with a specified order and number of stages.

The present work stems from a different motivation.  Monotonicity preservation
is not the only numerical property of interest in applications, and practitioners
may wish to use a particular integrator that has small or zero SSP coefficient.  
Our goal is then to perturb the prescribed method in order to achieve
larger monotonicity-preserving timesteps.
Little work has been done in this direction, because it is not known
how to find the best perturbation for a given method.  That problem 
is the main focus of this work.  
Most of our results concern only explicit methods,
although one major result (Theorem \ref{thm:feasible}) pertains
also to implicit methods.  

\subsection{Perturbed Runge--Kutta methods}
A Runge-Kutta method applied to the initial value problem \eqref{ode}
computes approximations $u_n \approx u(t_n)$ by
\begin{subequations} \label{Runge--Kutta}
\begin{align}
    Y & = u_n e + h K F\, ,  \\
    u_{n+1} & = Y_{s+1}.
\end{align}
\end{subequations}
Here $s$ is the number of stages, $e$ is a vector whose entries are equal to one, $Y$ is the vector containing the stage values and the numerical solution, $Y=(Y_1, \ldots, Y_s, Y_{s+1})^t$,  
$[F]_i = f(Y_i)$, and $K$ is the $(s+1)\times(s+1)$ matrix of Butcher
coefficients:
$$K=\begin{pmatrix} A & 0 \\ b^t & 0 \end{pmatrix}\, . 
$$
In this work we study perturbed Runge--Kutta methods:
\begin{subequations} \label{PRunge--Kutta}
\begin{align}
    Y & = u_n e + h K F + h\widetilde{\bbA}(F-\widetilde{F})\, ,  \\
    u_{n+1} & = Y_{s+1},
\end{align}
\end{subequations}
where $\widetilde{K}$ is given by
$$\widetilde{\bbA}=\begin{pmatrix} \tilde A & 0 \\ \tilde b^t & 0 \end{pmatrix}\, , 
$$
and $\widetilde{F}$ is defined analogously to $F$.  
We assume that the
perturbation $\tilde A$ has the same structure
(strictly lower-triangular, lower-triangular, or full)
as the matrix $A$.

Observe that the perturbed method \eqref{PRunge--Kutta} reduces to the RK method 
\eqref{Runge--Kutta} when $\widetilde{F} = F$.
Method \eqref{PRunge--Kutta} may be viewed as approximating the solution 
of the perturbed problem
$$u'(t) = f(u) + (f(u) - \tilde{f}(u)),$$
where $\tilde{f}\approx f$. 

We assume that $f$ and $\tilde{f}$ satisfy the explicit Euler assumptions 
\eqref{FEcond1} and \eqref{FEcond2}, respectively, with $\tilde{h}_0 = h_0$.

Most previous works, including \cite{ruuth2004,ruuth2006}, have focused
on methods with the following property
\begin{quote}
	{Property C}: 
    
    We say that a perturbation $\widetilde{K}$ to an RK method $K$ possesses
    {\em property C} if, for each value of $j$
    \begin{align}
    	\widetilde{K}_{ij}   \ne 0 \mbox{ (for some $i$)} \quad \implies  \quad
                  K_{ij}   = 0 \mbox{ (for all $i$).}
    \end{align}
\end{quote}
In words, property C means that in the $j$th column, only one of $K, \widetilde{K}$
has any nonzero entries.  Thus, only one of $f(y_j), \tilde{f}(y_j)$
need ever be evaluated, so only $s$ total function evaluations are required per step.
In \cite{gottlieb2006} it was shown that
for WENO discretizations, the cost of computing both
$f(y_j)$ and $\tilde{f}(y_j)$ is much less than twice the cost of computing $f(y_j)$
alone.  Therefore methods that without property C may also be of practical
interest.  In the present work, we do not assume property C.



\subsection{Scope and outline}
The central question of the present work is
\begin{itemize}
    \item Given a fixed RK method, what perturbation results in the largest
    		radius of absolute monotonicity for the perturbed method?
\end{itemize}
We investigate this question in the context of both linear problems
(Section \ref{sec:lin:prob}) and nonlinear problems (Section \ref{sec:nonlinear}).

For explicit methods applied to linear problems, the question above can be cast in terms of
absolute monotonicity of the (bivariate) stability polynomial.  In Section 
\ref{sec:upperbound:lin} we prove a general upper bound on the radius
of absolute monotonicity of the stability polynomial of an explicit perturbed RK method
with $s$ stages and linear order $p$.
 In Section \ref{sec:upperlin:num}, we provide an algorithm for
computing tighter bounds, and tabulate some of the resulting numerical values. 
Examples of optimal methods are given in \ref{sec:lin-examples}.

Section \ref{sec:nonlinear} is devoted to perturbed Runge--Kutta methods for
nonlinear problems.  
Theorem \ref{thm:feasible} states that a perturbation with positive
radius of absolute monotonicity exists for every Runge--Kutta method.
Theorems \ref{prop_vpos} and \ref{Pbound} give simple upper bounds on the optimal
perturbed radius for a given explicit method.
In Section \ref{sec:algorithms} we give two algorithms for computing optimal
perturbations.  The first is provably correct but approximate, while the second is heuristic but exact and agrees with the first in all cases we have tested.  Both are applicable only to explicit methods.
These algorithms have been implemented in the software package
Nodepy \cite{nodepy}.
We conclude Section \ref{sec:nonlinear} with an application of the theorems and
algorithms to optimal perturbations of some RK methods from the literature.
Among the results is the first truly optimal perturbation for the classical
4th-order method of Kutta.

Section \ref{sec:conclusions} contains some conclusions as well as some open questions to be studied in the future.

Finally, in the Appendix we give some details on perturbations for the family of second order 2-stage methods and the classical fourth order RK method.


\section{Explicit perturbed Runge--Kutta methods for linear problems}\label{sec:lin:prob}

To study the behavior of the perturbed Runge--Kutta method \eqref{PRunge--Kutta}
for linear problems, we apply it to a linear scalar test problem,
setting $f(u)=\lambda u$ and $\tilde{f}(u)=\tilde{\lambda} u$ in \eqref{PRunge--Kutta}.
This results in the iteration $$u_{n+1} = \phi_{(\bbA, \Apert)}(z, -\tilde{z}) \, u_n,$$
where $z=h\lambda$, $\tilde{z}=h\tilde{\lambda}$ and
\begin{align}
\phi_{(\bbA, \Apert)}(z, \tilde{z})= 1 +  \left(z b^t + (z+\tilde{z}) \tilde{b}^t\right)\left(I -  z A - (z+\tilde{z}) \tilde{A} \right)^{-1} e\, . \label{bistabfunc}
\end{align}
We refer to \eqref{bistabfunc} as the {\it stability function of the perturbed
Runge--Kutta method \eqref{PRunge--Kutta}}.

We say that a function $\psi:\mathbb{R}\to \mathbb{R}$  is absolutely monotonic (a.m.) at
$\xi$ if all derivatives at $\xi$ exist and they are non-negative. 
For a function $\psi:\mathbb{R}^2\to \mathbb{R}$ the concept of absolute
monotonicity can be defined in a similar way: $\psi$ is a.m.
at $(\xi, \tilde \xi)$ if all derivatives at $(\xi, \tilde \xi)$ exist and they
are non-negative. 



Given a function $\psi(z, \tilde{z})$, we define the {\it radius of absolute monotonicity}  as 
\begin{align}
R(\psi)=  \sup\left\{ r\in \mathbb{R} \, | \,  r=0, \hbox{ or }  r>0\, , \, \hbox{and }
\psi(z, \tilde{z}) \text{ is a.m. at } (-r, -r) \right\}\, . \label{threshold:pert}
\end{align}

For a perturbed Runge--Kutta method \eqref{PRunge--Kutta} with coefficients $(K, \widetilde{\bbA})$,
we write ${R}_{\hbox{\footnotesize Lin}}(K, \widetilde{\bbA})$ to denote the radius of absolute monotonicity of its stability
function:
$$ {R}_{\hbox{\footnotesize Lin}}(K, \widetilde{\bbA})=R(\phi_{(\bbA, \Apert)})\, .  $$
The quantity ${R}_{\hbox{\footnotesize Lin}}(K,\widetilde{\bbA})$ is referred to as the {\it threshold factor} due to its
role in the step size for monotonicity.
For a given Runge--Kutta method \eqref{Runge--Kutta} with coefficients $K$,
we are interested in determining perturbations $\widetilde{\bbA}$ that give the largest 
threshold factor.  The corresponding {\it  threshold factor of the optimal perturbation} is
denoted by
\begin{align} \label{opt-thresh-fact}
  {R}^\textup{opt}_{\hbox{\footnotesize Lin}}(K)=  \sup_{\widetilde{\bbA}} {R}_{\hbox{\footnotesize Lin}}(K, \widetilde{\bbA})\, .
\end{align}
The supremum in \eqref{opt-thresh-fact} is taken over all strictly lower
triangular matrices $\widetilde{K}$, in order to preserve the explicit nature of the method.
A perturbation ${\widetilde{\bbA}}$ such that 
$${R}_{\hbox{\footnotesize Lin}}(K, \widetilde{\bbA}) = {R}^\textup{opt}_{\hbox{\footnotesize Lin}}(K) $$
will be called an {\it optimal perturbation of the method $K$ for the linear
problem}. 

Taking $\widetilde{\bbA}=0$ gives a (not perturbed) Runge--Kutta method \eqref{Runge--Kutta} and a (not perturbed) stability function $\phi_K$.
In this case we denote the threshold factor ${R}_{\hbox{\footnotesize Lin}}(K, 0)$ simply by $R(\phi_K)$.
Clearly
\begin{align}R(\phi_{K}) \le {R}^\textup{opt}_{\hbox{\footnotesize Lin}}(K)\, .\label{bound:RRL}
\end{align}
In the next section, we give upper bounds on ${R}^\textup{opt}_{\hbox{\footnotesize Lin}}(K)$. 


  
\subsection{Upper bounds on the threshold factor for optimal perturbations \label{sect:opt-dw-tf}}
For any explicit perturbed RK method $(K,\widetilde{\bbA})$ of  linear order $p$,
it can be seen that $\phi_{(K,\widetilde{\bbA})} \in \widetilde{\Pi}_{s,p}$, where
$\widetilde{\Pi}_{s,p}$, with $p\le s$, denotes the set of bivariate polynomials with
the following properties:
\begin{enumerate}
    \item $ \psi(z,\tilde{z}) = \displaystyle \sum_{j=0}^p \frac{z^j}{j!} 
            + \sum_{j=p+1}^s \sigma_j z^j + (z+\tilde{z})\Psi(z,\tilde{z}) $;
    \item $\Psi$ is a polynomial of combined degree at most $s-1$.
\end{enumerate}
In this section we investigate
$$\widetilde\cR_{s,p}=\sup\left\{R(\psi) \, |\, \psi(z, \tilde{z})\in\widetilde{\Pi}_{s,p}\right\}\, ,$$
Clearly
\begin{align}{R}^\textup{opt}_{\hbox{\footnotesize Lin}}(K)\leq \widetilde\cR_{s,p}\, . \label{bound:RLRsp}
\end{align}
However, not all functions in $\widetilde{\Pi}_{s,p}$ can
be realized as the stability function of an $s$-stage perturbed Runge-Kutta
method \eqref{PRunge--Kutta}, so inequality \eqref{bound:RLRsp} is often strict (see Example \ref{ex:2eo2} below).
In case the optimal polynomial is
realizable, the corresponding method may be
of interest for the integration of linear systems.

In Subsection \ref{sec:upperbound:lin} we give an upper bound for
$\widetilde\cR_{s,p}$.  In Subsection \ref{sec:upperlin:num}, we give an algorithm to
compute, $\widetilde\cR_{s,p}$ for given $s$ and $p$, along with numerical values.

\subsubsection{Upper bound on $\boldsymbol{\widetilde\cR_{s,p}}$}
\label{sec:upperbound:lin}~\\


\begin{lemma} \label{lem:1sp}
Let $\varphi(z)$ be a polynomial satisfying 
\begin{align}\label{polpert} 
\varphi(z) & = 1 + \gamma_1 z +\cdots + \gamma_p z^p +\gamma_{p+1} z^{p+1} + \cdots + \gamma_s z^s\\
\gamma_j & \geq \frac{1}{j!}, \quad  j=1, \ldots, p\,  . \nonumber
\end{align}
Then the radius of absolute monotonicity of  $\varphi$ satisfies
\begin{align}R(\varphi) \le\sqrt[p]{s(s-1)\ldots (s-p+1)}\, . \label{R1sp}
\end{align}
\end{lemma}

\begin{proof}  If  $R(\varphi)=0$, inequality \eqref{R1sp}
is trivial.  Let $\varphi(z)$ satisfy \eqref{polpert}
and be absolutely monotonic at $-r$ with $r>0$.
Then it can be written as 
$$\varphi(z)=\sum_{j=0}^s \alpha_j \left(1 + \frac{z}{r}\right)^j=\sum_{j=0}^s \alpha_j \left(\sum_{\ell=0}^j  \frac{z^\ell}{r^\ell} {j \choose \ell}\right)=\sum_{\ell=0}^s \left(\sum_{j=\ell}^s \alpha_j  {j \choose \ell}\right ) \frac{z^\ell}{r^\ell}\,, 
$$
where $\alpha_j\geq 0$, and $\sum_{j} \alpha_j=1$. As $\varphi$ is of the form \eqref{polpert},   the coefficient of $z^p$ is larger than $1/p!$. Some computations give 
$$\frac{1}{p !} \leq \left(\sum_{j=p}^s \alpha_j  {j \choose p}\right ) \frac{1}{r^p} \leq  \left(\sum_{j=p}^s \alpha_j  \right ) {s \choose p} \, \frac{1}{r^p} \leq  {s \choose p} \frac{1}{r^p}=\frac{s \, (s-1) \cdots (s-p+1)}{ p ! \, \, r^p}. 
$$
Consequently, 
$$r\leq \sqrt[p]{s(s-1)\cdots (s-p+1)}\, .  $$
\end{proof}

We remark that equality in \eqref{R1sp} is obtained for the polynomial
\begin{equation}\varphi(z)=\left(1+ \frac{z}{r}\right)^s\, , \label{optTihi}
\end{equation}
where $r=\sqrt[p]{s(s-1)\cdots (s-p+1)}$.

With the previous results, we can prove an upper bound on $\widetilde\cR_{s,p}$.

\begin{theorem}
\begin{align}\widetilde\cR_{s,p}\leq \sqrt[p]{s(s-1)\cdots (s-p+1)}\, .  \label{incRtsp}
\end{align}
\end{theorem}
\begin{proof} 
If $R(\psi)=0$ for all $\psi \in \widetilde{\Pi}_{s,p}$, then $\widetilde\cR_{s,p}=0$ and
inequality \eqref{incRtsp} is true. Otherwise, there exists a function  $\psi
\in \widetilde{\Pi}_{s,p}$  a.m. at $(-r,-r)$ with $r>0$.  
By \cite[Lemmas 2.9 and 2.10]{higueras2006}, $\psi$ is a.m.
at the points $(\xi, \xi)$, with $\xi \in [-r, 0]$.  
Writing $\psi(z, \tilde{z}) = \sum \sum \mu_{jk} z^j \tilde{z}^k$
and differentiating shows that all coefficients $\mu_{jk}$ are non-negative
since $\psi$ is a.m. at $(0,0)$.  Thus  $\psi(z,z)$ (viewed as a function
of one variable) is of the form \eqref{polpert} and is a.m. at $-r$.
Application of Lemma \ref{lem:1sp} gives the desired result.
\end{proof}


\subsubsection{Numerical computation of bounds $\boldsymbol{\widetilde\cR_{s,p}}$}\label{sec:upperlin:num}
In this section we provide a means to compute tighter bounds on $\widetilde\cR_{s,p}$ for 
given $s,p$, using linear programming.
The material in this section closely follows \cite[Section 4.6.2]{ketchesonphdthesis}.
The stability function \eqref{bistabfunc} of an explicit $s$-stage perturbed Runge--Kutta method 
\eqref{PRunge--Kutta} with linear order $p$ can also be written in the form
\begin{align} \label{dwtaylor1}
\psi(z,\tilde{z})= \sum_{j=0}^s \sum_{\ell=0}^j \gamma_{j\ell} 
  \left(1+\frac{z}{r}\right)^{j-\ell} \left(1+\frac{\tilde{z}}{r}\right)^{\ell}
  \quad \quad \mbox{with } \gamma_{j\ell} = \frac{r^j}{j!}\frac{\partial^j\psi_i}{\partial z^{j-\ell} \partial \tilde{z}^{\ell}};
\end{align}
furthermore, $\psi(z,-z) =\phi_K(z)= \exp(z)+ {\mathcal O}(z^{p+1})$.
After considerable manipulation we find that
$\psi(z,-z)=\sum_{i=0}^s C_i z^i$ where
\begin{align*}
C_i(r,\gamma) = \sum_{j=i}^s \sum_{\ell=0}^j \gamma_{j\ell} \sum_{m=\max(0,i-\ell)}^{\min(i,j-\ell)} 
{j-\ell \choose m}{\ell \choose i-m} \frac{(-1)^{i-m}}{r^i}
\end{align*}

Hence we have the following problem for existence of a polynomial
\eqref{dwtaylor1} with perturbed threshold factor at least $r$ and order at least $p$:
\begin{subequations}
\label{optprob3}
\begin{align}
\mbox{Given } r \mbox{ find $\gamma$ such that} & \nonumber \\
& \gamma_{j\ell} \ge 0 & 0 \le \ell\le j \le s \\
& C_i(r,\gamma) = \frac{1}{i!} & 0 \le i \le p.
\label{eqconstr}
\end{align}
\end{subequations}
Since \eqref{eqconstr} is a system of linear equations (in $\gamma$)
then for any given value of $r$ \eqref{optprob3} represents a linear programming
feasibility problem.  Hence we can use bisection and an LP solver to find 
$\widetilde\cR_{s,p}$, as was done for similar problems in \cite{ketcheson2008,ketcheson2009a}.
Table \ref{dwtable} gives computed values of $ \widetilde\cR_{s,p}$ for $s$ and $p$ up to ten.

\begin{table}
\begin{center}
\caption{$\widetilde\cR_{s,p}$: upper bounds on the threshold factors for optimal perturbations\label{dwtable}}
\begin{tabular}{|l|rrrrrrrrrr|}
\hline
s p  & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & 9 & 10 \\ \hline
1  & 1.00  & &  &  &  &  &  &  &  &  \\
2  & 2.00  & 1.41 &  &  &  &  &  &  &  &  \\
3  & 3.00  & 2.45 & 1.60 &  &  &  &  &  &  &  \\
4  & 4.00  & 3.46 & 2.49 & 2.00 &  &  &  &  &  &  \\
5  & 5.00  & 4.47 & 3.20 & 2.94 & 2.18 &  &  &  &  &  \\
6  & 6.00  & 5.48 & 4.00 & 3.65 & 3.11 & 2.58 &  &  &  &  \\
7  & 7.00  & 6.48 & 4.86 & 4.45 & 3.88 & 3.55 & 2.76 &  &  &  \\
8  & 8.00  & 7.48 & 5.77 & 5.31 & 4.57 & 4.32 & 3.72 & 3.15 &  &  \\
9  & 9.00  & 8.49 & 6.62 & 6.22 & 5.24 & 5.02 & 4.52 & 4.14 & 3.33 &  \\
10 & 10.00 & 9.49 & 7.42 & 7.09 & 5.95 & 5.70 & 5.25 & 4.96 & 4.32 & 3.73 \\ 
\hline
\end{tabular} 
\end{center}
\end{table}

\subsection{Examples}\label{sec:lin-examples}  
We now give some examples of optimal polynomials and Runge--Kutta methods.

\subsubsection{Polynomials achieving $\boldsymbol{\widetilde\cR_{s,p}}$}
The algorithm just described also provides coefficients for an optimal
polynomial, which may or may not be realizable as the stability function of
a perturbed Runge--Kutta method.

The optimal first-order polynomial for any $s$ is just the stability polynomial
of a (not perturbed) Runge--Kutta method consisting of $s$ repeated forward Euler steps.  

The optimal order-two polynomial of degree $s$ also has a simple form:
\begin{align} \label{optpoly2}
\psi_{s,2} (z, \tilde{z})& = \frac{2(s+r)-1}{2(s+r)}\left(1+\frac{z}{r}\right)^s + 
\frac{1}{2(s+r)}\left(1+\frac{\tilde{z}}{r}\right)^s,
\end{align}
where $r=\widetilde\cR_{s,2}=\sqrt{s(s-1)}$. Observe that bound \eqref{incRtsp} is sharp for $p=2$.  

Some of the other optimal polynomials also have rational coefficients.
Two optimal degree-four fourth order polynomials we found are 
\begin{align*}
\psi_{4,4}^{1} (z, \tilde{z})& = \frac{1}{3}\left(1+\frac{z}{r}\right)^2 + 
\frac{17}{48}\left(1+\frac{z}{r}\right)^4 +
\frac{14}{48}\left(1+\frac{z}{r}\right)^2\left(1+\frac{\tilde{z}}{r}\right)^2 +
\frac{1}{48}\left(1+\frac{\tilde{z}}{r}\right)^4\, , 
\end{align*}
and
$$ \psi_{4,4}^{2}(z, \tilde{z})= \frac{7}{16} \left(1+\frac{z}{r}\right)^4+\frac{3}{8}  \left(1+\frac{z}{r}\right)^2 \left(1+\frac{\tilde{z}}{r}\right)^2+\frac{1}{6} 
   \left(1+\frac{z}{r}\right)^3 \left(1+\frac{\tilde{z}}{r}\right)+\frac{1}{48} \left(1+\frac{\tilde{z}}{r}\right)^4\, , $$
where $r=\widetilde\cR_{4,4}=2$.
Thus the optimal polynomial in $\widetilde{\Pi}_{s,p}$ is in general not unique. 

We have also computed (not shown in this paper) the exact values of
$\widetilde\cR_{s,p}$ and polynomials $\psi_{s,p}$ for $p=3$ and $s=3,4,5,6$. In each 
case,  we found an optimal polynomial of the form \eqref{dwtaylor1}
where $p$ of the coefficients $\gamma_{sj}$ are non-zero and the coefficients
$\gamma_{ij}$ with $i<s$ are all zero.

Exact values of  $\widetilde\cR_{s,p}$ can be found in a systematic way as follows.
First, a high-precision approximation can be computed using bisection and
linear programming as described above.  
In practice, this yields a set of coefficients $\gamma_{jl}$ in which only $p$ values are non-zero  \cite{hadjimichael_draft}.
Setting
the remaining values to zero {\em a priori} in \eqref{eqconstr} yields a 
system of $p+1$ equations in $p$ unknowns which, nevertheless, possesses
a solution.  The solution may be found using a symbolic linear algebra package.
We do not pursue this further in the present work.



\begin{remark}
As noted already, not all polynomials of the form \eqref{dwtaylor1} can be realized as 
the stability function of a perturbed Runge--Kutta method \eqref{PRunge--Kutta} with $s$ stages.
For example, the polyomial \eqref{optpoly2} with $s=2$ is not the stability
function of any two-stage method (i.e. using only evaluations of 
$f(u_n), \tilde{f}(u_n), f(y_1), \tilde{f}(y_1)$).  
It can be realized as the stability function of a method that has three stages,
using evaluations of $f(u_n), \tilde{f}(u_n), f(y_1), \tilde{f}(y_2)$.  The difference 
in cost between such methods depends on the nature of $f, \tilde{f}$; see \cite{gottlieb2006}.
For this reason, we stress that the values
in Table \ref{dwtable} are only {\em upper bounds} on what can be achieved.
We do not pursue the topic further here.   \hfill $\square$
\end{remark}


\subsubsection{Optimal threshold factors for perturbations of specified 2-stage 2nd-order methods}
We have no general method for finding ${R}^\textup{opt}_{\hbox{\footnotesize Lin}}(K)$ nor a corresponding
method.  In this section we report results of some symbolic searches.
In the case of the second-order methods, due to the small number of free
parameters, it is not difficult to prove that the results below are truly
optimal.


\begin{example}\label{ex:2eo2}
We consider explicit perturbed second-order 2-stage methods
\begin{align}\begin{tabular}{c|cc}
$0$ & $0$ & $0$ \\
 $\alpha$ & $\alpha$ & $0$ \\
\hline 
$K$ \phantom{\huge I}  & $1-\frac{1}{2 \alpha}$ & $\frac{1}{2 \alpha}$ 
\end{tabular}
\qquad \begin{tabular}{c|cc}
  & $0$ & $0$ \\
  & $\tilde a_{21}$ & $0$ \\
\hline
$\widetilde{\bbA}$ \phantom{\huge I}  & $\tilde{b}_1$ & $\tilde{b}_2$ 
\end{tabular}\label{2eo2OptLin}.
\end{align}
For these methods, function \eqref{bistabfunc} can be expanded as
\begin{align}\phi_{(\bbA, \Apert)}(z,\tilde{z})= 1 + z + \frac{1}{2} z^2 + \beta_{11} z (z+\tilde{z}) + \beta_{1} (z+\tilde{z})+ \beta_2 (z+\tilde{z})^2\, , \label{2estord2}
\end{align}
where
\begin{align}\label{betas}
\beta_{11}= b^t \tilde{A} e + \tilde{b}^t A e= \tilde{b}_2 a_{21} + b_2 \tilde a_{21} \, , \qquad  \beta_1= \tilde{b}^t e= \tilde{b}_1+ \tilde{b}_2 \, , \qquad \beta_2=\tilde{b}^t \tilde{A} e=\tilde{b}_2 \tilde a_{21}
\end{align}
The polynomial \eqref{2estord2} is realizable (in the sense that it corresponds to a 2-stage Runge--Kutta method \eqref{2eo2OptLin}) if the first and last equations in \eqref{betas} can be solved for $\tilde a_{21}$ and $\tilde{b}_2$ in $\mathbb{R}$. A simple computation gives that a necessary condition is 
$\beta_{11}^2-2 \beta_2\geq 0$.  Note that the stability function is independent of $\alpha$.

With the help of Maple, we have computed the largest $r$ such that \eqref{2estord2} is a.m. at $(-r, -r)$ and the polynomial is realizable.  We have obtained that the optimal perturbation, denoted by $\widetilde{\bbA}_L$, satisfies $\tilde{b}_2=\tilde a_{21}=0$ and $\tilde{b}_1=\frac{1}{3} \left(\sqrt{7}-2\right)$, and \begin{equation}{R}^\textup{opt}_{\hbox{\footnotesize Lin}}(K)=\frac{1}{3} \left(1+\sqrt{7}\right)\approx 1.21525\, . \label{RLinopt}
\end{equation}
Observe that ${R}^\textup{opt}_{\hbox{\footnotesize Lin}}(K)< \widetilde\cR_{2,2}=\sqrt{2}$. 
The stability function \eqref{bistabfunc} for the  optimal perturbed method is 
$$\tilde{\phi}_{(K, \widetilde{\bbA}_L)}(z,\tilde{z})= \frac{1}{9} \left(4+\sqrt{7}\right) \left(1+ \frac{z}{r}\right)^2+ \frac{1}{9} \left(5-\sqrt{7}\right) \left(1+ \frac{\tilde{z}}{r}\right) \, , \label{2estord2opt}
$$
where $r={R}^\textup{opt}_{\hbox{\footnotesize Lin}}(K)$. 
\end{example}

\begin{example} \label{ex:4eo4} 
 We consider perturbations of the classical four stage, order four method, of the form
\begin{equation}
\begin{tabular}{c|cccc}
$0$ & $0$ & $0$ & $0$ & $0$ \\[1ex]
$\frac{1}{2}$ & $\frac{1}{2}$ & $0$ & $0$ & $0$ \\[1ex]
$\frac{1}{2}$ & $0$ &  $\frac{1}{2}$ &$0$ & $0$ \\[1ex]
$1$ & $0$ & $0$ & $1$ & $0$ \\[1ex]
\hline
$K$ \phantom{\huge I} & $\frac{1}{6}$  & $\frac{1}{3}$& $\frac{1}{3}$ & $\frac{1}{6}$   
\end{tabular} \qquad \begin{tabular}{c|cccc}
 & $0$ & $0$ & $0$ & $0$ \\[1ex]
 & $0$ & $0$ & $0$ & $0$ \\[1ex]
 & $\tilde a_{31}$ & $0$ & $0$ & $0$ \\[1ex]
& $\tilde a_{41}$ & $\tilde a_{42}$ & $0$ & $0$ \\[1ex]
\hline
$\widetilde{\bbA}$ \phantom{\huge I}  & $\tilde{b}_{1}$& $\tilde{b}_{2}$ & $0$ & $0$ 
\end{tabular}
\label{4eo4}
\end{equation}
We consider these perturbations because,
 in order to obtain a nonzero  SSP coefficient for nonlinear problems, the analysis done in \cite{higueras2005a} shows that only the entries $\tilde a_{31}$, $\tilde a_{41}$, $\tilde a_{42}$, $\tilde{b}_1$ and $\tilde{b}_2$ in $\widetilde{\bbA}$ need be nonzero. To study SSP coefficients for the linear case, we have to analyze the perturbed stability function, that in this case  is of the form 
\begin{align}\phi_{(\bbA, \Apert)}(z,\tilde{z}) = 1 + z + \frac{1}{2} z^2 +   \frac{1}{6} z^3+  \frac{1}{24} z^4+   \beta_{1} (z+\tilde{z})   + \beta_{11} z (z+\tilde{z})+ \beta_{21} z^2 (z+\tilde{z})\,  \label{esfun:4eo4}
\end{align}
where
$$\beta_{1}= \tilde{b}_1+\tilde{b}_2\, , \qquad \beta_{11}= \frac{1}{6} \left(3 \,\tilde{b}_2+2 \,\tilde a_{31}+\tilde a_{41}+\tilde a_{42}\right)\, , \qquad \beta_{21}=  \frac{1}{12} \left(2 \,\tilde a_{31}+\tilde a_{42}\right)  \, . 
$$
After some computations with \eqref{esfun:4eo4}, we obtain a coefficient
${R}_{\hbox{\footnotesize Lin}}(K,\widetilde{\bbA})\approx 1.66728$ that is the positive root of  the polynomial $15 \, x^4-4 \,x^3-12 \,x^2-24 \,x-24=0$.  The coefficients are 
$$\beta_{1}=\frac{ 7 \,r^3-2\, r^2-6 \,r-12 }{12} \, , \qquad  \beta_{11}=\frac{ 5\, r^2-2 \,r-6 }{12} \, , \qquad   \beta_{21}=\frac{r-1}{6}\, . 
$$
where $r={R}_{\hbox{\footnotesize Lin}}(K,\widetilde{\bbA})$. With these values, the perturbed stability function can be written as
$$\phi_{(\bbA, \Apert)}(z,\tilde{z}) =
  \gamma_{01} \left(1+\frac{\tilde{z}}{r}\right) +\gamma_{11} 
   \left(1+\frac{z}{r}\right)\left(1+\frac{\tilde{z}}{r}\right)+
  \gamma_{21}
  \left(1+\frac{z}{r}\right)^2 \left(1+\frac{\tilde{z}}{r}\right) + \gamma_{40} \left(1+\frac{z}{r}\right)^4\, , 
$$
where
$$\gamma_{01}=\frac{r \left(2\, r^3-r^2-6\right)}{6} \, ,    \quad \gamma_{11}=\frac{r^2 \left(r^2+2 \,r-6\right)}{12} \, , \quad \gamma_{21}=\frac{r^3\, (r-1) }{6}  \, ,  \quad \gamma_{40}= \frac{r^4}{24} \, . $$
 This perturbed stability function can be realized with the family of perturbations
\begin{subequations}\label{coef4eo4L}
\begin{align*}
\tilde a_{31}&= \frac{1}{2} \left(2 \, r-2-\tilde a_{42}\right)\, ,\, \\ \tilde a_{41}&=\frac{1}{2} \left(5  \,r^2-6 \,
   r-2-6  \,\tilde{b}_2\right)\,, \\ \tilde{b}_1&= \frac{1}{12} \left(7  \,r^3-2  \,r^2-6  \,r-12-12  \,\tilde{b}_2\right). 
\end{align*}
\end{subequations}
\end{example}

\section{Perturbed Runge--Kutta methods for nonlinear problems}\label{sec:nonlinear}

In this section we seek to answer the question posed in the
introduction: for a given method $K$, what perturbation $\widetilde{K}$ gives
the largest value of $R(K,\widetilde{K})$?  We begin by providing
some upper bounds.


It is convenient to write  
scheme \eqref{Runge--Kutta} in canonical Shu-Osher form \cite{SSPbook}
\begin{align} \label{cso}
    Y & = v_r u_n + \alpha_r \left( Y + \frac{h}{r} F\right)
\end{align}
where
\begin{align} \label{canonical-coefficients}
   v_r       = (I + rK)^{-1} e\, , \qquad 
     \alpha_r   = r(I + rK)^{-1} K\, .
\end{align}
Observe that matrices $K$ and $\alpha_r$ have the same structure (strictly lower triangular, lower triangular or full).

A Runge--Kutta method \eqref{Runge--Kutta} is said to be absolutely monotonic at $r$ if $(I + rK)^{-1}$
exists and all entries of $\alpha_r, v_r$ are non-negative. For a given method
$K$, the largest $r$ such that the method is absolutely monotonic at $r$ is
known as the SSP coefficient, Kraaijevanger coefficient, 
or radius of absolute
monotonicity \cite{Kraaijevanger1991}, and it is denoted herein by $R(K)$:
\begin{align}R(K)=\sup\left\{r\, | \, r=0 \text{ or } r>0, (I + rK )^{-1} \hbox{ exists, and }   \alpha_r,v_r\ge 0 \right\}\, . \label{ramdef0}
\end{align}
As usual, the inequalities above should be understood component--wise.

Next, we consider perturbed Runge--Kutta methods \eqref{PRunge--Kutta}. 
To study absolute monotonicity of perturbed Runge--Kutta
methods, we write method \eqref{PRunge--Kutta} also in a canonical
Shu-Osher-like form
\begin{align} \label{dwRunge--Kutta_cso}
Y & = \gamma_r u_n + \alpha^\textup{up}_r  \left(Y + \frac{h}{r}F\right) + 
                    \alpha^\textup{down}_r \left(Y - \frac{h}{r}\widetilde{F}\right)\, , 
\end{align}
where
\begin{subequations} \label{pert-cso-coeff}
\begin{align}
\gamma_r & = (I + rK + 2r\widetilde{\bbA})^{-1} e\,  ,\label{pert-cso-coeffa}\\[0.5ex]
\alpha^\textup{up}_r & = r(I + rK + 2r\widetilde{\bbA})^{-1} (K + \widetilde{\bbA})\, ,  \\[0.5ex] 
\alpha^\textup{down}_r & = r(I + rK + 2r\widetilde{\bbA})^{-1} \widetilde{\bbA}\, .
\end{align}
\end{subequations}
Observe that method  \eqref{dwRunge--Kutta_cso}, with $ \gamma_r=(I-\alpha^\textup{up}_r-\alpha^\textup{down}_r)e$, is a perturbed Runge--Kutta scheme with Butcher coefficients 
\begin{align}K=\frac{1}{r}(I-\alpha^\textup{up}_r-\alpha^\textup{down}_r)^{-1} (\alpha^\textup{up}_r-\alpha^\textup{down}_r)\, , \qquad \widetilde{\bbA}=\frac{1}{r}(I-\alpha^\textup{up}_r-\alpha^\textup{down}_r)^{-1} \alpha^\textup{down}_r\, , \label{AApert} \end{align}
provided that $(I-\alpha^\textup{up}_r-\alpha^\textup{down}_r)^{-1}$ exists. 

The {\em radius of absolute monotonicity} of a perturbed Runge--Kutta method  $(K,\widetilde{\bbA})$ is
the largest $r$ such that $\gamma_r$, $\alpha^\textup{up}_r$ and $\alpha^\textup{down}_r$ in \eqref{pert-cso-coeff} exist and are non-negative
\cite[Definition 3.1]{higueras2005a}:
\begin{align}R(K,\widetilde{\bbA})=\sup\left\{r\, | \, r=0 \text{ or } \, r>0, \,(I + rK + 2r\widetilde{\bbA})^{-1} \hbox{ exists, and }   \gamma_r,\, \alpha^\textup{up}_r,\, \alpha^\textup{down}_r\ge 0 \right\}\, . \label{ramdef}
\end{align}


\subsection{Zero-well-defined perturbations}
Regularity of $(I-\alpha^\textup{up}_r-\alpha^\textup{down}_r)$ is evidently important
in our study. Observe that from \eqref{pert-cso-coeff} we have
\begin{equation}\label{zerowell}
(I-\alpha^\textup{up}_r-\alpha^\textup{down}_r)(I+r K)=(I-2 \alpha^\textup{down}_r)\, . 
\end{equation}
Consequently, if $I+r K$ is regular for some $r$, then $(I-\alpha^\textup{up}_r-\alpha^\textup{down}_r)$ 
is regular if and only if $(I-2 \alpha^\textup{down}_r)$ is regular.

If $I-2\alpha^\textup{down}_r$ is singular, then the stage equations
do not have a unique solution even for the trivial ODE given
by $f=0$.  Hence we say that methods for which
$I-2\alpha^\textup{down}_r$ is singular are not {\em zero-well-defined}.  See \cite[Chap. 3]{SSPbook}
for the analogous definition in the context of traditional Runge--Kutta methods.


\subsection{Optimal perturbations}
The {\em optimal perturbed SSP coefficient} of a Runge--Kutta method $K$ is
denoted by
$$R^\textup{opt}(K)=\sup_{\widetilde{\bbA}} R(K,\widetilde{\bbA})\, . 
$$
For a given method $K$ that is (explicit/diagonally implicit/fully implicit), 
we consider the supremum  over perturbations $\widetilde{\bbA}$ that are zero-well-defined and
correspond to the same class of methods.
A matrix $\widetilde{\bbA}$ such that $R(K,\widetilde{\bbA})=R^\textup{opt}(K)$ is called an
optimal perturbation. For some methods, the optimal perturbation is not unique. 

The following result shows that every method can be perturbed so as to 
give a method with strictly positive SSP coefficient.
\begin{theorem}\label{thm:feasible}
Let $K$ be a Runge-Kutta method that belongs to a specified class of methods 
(explicit, diagonally implicit, or fully implicit).  Then it is always possible to find a perturbation $\widetilde{\bbA}$ within the same class such that $R(K,\widetilde{\bbA})>0$.
\end{theorem}
\begin{proof}From \cite[Proposition 3.7]{higueras2005a}, we have
$R(K,\widetilde{\bbA})>0$ if and only if the Butcher coefficients satisfy
\begin{align}\label{ineqa}
    K + \widetilde{\bbA} \geq 0 \, , \qquad \widetilde{\bbA} \geq 0\, ,
\end{align}
and the following inequalities  hold,
\begin{subequations} \label{ineqb}
\begin{eqnarray}
&& \hbox{Inc }((K+ 2 \widetilde{\bbA})  (K+ \widetilde{\bbA})) \leq 
\hbox{Inc }(K+ \widetilde{\bbA})\, , \label{3.6} \\[1ex]
&& \hbox{Inc }((K+ 2 \widetilde{\bbA}) \widetilde{\bbA} ) \leq 
\hbox{Inc }(\widetilde{\bbA})\, , 
\end{eqnarray}
\end{subequations}
where $\hbox{Inc } (F)$ denotes the incidence matrix of matrix $F$ defined as
$\hbox{Inc } (F)=(g_{ij})$ where $g_{ij}=1$ if $f_{ij}\neq 0$, and $g_{ij}=0$ 
if $f_{ij}=0$. 

Consider first the implicit case.  By making all entries of $\widetilde{\bbA}$ positive
we can satisfy \eqref{ineqb}, and by making them large enough we can satisfy
\eqref{ineqa}.  For the explicit and diagonally implicit cases, note that if $K, \widetilde{\bbA}$ are (strictly) lower-triangular, then the left-hand sides
of \eqref{ineqb} are also.  Thus by making all the (strictly) lower-triangular
entries of $\widetilde{\bbA}$ positive, and by taking them large enough, we can 
satisfy the above inequalities.
\end{proof}       
       
Observe that a perturbed Runge--Kutta method $(K,\widetilde{\bbA})$ can be interpreted as an
additive Runge--Kutta method $(K+\widetilde{\bbA},\widetilde{\bbA})$ for functions $(f, \tilde{f})$,
and conditions \eqref{pert-cso-coeff} are the ones required for the absolute
monotonicity of this additive scheme at $(z, \tilde{z})=(-r, -r)$ (see
\cite{higueras2006}). From Lemma 2.8 in \cite{higueras2006}, we obtain that
the stability function $\phi_{(\bbA, \Apert)}$ defined by \eqref{bistabfunc},
is  absolutely monotonic at $(\xi, \tilde \xi)=(-r, -r)$. Consequently,
\begin{equation} R(K,\widetilde{\bbA})\leq {R}_{\hbox{\footnotesize Lin}}(K,\widetilde{\bbA})\leq {R}^\textup{opt}_{\hbox{\footnotesize Lin}}(K) \, . \label{bound1}
\end{equation}
Furthermore, from  SSP theory and inequality \eqref{bound:RRL},  we have
\begin{equation}R(K)\leq R(\phi_K)\leq {R}^\textup{opt}_{\hbox{\footnotesize Lin}}(K)\, , \qquad  R(K)\leq  R^\textup{opt}(K)\leq {R}^\textup{opt}_{\hbox{\footnotesize Lin}}(K)\, . \label{bound2}
\end{equation}
The following example illustrates that $R(\phi_K)$ can be either larger
or smaller than $R^\textup{opt}(K)$.

\begin{example}\label{ex_2stage2order}
We consider the family of second order 2-stage Runge-Kutta methods \eqref{2eo2OptLin} for $\alpha\in\mathbb{R}$. For this family we have
\begin{align*}
v_r\geq 0 \quad & \Longleftrightarrow \qquad \alpha>0 \qquad \hbox{and}\quad 0\leq r\leq\frac{1}{\alpha}\, , \\[0.5ex]
\alpha_r\geq 0 \qquad & \Longleftrightarrow \qquad \alpha\geq\frac{1}{2} \quad \hbox{and}\quad 0\leq r\leq\frac{2\alpha-1}{\alpha}\, . 
\end{align*}
Thus
\begin{equation}\label{R2stageorder2}
R(K)=\begin{cases} 0, & \hbox{if} \quad \displaystyle  \alpha\leq \frac{1}{2}\, , \\[1ex]
\displaystyle  \frac{2 \, \alpha -1}{\alpha }\, ,  & \displaystyle \hbox{if}   \quad\frac{1}{2}<\alpha \leq 1\, ,\\[1ex]
\displaystyle   \frac{1}{\alpha }\, ,  & \hbox{if}   \quad 1<\alpha\, . 
\end{cases}
\end{equation}
In Figure \ref{f22} we show the threshold factor $R(\phi_K)$ (thin solid blue line) and the SSP coefficient  $R(K)$ (thick solid black line). We also show the corresponding optimal coefficients for perturbed methods, namely,  the optimal threshold factor ${R}^\textup{opt}_{\hbox{\footnotesize Lin}}(K)$  (thin dashed blue line) and the optimal SSP coefficient $R^\textup{opt}(K)$ for the perturbed method (thick dashed black line).

We see that for optimal SSP method ($\alpha=1$) it is not possible to increase the SSP coefficient by means of perturbations. However, for $\alpha=(\sqrt{7}-1)/2$ it is possible to obtain a perturbation that raises the SSP coefficient to $R^\textup{opt}(K)={R}^\textup{opt}_{\hbox{\footnotesize Lin}}(K)=\left(1+\sqrt{7}\right)\approx 1.21525$ (see \eqref{RLinopt}). 

 We also see that for  $2/3< \alpha<1$ we obtain 
$R(\phi_K) <R^\textup{opt}(K) 
$, whereas for $0<\alpha < 2/3$  and for $1< \alpha$  
we have $R^\textup{opt}(K)  < R(\phi_K) $.

Coefficients of the perturbations that give rise to these values are given in
Appendix \ref{sec:2stage-details}. 
\end{example}

\begin{figure}\begin{center}\includegraphics[scale=0.45]{Figures/2stage2order}
\caption{Family of second order 2-stage methods: SSP coefficients for unperturbed methods and optimal SSP coefficients for perturbed methods.
}
\label{f22}
\end{center}
\end{figure}


\subsection{Upper bounds on the SSP coefficient for perturbed RK methods\label{sec:nonlinear-bounds}}
In this section, we explore some upper bounds on the SSP coefficient $R^\textup{opt}(K)$ where $K$ is an $s$-stage order $p$ method. 
A straightforward upper bound is obtained from inequality \eqref{bound2} and Theorem \ref{incRtsp}:
\begin{align} R^\textup{opt}(K)\leq \sqrt[p]{s(s-1)\ldots (s-p+1)}\,  .\label{bound4}
\end{align}

As the next Theorem shows, the largest positive value such that vector $v_r$ in \eqref{canonical-coefficients} is non-negative is also an upper bound for $R^\textup{opt}(K)$.
\begin{theorem}\label{prop_vpos}
Consider an explicit Runge-Kutta method $K$  and let $r_e$ be the largest positive value such that vector $v_r$ in \eqref{canonical-coefficients} is non-negative. Then
\begin{equation}
R^\textup{opt}(K)\leq r_e\, . 
\end{equation}
\end{theorem}
\begin{proof}
Let $r=R^\textup{opt}(K)$.  Then $\gamma_r=(I-\alpha^\textup{up}-\alpha^\textup{down})e\geq 0$, and thus
from \eqref{canonical-coefficients}  and \eqref{zerowell}   we get
\begin{equation}(I-2\alpha^\textup{down}_r)v_r   \geq 0\, . \label{vpos}
\end{equation}
As $\alpha^\textup{down}_r\geq 0$, and since we consider only explicit, zero-well-defined perturbations,
$I-2\alpha^\textup{down}_r$ is an $M$ matrix.  Thus $(I-2\alpha^\textup{down}_r)^{-1}\geq 0$. If we multiply \eqref{vpos} by  $(I-2\alpha^\textup{down}_r)^{-1}$ we obtain that $v_r\geq 0$. 
\end{proof}

From Theorem \ref{prop_vpos} we obtain that
\begin{equation}
{R}(K)\leq R^\textup{opt}(K)\leq r_e\, . \label{in_re}
\end{equation}
Consequently, for those methods such that ${R}(K)=r_e$, the SSP coefficient cannot be increased by perturbation. This is the case for the family of second-order two-stage methods. 
For $\alpha\geq 1$, ${R}(K)= r_e=1/\alpha$ (see Example \ref{ex_2stage2order}). 

On the other hand, if ${R}(K)<r_e$ one can try to find a perturbation to increase the SSP coefficient. This is the case for the classical 4-stage order 4 method for which ${R}(K)=0$ and  $r_e\approx 1.2956$, the real root of $x^3 - 2 x^2 + 4 x - 4=0$. 


Another interesting bound, for explicit methods only, can be obtained in terms of the 
Butcher coefficients of the Runge-Kutta method $K$. 
\begin{theorem} \label{Pbound} Consider an explicit Runge-Kutta method $(K, \widetilde{\bbA})$ with perturbed SSP coefficient $R^\textup{opt}(K)>0$. Let $K=(a_{ij})$.  Then 
\begin{align}R^\textup{opt}(K)\leq  \frac{1}{\max_{ij} |a_{ij}|}\label{Rpert}
\end{align}
\end{theorem}
\begin{proof} The proof is similar to that of \cite[Lemma 3.2]{ruuth2004}.
Consider an optimal perturbation $\widetilde{\bbA}$ and set $r=R^\textup{opt} (K, \widetilde{\bbA})>0$;  consider too the canonical representation
\eqref{dwRunge--Kutta_cso}. Let
$\Lambda=\alpha^\textup{up}_r+\alpha^\textup{down}_r=(\alpha_{ij})$, $\Gamma=
\alpha^\textup{up}_r/r=(\beta_{ij})$, $\tilde \Gamma= \alpha^\textup{down}_r/r=(\tilde\beta_{ij})$;
observe that $\Lambda, \Gamma, \tilde \Gamma\geq 0$, and that $\Lambda=r
(\Gamma+\tilde \Gamma)$.  As $(I-\Lambda)e=\gamma_r\geq 0$ and $\alpha_{ik}\geq
0$, we have $\alpha_{ik}\leq 1$; as $(I-\Lambda)K=\Gamma-\tilde \Gamma$, we
have
\begin{align}a_{ik}= \beta_{ik}-\tilde \beta_{ik} + \sum_{j=k+1}^{i-1} \alpha_{ij} a_{jk}\, . \label{Acoef}
\end{align}
As $\alpha_{ik}= r (\beta_{ik} + \tilde \beta_{ik})$, then $\beta_{ik} + \tilde \beta_{ik}= \alpha_{ik}/r\leq 1/r$. In particular, from \eqref{Acoef}, 
$$|a_{21}|=\left|\beta_{21} - \tilde \beta_{21}\right|\leq \beta_{21} + \tilde \beta_{21}\leq \frac{1}{r}\, . $$
We proceed by induction on row $\ell$ of $K$. Assume that $|a_{ij}|\leq 1/r$, for $i= 2, \ldots, \ell$,  $j=1, \ldots , \ell-1$, and consider row $\ell+1$. Then, from  \eqref{Acoef}, 
\begin{align*}|a_{\ell+1, 1}|&=\left|\beta_{\ell+1, 1}-\tilde \beta_{\ell+1, 1} + \sum_{j=2}^{\ell} \alpha_{\ell+1, j} \, a_{j, 1}\right|\leq \beta_{\ell+1, 1}+\tilde \beta_{\ell+1, 1} + \sum_{j=2}^{\ell} \alpha_{\ell+1, j} |\, a_{j, 1}|\\
&\leq  \frac{1}{r} \alpha_{\ell+1, 1}+ \frac{1}{r} \sum_{j=2}^{\ell} \alpha_{\ell+1, j} \leq  \frac{1}{r} \sum_{j=1}^{\ell} \alpha_{\ell+1, j}  \leq \frac{1}{r}  \, . 
\end{align*}
A similar argument can be used to show that $|a_{\ell+1,j}|\leq 1/r$, $j= 2, \ldots, \ell$. The Theorem follows by induction. 
\end{proof}


Consequently, $$R(K)\leq  R^\textup{opt} (K)\leq \frac{1}{\max_{ij} |a_{ij}|}\,  .$$ 
For those methods such that
$R (K)= 1/{\max_{ij} |a_{ij}|}
$
it is not possible to increase the SSP coefficient by perturbing the method.
This is the case for all known optimal explicit SSP RK methods of orders one through
four, with any number of stages \cite{ketcheson2008}.


For the restricted class of perturbations considered in \cite{ruuth2004},
similar results were obtained in  \cite[Theorems 3.1, 3.4, 3.5 and
3.6]{ruuth2004}. Theorem \ref{Pbound} extends those results, showing that
no improvement in the radius of absolute monotonicity is possible for many 
optimal SSP methods, even when more general perturbations are considered.


\subsection{Relations among the Butcher and canonical Shu-Osher representations}


For perturbations of explicit Runge--Kutta methods (or of other methods with
one stage equal to $u_n$) there exists a certain simple transformation that may yield a larger
value of $R(K,\widetilde{\bbA})$ -- and that never yields a smaller value.
\begin{proposition} \label{propfirstcol}
	Let an $s$-stage explicit perturbed Runge--Kutta method be given with coefficients 
    $\gamma_r, \alpha^\textup{up}_r, \alpha^\textup{down}_r \ge 0,$ where $r$ is the radius of absolute
    monotonicity of the method.  Consider the perturbed method with coefficients
    \begin{subequations} \label{transform}
    \begin{align}
    	\widehat{\gamma} & = (1, 0 , \ldots, 0)^t\, ,  \\
        \widehat{\alpha}^\textup{up}_{i,1} & = \alpha^\textup{up}_{i,1} + (\gamma_r)_{i}/2 & 2 \le i \le s  \\
        \widehat{\alpha}^\textup{down}_{i,1} & = \alpha^\textup{down}_{i,1} + (\gamma_r)_{i}/2 & 2 \le i \le s .
    \end{align}
    \end{subequations}
    Then the perturbed method with coefficients $(\gamma_r, \alpha^\textup{up}_r, \alpha^\textup{down}_r)$
    and the modified perturbed method with coefficients $(\widehat{\gamma}, 
    \widehat{\alpha}^\textup{up}, \widehat{\alpha}^\textup{down})$
    correspond to the same RK method $K$.  The modified perturbed method has radius of absolute
    monotonicity at least equal to $r$.
\end{proposition}
\begin{proof}
It is easily seen that the modified method is equivalent to the original
one when $f=\tilde{f}$, so they correspond to the same unperturbed method.
Meanwhile, the transformation never leads to negative coefficients,
so the modified method is a.m. at $r$.
\end{proof}
\begin{remark}Proposition \ref{propfirstcol} is also valid for Runge-Kutta methods whose first row is equal to zero.
\end{remark}

In the Butcher form \eqref{PRunge--Kutta} it is obvious which
perturbed methods $(K,\widetilde{K})$ correspond to a given method $K$.
In the canonical Shu-Osher form it is less obvious.
The following lemma characterizes which methods of the form \eqref{dwRunge--Kutta_cso}
are perturbations of a given method \eqref{cso}.
\begin{lemma}\label{equivprop}
If method \eqref{dwRunge--Kutta_cso} is a perturbation of method \eqref{cso}, then 
their coefficients are related as follows:
\begin{subequations} \label{alpharel1}
\begin{align} 
 (I-2\alpha^\textup{down}_r)\alpha_r & = (\alpha^\textup{up}_r - \alpha^\textup{down}_r) \label{alpharel1a} \\
 (I-2\alpha^\textup{down}_r)v_r & = \gamma_r.
\end{align}
\end{subequations}
Furthermore, if \eqref{alpharel1} holds and 
the perturbation is zero-well-defined,
then \eqref{dwRunge--Kutta_cso} is a perturbation of \eqref{cso}.
\end{lemma}
\begin{proof}
To prove the first part,
take $\tilde{f} = f$ in \eqref{dwRunge--Kutta_cso} to obtain:
\begin{align*}
Y & = \gamma_r u_n + (\alpha^\textup{up}_r + \alpha^\textup{down}_r) Y
                    + (\alpha^\textup{up}_r - \alpha^\textup{down}_r) \frac{h}{r}F.
\end{align*}
Subtract $2\alpha^\textup{down}_r Y$ from both sides to get
\begin{align} \label{zwd}
(I-2\alpha^\textup{down}_r) Y & = \gamma_r u_n + (\alpha^\textup{up}_r - \alpha^\textup{down}_r) 
                    \left(Y +  \frac{h}{r}F\right).
\end{align}
Substituting \eqref{cso} in the above gives
\begin{align*}
(I-2\alpha^\textup{down}_r) v_r u_n + (I-2\alpha^\textup{down}_r)\alpha_r \left(Y +  \frac{h}{r}F\right)
    & = \gamma_r u_n + (\alpha^\textup{up}_r - \alpha^\textup{down}_r) \left(Y +  \frac{h}{r}F\right),
\end{align*}
Equating coefficients yields \eqref{alpharel1}.

To prove the second part, assume $I-2\alpha^\textup{down}_r$ is invertible and
write \eqref{alpharel1} as
\begin{subequations} \label{alpharel2}
\begin{align}
 \alpha_r & = (I-2\alpha^\textup{down}_r)^{-1}(\alpha^\textup{up}_r - \alpha^\textup{down}_r) \\
 v_r & = (I-2\alpha^\textup{down}_r)^{-1}\gamma_r.
\end{align}
\end{subequations}
Substitute \eqref{alpharel2} in \eqref{cso}, multiply on the left by
$(I-2\alpha^\textup{down})^{-1}$, and follow the steps above in reverse.
\end{proof}

 
Lemma \ref{equivprop} does {\em not} imply that the perturbation \eqref{dwRunge--Kutta_cso}
is unique for a given $r$; see Proposition \ref{propfirstcol}. 


\begin{remark} \label{invertibility-remark}
The necessity of the zero-well-defined condition in the second part of
Proposition \ref{equivprop} can be seen from the following example.
We take the implicit trapezoidal Runge-Kutta method
\begin{align*}
 \begin{array}{c|cc}
0 & 0 & 0 \\
1/2 & 1/2 & 1/2 \\ \hline
    & 1/2 & 1/2
\end{array}.
\end{align*}
The canonical form \eqref{cso} is then
\begin{align*}
\alpha_r & = \begin{pmatrix}
0 & 0 & 0 \\[1ex]
\frac{r}{r+2} & \frac{r}{r+2} & 0 \\[1ex]
\frac{r}{r+2} & \frac{r}{r+2} & 0
\end{pmatrix}\, ,
& v_r & = \begin{pmatrix} 1 \\[1ex] \frac{2-r}{r+2} \\[1ex] \frac{2-r}{r+2} \end{pmatrix}.
\end{align*}
Then \eqref{alpharel1} is satisfied -- for any $r$ -- by
\begin{align*}
\alpha^\textup{up} = \alpha^\textup{down} & = \begin{pmatrix} 
1/3 & 0 & 0 \\
0 & 1/2 & 0 \\
0 & 1/2 & 0
\end{pmatrix}\, , &
\gamma & = \begin{pmatrix} 1/3 \\ 0 \\ 0 \end{pmatrix}.
\end{align*}
However, this method -- which involves a perturbation that is not zero-well-defined 
-- is not a perturbation of the original method. \hfill $\square$
\end{remark}



\subsection{Computing optimal perturbations\label{sec:algorithms}}
In this section we present two algorithms to symbolically
or numerically find the optimal perturbed SSP coefficient and a corresponding perturbation of a given RK method.  The first algorithm is proven to approximate the
optimal value to any accuracy, contingent on the computational solution of
linear program subproblems.  It is only valid for explicit perturbations.  The second
algorithm is analytical and exact, and valid for both explicit and implicit
methods, but it is not proven to give the optimal value.  The results of the
two algorithms coincide (to high precision) for all explicit methods on which 
we have tested them.

\subsubsection{Provably optimal algorithm for explicit perturbations}
In the foregoing, we have shown that finding an optimal perturbation
consists of determining the largest $r$ such that there exists a splitting
satisfying \eqref{alpharel1} with positive coefficients.
Note that the range of
values for which a method $(K,\widetilde{K})$ is absolutely monotonic is
always the interval $[0, R(K,\widetilde{K})]$.  Therefore, one way to
find the largest $r$ is to devise a method for testing for
a given $r$ whether there exists a perturbation $\widetilde{K}$ such that $R(K,\widetilde{K})\ge r$.
For given method \eqref{Runge--Kutta} and value of $r$,
the system of equations \eqref{alpharel1} together with the inequalities
$\alpha^\textup{up}_r, \alpha^\textup{down}_r \ge 0$
constitutes a linear programming (LP) feasibility problem.
The following theorem is an immediate consequence of Lemma \ref{equivprop}.

\begin{theorem}
Let an $s$-stage RK method $K$ and a positive number $r$ be given.  There exists a 
perturbation $\widetilde{\bbA}$ with $R(K,\widetilde{\bbA})\ge r$ if and only if there exists
an $(s+1)\times(s+1)$ matrix $\alpha^\textup{down}_r$ such that $(I-2\alpha^\textup{down}_r)$ is regular and
the following componentwise inequalities hold:
\begin{subequations} \label{LP}
\begin{align} 
 (I-2\alpha^{down}_r)\alpha_r + \alpha^{down}_r & \ge 0 \\
 (I-2\alpha^{down}_r)v_r & \ge 0 \\
 \alpha^\textup{down}_r \ge 0.
\end{align}
\end{subequations}
\end{theorem}
The linear program \eqref{LP} can be solved by standard LP solvers.
By embedding this solution in a one-dimensional root-finding algorithm,
optimal perturbations can be found.  An algorithm based on bisection follows. 
\begin{algorithm}\caption{Optimal explicit perturbation}
\label{alg:explicit-LP}
\begin{algorithmic}
\INPUT $K$
\State $r_\textup{max} := 1/\max|a_{ij}|, r_\textup{min} := 0$.
\While {$r_\textup{max} - r_\textup{min} > \epsilon$}
	\State $r=\frac{r_\textup{max} + r_\textup{min}}{2}$.
	\State Compute the coefficient matrices $\alpha_r, v_r$ using \eqref{canonical-coefficients}.
	\State Solve the LP given by \eqref{LP}.
    \If {it is feasible} 
    	\State $r_\textup{min}:=r$  
    \Else
    	\State $r_\textup{max}:=r$
    \EndIf
\EndWhile \\
\Return $r_\textup{min}$
\end{algorithmic}
\end{algorithm}

Assuming the solution of the LP is correct, the algorithm provably
finds an optimal explicit perturbation.  However,
for implicit perturbations the LP solver may converge to a solution
(like the method in Remark \ref{invertibility-remark} above)
for which $I-2\alpha^\textup{down}_r$ is singular.

\subsubsection{Iterated splitting algorithm}
We next investigate how to choose $\alpha^\textup{up}_r, \alpha^\textup{down}_r$ directly so as
to find a perturbation with radius of a.m. at least $r$.
The following result suggests an approach.
\begin{lemma} \label{signlemma}
Given an explicit Runge--Kutta method \eqref{cso},
let $\alpha^\textup{up}_r\ge0, \alpha^\textup{down}_r\ge0$ denote coefficients of a 
zero-well-defined perturbation of \eqref{cso}.
Then there exist matrices $\alpha^+\ge 0$, $\alpha^-\ge 0$ such that
\begin{subequations} \label{alphar}
\begin{align}
    \alpha^\textup{up}_r     & = (I+2\alpha^-)^{-1} \alpha^+ \\
    \alpha^\textup{down}_r   & = (I+2\alpha^-)^{-1} \alpha^-,
\end{align}
\end{subequations}
and $\alpha_r=\alpha^+ - \alpha^-$.
\end{lemma}
\begin{proof}
Since the perturbation is zero-well-defined, we can define
\begin{subequations} \label{a1a2}
\begin{align}
    \alpha^+ & = (I-2\alpha^\textup{down})^{-1} \alpha^\textup{up}_r \\
    \alpha^- & = (I-2\alpha^\textup{down})^{-1} \alpha^\textup{down}_r.
\end{align}
\end{subequations}
Then, by \eqref{alpharel1a}, $\alpha_r=\alpha^+-\alpha^-$.  Furthermore, since $I-2\alpha^\textup{down}_r$ is an $M$-matrix,
we have $\alpha^+\ge0$ and $\alpha^-\ge0$.
Solving \eqref{a1a2} for $\alpha^\textup{up}_r,\alpha^\textup{down}_r$ gives \eqref{alphar}.
\end{proof}

In the next algorithm we use the following notation:
\begin{align*}
((x)^+)_{ij} & = \begin{cases}
x_{ij} & \mbox{ if } x_{ij}\ge0 \\
0 & \mbox{ if } x_{ij}<0.  \end{cases} &
((x)^-)_{ij} & = \begin{cases}
0 & \mbox{ if } x_{ij}\ge0 \\
-x_{ij} & \mbox{ if } x_{ij}<0\, ,   \end{cases}
\end{align*}
and thus $x=(x)^+-(x)^-$ is a sign splitting of matrix $x$, with $(x)^+\geq 0$, $(x)^-\geq 0$. 

Given a perturbed Runge-Kutta method \eqref{dwRunge--Kutta_cso} with $\gamma_r=e_1$, where $e_1=(1, 0, \ldots, 0)^t$, and $\alpha^\textup{up}$ or $\alpha^\textup{down}$ containing negative values,  
 we construct 
\begin{subequations}\label{newpert}
\begin{align}
\tilde \gamma_r&=\left(I+2 \, (\alpha^\textup{up}_r)^-+ 2 \,(\alpha^\textup{down}_r)^-\right)^{-1}  e_1\, ,\\[0.5ex]
\tilde\alpha_r^\textup{up}&=\left(I+2 \,(\alpha^\textup{up}_r)^-+ 2 \,(\alpha^\textup{down}_r)^-\right)^{-1} \left((\alpha^\textup{up}_r)^++ (\alpha^\textup{down}_r)^-\right)\, ,   \\[0.5ex]
\tilde\alpha_r^\textup{down}&=\left(I+2 \,(\alpha^\textup{up}_r)^-+ 2 \,(\alpha^\textup{down}_r)^-\right)^{-1}  \left((\alpha^\textup{up}_r)^-+(\alpha^\textup{down}_r)^+\right)\, ,  
\end{align}
\end{subequations}
where 
$\alpha^\textup{up}=(\alpha^\textup{up}_r)^+-(\alpha^\textup{up}_r)^-$, $\alpha^\textup{down}_r=(\alpha^\textup{down}_r)^+-(\alpha^\textup{down}_r)^-$, provided that $I+2 (\alpha^\textup{up}_r)^-+ 2 (\alpha^\textup{down}_r)^-$ exists. Using Lemma \ref{equivprop}, it is straightforward to prove that, if method $\alpha^\textup{up}_r$, $\alpha^\textup{down}_r$ is a perturbation of   \eqref{cso}, then \eqref{newpert} is also perturbation of \eqref{cso}. 
Next, for explicit methods, we perform transformation \eqref{transform}. 
In this way, \eqref{newpert} followed by transformation \eqref{transform} give a perturbation of the form \eqref{dwRunge--Kutta_cso} with $\gamma_r=e_1$,  that we denote by  $\hat \alpha_r^\textup{up}$, $\hat \alpha_r^\textup{down}$. If $\hat \alpha_r^\textup{up}\geq 0$,  $\hat \alpha_r^\textup{down}\geq 0$, then $r$ is an SSP coefficient; otherwise, we can repeat the above process.  

The following lemma studies the sign of $(\hat \alpha_r^\textup{up})_{ij}$, $(\hat \alpha_r^\textup{down})_{ij}$ when   $(\alpha^\textup{up}_r)_{ij}<0$ or $(\alpha^\textup{down}_r)_{ij}<0$. For the sake of clarity, we drop the index $r$. 

\begin{lemma}\label{lemma_row} We consider a perturbed explicit Runge-Kutta method  with coefficients $\gamma=e_1$, $\alpha^\textup{up}$, $\alpha^\textup{down}$, and 
the perturbation   $\hat \alpha^\textup{up}$, $\hat \alpha^\textup{down}$ obtained by computing \eqref{newpert} followed by transformation \eqref{transform}.  Assume that  $j_0\geq 2$ is the first row with negative terms in $\alpha^\textup{up}$ or $\alpha^\textup{down}$. Let $m_0$ be the largest index $m_0\geq 1$ such that  $\alpha^\textup{up}_{j_0,m_0}<0$ or $\alpha^\textup{down}_{j_0,m_0} <0$. Then
\begin{enumerate}
\item For first to $(j_0-1)$-th row, we have: $\hat \alpha^\textup{up}_{i,j}=\alpha^\textup{up}_{i,j}$ and $\hat \alpha^\textup{down}_{i,j}=0$ for $1\leq i\leq j_0-1$, $1\leq j\leq j_0-2$. 
\item For the $j_0$-th row, we have:
\begin{enumerate}
\item If  $m_0=1$, then $\hat \alpha^\textup{up}_{j_0,1}<0$ or $\hat \alpha^\textup{down}_{j_0,1}<0$. 
\item If $m_0\geq 2$, then,  $\hat \alpha^\textup{up}_{j_0,m_0}\geq 0$  and $\hat \alpha^\textup{down}_{j_0,m_0}\geq 0$. 
\item For $1\leq m_0\leq j_0-2$, we have $\hat \alpha^\textup{up}_{j_0,\ell}\geq 0$  and $\hat \alpha^\textup{down}_{j_0,\ell}\geq 0$ for $\ell=m_0+1, \ldots, j_0-1$. 
\end{enumerate}
\end{enumerate}
 
 \end{lemma}

 \begin{proof}  If $j_0$ is the first row with negative terms in $\alpha^\textup{up}$ or $\alpha^\textup{down}$, straightforward computations give that $\hat \alpha^\textup{up}_{i,j}=\alpha^\textup{up}_{i,j}$ and $\hat \alpha^\textup{down}_{i,j}=0$ for $1\leq i\leq j_0-1$, $1\leq j\leq j_0-2$, and
\begin{subequations} \begin{align}
\hat \alpha^\textup{up}_{j_0,1}&=  \alpha^\textup{up}_{j_0,1} -2  \sum_{i=2}^{j_0-1} \left((\alpha^\textup{up}_{j_0,i})^-+ (\alpha^\textup{down}_{j_0,i})^-\right) \,  \alpha^\textup{up}_{i,1}\, , \label{d1} \\\quad
\hat \alpha^\textup{up}_{j_0,\ell}&=  (\alpha^\textup{up}_{j_0,\ell})^++ (\alpha^\textup{down}_{j_0,\ell})^- -2  \sum_{i=\ell+1}^{j_0-1} \left((\alpha^\textup{up}_{j_0,i})^-+ (\alpha^\textup{down}_{j_0,i})^-\right) \, \alpha^\textup{up}_{i,\ell} \, , \quad \ell=2, \ldots, j_0-1\, .\label{dm} 
\end{align}
\end{subequations}
and
\begin{subequations}
\begin{align}
\hat \alpha^\textup{down}_{j_0,1}&=  \alpha^\textup{down}_{j_0,1} -2  \sum_{i=2}^{j_0-1} \left((\alpha^\textup{up}_{j_0,i})^-+ (\alpha^\textup{down}_{j_0,i})^-\right) \, \alpha^\textup{down}_{i,1}\, , \label{d11} \\\quad
\hat \alpha^\textup{down}_{j_0,\ell}&=  (\alpha^\textup{up}_{j_0,\ell})^-+ (\alpha^\textup{down}_{j_0,\ell})^+  -2  \sum_{i=\ell+1}^{j_0-1} \left( (\alpha^\textup{up}_{j_0,i})^-+ (\alpha^\textup{down}_{j_0,i})^-\right) \, \alpha^\textup{down}_{i,\ell} \, , \quad \ell=2, \ldots, j_0-1\, .\label{dm1} 
\end{align}
\end{subequations}
Let $m_0$ be the largest index $m_0\geq 1$ such that  $\alpha^\textup{up}_{j_0,m_0}<0$ or $\alpha^\textup{down}_{j_0,m_0} <0$. In this case,
$\alpha^\textup{up}_{j_0,i}\geq 0$, $\alpha^\textup{down}_{j_0, i}\geq 0$ for $i=m_0+1, \ldots, j_0-1$, and  thus 
$$(\alpha^\textup{up}_{j_0,i})^+=\alpha^\textup{up}_{j_0,i}\, , \quad (\alpha^\textup{down}_{j_0,i})^+=\alpha^\textup{down}_{j_0,i}\, , \quad (\alpha^\textup{up}_{j_0,i})^-=(\alpha^\textup{down}_{j_0,i})^-=0\, ,\qquad i=m_0+1, \ldots, j_0-1\, . $$ 
 If $m_0=1$, from \eqref{d1} and \eqref{d11} we get  $\hat \alpha^\textup{up}_{j_0,1} = \alpha^\textup{up}_{j_0,1}$ and $\hat \alpha^\textup{down}_{j_0,1}= \alpha^\textup{down}_{j_0,1}$,  and thus $\hat \alpha^\textup{up}_{j_0,1}<0$ or $\hat \alpha^\textup{down}_{j_0,1}<0$.  If $m_0\ge  2$, from \eqref{dm} and \eqref{dm1} we get
 $$\hat \alpha^\textup{up}_{j_0,m_0} = (\alpha^\textup{up}_{j_0,m_0})^++ (\alpha^\textup{down}_{j_0,m_0})^- \geq 0\, , \qquad \hat \alpha^\textup{down}_{j_0,m_0} =  (\alpha^\textup{up}_{j_0,m_0})^-+(\alpha^\textup{down}_{j_0,m_0})^+ \geq 0\, . 
 $$
Finally, for $1\leq m_0\leq j_0-2$, from  \eqref{dm} and \eqref{dm1} we get that, for $\ell=m_0+1, \ldots, j_0-1$, we have
$$\hat \alpha^\textup{up}_{j_0,\ell} = (\alpha^\textup{up}_{j_0,\ell})^+ \geq 0\, ,  \qquad \hat \alpha^\textup{down}_{j_0,\ell} =  (\alpha^\textup{down}_{j_0,\ell})^+ \geq 0\,.   
$$
 \end{proof}
 
Consequently, if matrices $\alpha^\textup{up}_r$ and $\alpha^\textup{down}_r$  contain negative elements in the second or later columns, an iterated construction of perturbations $\hat \alpha_r^\textup{up}$, $\hat \alpha_r^\textup{down}$ removes these negative values obtaining a perturbation with non-negative elements from second column on. However, if in a row $j_0$ we have:
\begin{equation}\label{stop}
  \alpha^\textup{up}_{j_0,1}<0 \, \qquad \hbox{and} \qquad \alpha^\textup{up}_{j_0,\ell}\geq 0 \qquad \ell=2, \ldots, j_0-1\, ,
\end{equation}
or
\begin{equation}\label{stop2}
  \alpha^\textup{down}_{j_0,1}<0 \, \qquad \hbox{and} \qquad \alpha^\textup{down}_{j_0,\ell}\geq 0 \qquad \ell=2, \ldots, j_0-1\, ,
\end{equation}
 the new perturbation $\hat \alpha_r^\textup{up}$, $\hat \alpha_r^\textup{down}$ will also contain negative elements in the first column.    


We now give Algorithm \ref{alg:2} to determine whether there exists a perturbation 
with a.m. radius $r$ for a given method.
\begin{algorithm}\caption{Existence of a perturbation with radius $r$}
\label{alg:2}
\begin{algorithmic}
\INPUT $r, K$
    \State Compute the coefficient matrices $\alpha_r, v_r$ using \eqref{canonical-coefficients}.

    \State Set $\alpha^\textup{up} = \alpha_r$ and $\alpha^\textup{down} = 0$.
	\While {$\alpha^\textup{up}$ or $\alpha^\textup{down}$ has any negative entries}
	    \State If $K$ has a zero row, perform the transformation \eqref{transform}.
	    \State If $\alpha^\textup{up}, \alpha^\textup{down} \ge 0$, stop.  This is a feasible perturbation.
        \State If condition \eqref{stop} or \eqref{stop2} hold, stop. A feasible perturbation cannot be found.   
        \State Set $\alpha^- = (\alpha^\textup{up})^- + (\alpha^\textup{down})^+$ and 
    		  $\alpha^+ = (\alpha^\textup{up})^+ + (\alpha^\textup{down})^-$
        \State Compute a new splitting:
    		\begin{align*}
            	\alpha^\textup{up} & = \left(I + 2 ((\alpha^\textup{up})^- + (\alpha^\textup{down})^-)\right)^{-1}\alpha^+ \\
                \alpha^\textup{down} & =\left(I + 2((\alpha^\textup{up})^- + (\alpha^\textup{down})^-)\right)^{-1}\alpha^-
            \end{align*}
\EndWhile \\
\end{algorithmic}
\end{algorithm}

\begin{remark}
This approach seems to lead to optimal splittings for all the explicit
methods on which we have tested it.  However, for all implicit methods
we have tested, it fails to increase the radius of absolute monotonicity at all.
Even for explicit methods, we have no proof that it's optimal because one could use
$(\alpha_r)^+ + \delta, (\alpha_r)^- + \delta,$
in place of $(\alpha_r)^+, (\alpha_r)^-$, where $\delta$ is any non-negative matrix.  \hfill $\square$
\end{remark}


\subsection{Examples}\label{sec:examples}
In this section we compute optimal perturbations of some existing methods,
using the algorithms described in the last section.

\begin{table}
\begin{tabular}{ll|l|llllll}
\toprule
Order & Stages & Method & $R(K)$ & $R^\textup{opt}(K)$ & Bound \eqref{Rpert} & Bound \eqref{bound4} & Property C \\
\midrule
1 & 1 & Forward Euler & 1 &          1 &      1 &       1 &                  True \\
\midrule
2 & 2 & Midpoint         & 0 &      0.732 &      1 &   1.414 &                   True \\
  & 2 & Minimal trunc. error & 0.5 &          1 &  1.333 &   1.414 &                   True \\
  & 2 & SSP22 \cite{shu1988}          & 1 &          1 &      1 &   1.414 & True \\
  & 2 & SSP22* \cite{gottlieb2006} & 0.784 &      1.215 &  1.215 &   1.414 & True \\
\midrule
3 & 3 & Heun33 \cite{heun1900neue} & 0 &      0.776 &  1.333 &   1.817 &  False \\
  & 3 & SSP33 \cite{shu1988}   & 1 &          1 &      1 &   1.817 & True \\
\midrule
4 & 4 & RK44 (Kutta)          & 0 &      0.685 &      1 &   2.213 &                  False \\
  & 5 & Merson \cite{merson1957operational} & 0 &      0.242 &    0.5 &   3.309 & False \\
  & 10 &SSP104 \cite{ketcheson2008}       & 6 &          6 &      6 &   8.425 & False \\
\midrule
5 & 6 & Fehlberg \cite{fehlberg1969klassische}     & 0 & 0.057 &  0.125 & 3.727 & False \\
  & 7 & Dormand-Prince \cite{Dormand1980} & 0 &      0.040 &  0.086 &   4.789 & False \\
  & 8 & Bogacki  \cite{Bogacki1996} & 0 &      0.313 &  0.859 &   5.827 & False \\
  & 7 & SSP75 \cite{ruuth2004} & 0 &      1.396 &  1.792 &   4.789 & False \\
  & 8 & SSP85 \cite{ruuth2004} & 0 &      1.875 &  1.919 &   5.827 & True \\
  & 9 & SSP95 \cite{ruuth2004} & 0 &      2.738 &  3.198 &   6.853 & False \\
\midrule
6 & 9 & Calvo \cite{calvo1990} & 0 & 0.021 &  0.059 &   6.265 & False \\
\midrule
8 & 13 & Prince-Dormand \cite{Prince1981} & 0 &      0.013 &  0.059 &   9.212 & False \\
\bottomrule
\end{tabular}
\caption{Properties of some RK methods and their optimal perturbations.  The optimal perturbed
radius of absolute monotonicity was computed by both the linear programming algorithm
and the iterated splitting algorithm; in every case they gave identical results (up to 
roundoff errors).  Decimal values have been truncated to the number of digits shown.
\label{tbl:explicit-methods}}
\end{table}


We have computed optimal perturbations for several known explicit methods using the two 
algorithms described above.  In all cases, the two algorithms gave the same
values.  It thus seems possible that Algorithm 2
also gives truly optimal results in general, but we do not have a proof.  Properties
of the methods studied are given in Table \ref{tbl:explicit-methods}.  Several interesting
facts are evident:
\begin{itemize}
	\item For all optimal SSP methods (up to order four), perturbation cannot yield a 
          larger coefficient.  This is evident already from the bound \eqref{Rpert}.
          For all other methods, some improvement is achieved.
	\item Consistent with Theorem \ref{Pbound}, for every method considered, 
    	  it is possible to achieve $R^\textup{opt}>0$ by some perturbation.          
    \item The simple bound \eqref{Rpert} predicts the optimal coefficient to within a factor of three in
    		every case.
    \item The methods SSP75, SSP85, and SSP95 are optimal methods found in \cite{ruuth2004},
    		with property C.  By considering methods without property C, we obtain slightly
            larger coefficients for perturbations of SSP75 and SSP95.  On the other hand,
            relaxing the column assumption gives no benefit in the case of the SSP85 method.
    \item The values found have been truncated to three decimal places but are known to greater precistion.  For the 4-stage, order-four method of Kutta, the three-digit value of $R^\textup{opt}(K)$ given in the table matches the value found by Shu and Osher.  However, the exact (irrational) value is slightly larger and is given in the appendix.
\end{itemize}


\section{Conclusions}\label{sec:conclusions}
In this work we have studied SSP coefficients for perturbations of a given explicit Runge-Kutta method. We have considered both the linear and the nonlinear case, and have obtained useful bounds on the threshold factor and on the radius of absolute monotonicity
for perturbed Runge--Kutta methods.  
 We have also provided an algorithm for computing
optimal perturbations of explicit Runge-Kutta methods, and given optimal perturbations for many methods from the literature.

This work seems to provide a complete picture for the case of most interest: explicit
methods applied to nonlinear problems.
Nevertheless, some other interesting issues remain unsolved.  These include:
\begin{itemize}
	\item A method to compute optimal perturbations for linear problems.
	\item An algorithm for obtaining optimal splittings of implicit methods.
\end{itemize}
These may be a starting point for future work.


\section{Appendix}

In this section we give additional details on SSP coefficients and optimal perturbations of
second order 2-stage Runge--Kutta methods and the classical 4-stage fourth order Runge--Kutta method. 

\subsection{Second order 2-stage methods} \label{sec:2stage-details}

We consider the family of 2-stage second order methods \eqref{2eo2OptLin}. In example \ref{ex:2eo2}   we studied perturbations that increase the SSP coefficient for the linear case. For nonlinear problems, in example \ref{ex_2stage2order},  figure \ref{f22} shows the values of $R^\textup{opt}(K)$ for $\alpha\in[-3,3]$. 

In this section,  for each $\alpha$,  we give the expressions for $R^\textup{opt}(K)$ and we show optimal perturbations $\widetilde{\bbA}_{NL}$ such that $R(K, \widetilde{\bbA}_{NL})=R^\textup{opt}(K)$.  It is important to point out the convenience of choosing 
$\widetilde{\bbA}_{NL}= \widetilde{\bbA}_{L}$, where  $\widetilde{\bbA}_{L}$ denotes the optimal perturbation for the linear case. In this case, we have not only $R(K, \widetilde{\bbA}_{L})=R^\textup{opt}(K)$ but also ${R}_{\hbox{\footnotesize Lin}}(K,  \widetilde{\bbA}_{NL})={R}^\textup{opt}_{\hbox{\footnotesize Lin}}(K)$. The computations required to obtain the results in this section have been done with the symbolic computation program \emph{Mathematica}.

If we denote by $r=R^\textup{opt}(K)$,  we have that
\begin{equation}
r= \begin{cases}
\displaystyle \frac{1}{|\alpha|}\, ,  &  \qquad \displaystyle  \hbox{if} \quad  \alpha \in\left(-\infty, -\frac{1}{2} \left(1+\sqrt{7}\right)\right]\bigcup  \left[\frac{1}{2} \left(-1+\sqrt{7}\right), \infty\right)\,,  \\[3ex]
\displaystyle   \frac{-1+\alpha+\sqrt{ 3 \alpha ^2-2 \alpha +1} 
   }{|\alpha| }\, ,  &   \qquad \displaystyle  \hbox{if} \quad  \alpha\in\left(-\frac{1}{2} \left(1+\sqrt{7}\right), 0\right) \bigcup \left(0, \frac{1}{2} \left(-1+\sqrt{7}\right)\right)\, . \end{cases}\label{radiopert22}
\end{equation}
Next we give optimal perturbations $\widetilde{\bbA}_{NL}$. 

For $\alpha<0$, we obtain that it is not possible to obtain a perturbation of the form \eqref{2eo2OptLin}  with $\tilde b_2=0$ and $\tilde a_{21}=0$. Consequently, $\widetilde{\bbA}_{NL}\neq \widetilde{\bbA}_{L}$ and we always have that ${R}_{\hbox{\footnotesize Lin}}(K,  \widetilde{\bbA}_{NL})<{R}^\textup{opt}_{\hbox{\footnotesize Lin}}(K)$. Optimal perturbations of the form \eqref{2eo2OptLin} for different values of $\alpha<0$ must satisfy the following conditions.
\begin{itemize}   
\item For $-\frac{1}{2} \left( 1+\sqrt{7}\right)\leq\alpha <0$, the coefficients $\tilde a_{21}$, $\tilde b_1$ and $\tilde b_2$ in $\widetilde{\bbA}_{NL}$ must satisfy
   $$-\alpha \leq \tilde a_{21}\leq \frac{1-r \, \alpha}{2 \, r}\, , \qquad \tilde b_1=-\frac{r\, \tilde a_{21}}{2\, \alpha } \, , \qquad \tilde b_2=-\frac{1}{2\, \alpha}\, , 
   $$
   where $r=R^\textup{opt}(K)$.
\item For $\alpha\leq-\frac{1}{2} \left( 1+\sqrt{7}\right)$,    
   we should have 
   $$\tilde a_{21}=-\alpha,\qquad  -\frac{1}{2 \,\alpha }\leq\tilde b_1\leq \frac{-2 \,\alpha ^2-2
 \,  \alpha +1}{4\, \alpha }\, , \qquad  -\frac{1}{2 \,\alpha }\leq
  \tilde b_2\leq \frac{2\, \alpha  \,\tilde b_1-1}{4\, \alpha
   }\, . $$
\end{itemize}
For  $\alpha>0$ we can find optimal perturbations with $\tilde b_2=0$ and $\tilde a_{21}=0$. Coefficient $\tilde{b}_1$ must satisfy the following conditions.
\begin{itemize}
\item For $0<\alpha \leq  
   \left(-1+\sqrt{7}\right)/2$,  we have that   
   \begin{align}  \qquad \tilde{b}_1=\frac{\sqrt{3 \alpha ^2-2
   \alpha +1}-\alpha}{2 \alpha }  \, . \label{tb1NL}
   \end{align}   
Thus there is a unique $\widetilde{\bbA}_{NL}$ of the form \eqref{2eo2OptLin}. In this case, we   have  $R(K)<R(K, \widetilde{\bbA}_{NL})=R^\textup{opt}(K)$.

   \item For $\left(-1+\sqrt{7}\right)/2<\alpha< 1$, we also get $R(K)<R^\textup{opt}(K)$,  but in this case the optimal perturbation  $\widetilde{\bbA}_{NL}$ is not unique. All the perturbations with $\tilde{b}_1$ satisfying
$$   \frac{1-\alpha }{\alpha } \leq \tilde{b}_1\leq
   \frac{2 \alpha ^2-2 \alpha +1}{4 \alpha } \, , 
$$
are optimal. In particular, we can take $\widetilde{\bbA}_{NL}=\widetilde{\bbA}_{L}$. With this choice, $R(K, \widetilde{\bbA}_{L})=R^\textup{opt}(K)=1/\alpha$ and ${R}_{\hbox{\footnotesize Lin}}(K,  \widetilde{\bbA}_{L})={R}^\textup{opt}_{\hbox{\footnotesize Lin}}(K)\approx 1.22$. Furthermore,   $\alpha=\left(-1+\sqrt{7}\right)/2$  provides the largest  SSP coefficient within the family of 2-stage second order method (see figure \ref{f22}).

\item For $1\leq \alpha$, we have $R(K)=R^\textup{opt}(K)=1/\alpha$ and the optimal perturbation $\widetilde{\bbA}_{NL}$  is not unique. All the values
$$0\leq \tilde{b}_1\leq
   \frac{2 \alpha ^2-2 \alpha +1}{4 \alpha } 
$$
give  optimal perturbations. We can take $\widetilde{\bbA}_{NL}=0$,   but in this case ${R}_{\hbox{\footnotesize Lin}}(K, 0)< {R}^\textup{opt}_{\hbox{\footnotesize Lin}}(K)$. A better choice is $\widetilde{\bbA}_{NL}=\widetilde{\bbA}_{L}$. Observe that, for $\alpha=1$, we get the optimal SSP coefficient $R(K)=1$ that cannot be increased by perturbations.
  \end{itemize}

Next, we consider some concrete values of $\alpha$ to show the the expressions of the perturbations. For each value, we give the Butcher tableau of the perturbation and matrices  $\alpha^\textup{up}$ and $\alpha^\textup{down}$ in \eqref{dwRunge--Kutta_cso}. 
\begin{itemize}
\item For $\alpha=1/2$ we get method RK2a in \cite{HuKoVa} with $R(K)=0$. With perturbation
\begin{align*}
\widetilde{\bbA}= \begin{pmatrix}
0 & 0 & 0 \\
0 & 0 & 0\\
\tilde{b}_1 & 0 & 0
\end{pmatrix} , \, \alpha^\textup{up}= \begin{pmatrix}
0 & 0 & 0 \\
\tilde{b}_1 & 0 & 0\\
0 & 2 \tilde{b}_1  & 0
\end{pmatrix}
  , \,  \alpha^\textup{down} = \begin{pmatrix}
0 & 0 & 0 \\
0 & 0 & 0\\
1-2\tilde{b}_1  & 0 & 0
\end{pmatrix}
 , \, \gamma = \begin{pmatrix}
1\\
1-\tilde{b}_1\\
0
\end{pmatrix}, 
\end{align*}
where $\tilde{b}_1=\frac{1}{2} \left(\sqrt{3}-1\right)$, we get  $R(K,\widetilde{\bbA}) ={R}_{\hbox{\footnotesize Lin}}(K, \widetilde{\bbA})= \sqrt{3}-1$. 
\item For $\alpha=2/3$,  we have a nontrivial SSP coefficient $R^\textup{opt}(K)=1/2$, but we can increase this value to $R(K,\widetilde{\bbA}_1)=R^\textup{opt}(K)=1$ with perturbation
\begin{align*}
\widetilde{\bbA}_1= \begin{pmatrix}
0 & 0 & 0 \\
0 & 0 & 0\\
1/4 & 0 & 0
\end{pmatrix}
\, , \, \alpha^\textup{up}= \begin{pmatrix}
0 & 0 & 0 \\
2/3 & 0 & 0\\
0 & 3/4 & 0
\end{pmatrix}
\, , \,  \alpha^\textup{down} = \begin{pmatrix}
0 & 0 & 0 \\
0 & 0 & 0\\
1/4 & 0 & 0
\end{pmatrix}
\, , \, \gamma  = \begin{pmatrix}
1\\
1/3\\
0
\end{pmatrix}\, . 
\end{align*}
For this perturbation, $R(\phi_K)={R}_{\hbox{\footnotesize Lin}}(K, \widetilde{\bbA}_1)=1$. We can take $\gamma=(1,0,0)^t$ by modifying the first column of $\alpha^\textup{up}$ and $\alpha^\textup{down}$ according to \eqref{transform},
\begin{align*}
\widetilde{\bbA}_2= \begin{pmatrix}
0 & 0 & 0 \\
1/6 & 0 & 0\\
3/8 & 0 & 0
\end{pmatrix}\, , \, \alpha^\textup{up}= \begin{pmatrix}
0 & 0 & 0 \\
5/6 & 0 & 0\\
0 & 3/4 & 0
\end{pmatrix}
\, , \,  \alpha^\textup{down} = \begin{pmatrix}
0 & 0 & 0 \\
1/6 & 0 & 0\\
1/4 & 0 & 0
\end{pmatrix}
\, , \, \gamma = \begin{pmatrix}
1\\
0\\
0
\end{pmatrix}. 
\end{align*}

\item As it has been pointed out above, the largest value in the $\alpha$-family of 2-stage second order schemes is $R^\textup{opt}(K)=(1+\sqrt{7})/3$ and it is obtained for $\alpha=(\sqrt{7}-1)/2$. The perturbation is  of the form \eqref{2eo2OptLin} with $\tilde{b}_1= \left(\sqrt{7}-2\right)/2$, and
\begin{align*}
\alpha^\textup{up}= \begin{pmatrix}
0 & 0 & 0 \\
1  & 0 & 0\\
0 & \frac{1}{9} \left(4+\sqrt{7}\right)  & 0
\end{pmatrix}
\, , \,  \alpha^\textup{down} = \begin{pmatrix}
0 & 0 & 0 \\
0 & 0 & 0\\
\frac{1}{9} \left(5-\sqrt{7}\right)  & 0 & 0
\end{pmatrix}
\, , \, \gamma_r = \begin{pmatrix}
1\\
0 \\
0
\end{pmatrix}
\end{align*}
This is the perturbation obtained in \cite[Table V]{gottlieb2006} by numerical search in the class of perturbations considered in \cite{gottlieb2006}. 
\end{itemize}


\subsection{Classical fourth order 4-stage method}\label{4eo4L}


For nonlinear problems, applying the analysis above, we find that the optimal 
perturbation of the classical method has SSP coefficient given by
the real root of $x^3+2 x^2+4 x-4=0$, which is approximately
$R^\textup{opt}(K)\approx 0.685016$.  The corresponding perturbation is
{\em not unique.}  For instance, we can take $\gamma_r=(1,0,0,0,0)$,
and all entries of $\alpha^\textup{down}_r$ equal to zero except
\begin{align}
	(\alpha^\textup{down}_r)_{31} & = \frac{r^2}{4}   & (\alpha^\textup{down}_r)_{42} & = \frac{r^2}{2},  
\end{align}
where $r=R^\textup{opt}(K)$.
However, there exist other optimal perturbations with additionally $(\alpha^\textup{down}_r)_{42}=\epsilon$ where $0\le \epsilon \le 0.782$.  

We remark that nearly-optimal perturbations for this method are given
in \cite[p. 448]{shu1988} and \cite{Higueras2010}.
Interestingly, these different perturbed methods have different values of
${R}_{\hbox{\footnotesize Lin}}(K, \widetilde{\bbA})$.



