\section{Introduction}

Microarray gene expression experiments help in the measurement of expression levels of thousands of genes simultaneously. Such data help in diagnosing various types of tumors with better accuracy. The fact that this process generates a lot of complex data happens to be its major limitation. Normally the number of genes (features) is much greater than the number of samples (instances) in a microarray gene expression dataset. Such structures pose problems to machine learning and make the problem of classification difficult to solve. This is mainly because, out of thousands of genes, most of the genes do not contribute to the classification process. As a result gene subset selection acquires extreme importance towards the construction of efficient classifiers with high predictive accuracy.

To overcome this problem, one way is to select a small subset of informative genes from the data. This technique which is known as \textit{gene selection} or \textit{feature selection} helps in tackling overfitting by getting rid of noisy genes,  reducing the computational load and in increasing the overall classification performance of the learning models.

Gene selection algorithms are mainly categorized as : wrappers and filters. Wrappers make use of learning algorithms to estimate the quality or suitability of genes to the modelling problem. Optimization algorithms in combination with various classifiers fall into this category as described in \cite{vkj-aco, aco-rf, vkj-book}. On the other hand, filters  \cite{john-feature} evaluate the genes considering their inherent characteristics without making use of a learning algorithm. Filters, therefore, give an insight into the properties of the dataset we use. Algorithms based on statistical tests and mutual information are some examples of filters.

This paper presents hybrid BBO \textendash \space RF and hybrid BBO \textendash \space SVM approaches for simultaneous informative gene selection and high performance classification. Additionally, for enhancing performance we provide information gain gene ranking as heuristic knowledge to our BBO algorithm. It traverses the enormously large search space by using this ranking information to iteratively obtain informative gene subsets. The selected subsets of genes (candidate solutions) in each generation are subsequently evaluated by SVM and Random Forests CV (cross validation) accuracies.

\section{Methodology}
\subsection{Biogeography\textendash based Optimization}
Biogeography is the study of geographical distribution of species over geological period of time. Biological literature on the same is massive. In 2008, for the first time, Simon \cite{simonbbo} applied the biogeography analogy to the idea of engineering optimization and thus introduced the Biogeography\textendash based Optimization (BBO) technique. It is a  population based method that works with a collection of candidate solutions over generations. It attempts to explore the  combinatorially large solution spaces with a stochastic approach like many other evolutionary algorithms \cite{pso-eberhart, ea-de}. It mimics the geographical distribution of species to represent the problem and its candidate solutions in the search space, subsequently using the process of species migration and mutation to redistribute solution instances across the search space in quest of globally optimal or near optimal solutions.

BBO, as is or in variations, has been explored for various combinatorial and constrained/unconstrained optimization problems \cite{bboblended} including the likes of the Traveling Salesman Problem \cite{bbotspMoXu,bbotspSong}, satellite image classification \cite{panchal} and sensor selection \cite{simonbbo} among others. But as of 2012, no work is reported of using BBO as a gene selection technique for microarray gene expression data analysis. We attempt to study BBO for gene selection and classification in this work.

In BBO, there exists an ecosystem (population) which in turn consists of a number of habitats (islands). Each habitat has a \textit{habitat suitability index} (HSI), which is similar to a fitness function and depends on many features/attributes  of the island. If a value is assigned to each feature, then the HSI of a habitat \textit{H} is a function of these values. These variables characterizing a habitat's suitability collectively form the \textit{\textquoteleft suitability index variables\textquoteright} \space (SIVs). Thus,
\begin{displaymath}
 HSI (Habitat_{i}) \space \rightarrow \space f(SIV_{1}, SIV_{2}, \ldots , SIV_{m})
\end{displaymath}

For the problem of gene selection, the SIVs of a habitat (candidate solution) are the selected subsets of genes out of the
set of all genes. The ecosystem is therefore a random collection of candidate gene subsets.

A good solution is thus analogous to a good HSI and \textit{vice versa}. Good HSI solutions tend to share SIVs with poor HSI solutions. This form of sharing, termed as migration, is controlled by emigration and immigration rates of the habitats. We have purposely kept the model simple and have obeyed the original simple linear model for migration as shown in Figure \ref{fig_sim}.

\begin{figure}[h]
\centering
\includegraphics[width=2.0in]{Screenshot-1.pdf}
\caption{Migration rate vs. No. of species}
\label{fig_sim}
\end{figure}

where E and I are the maximum emigration and immigration rates, both typically set to 1. The individual immigration and emigration rates ($\lambda$ and $\mu$ respectively) are calculated by the same formulae for this simple linear model as in \cite{simonbbo}.

\begin{displaymath}
 \lambda_{k} = I\left( 1 - \frac{k}{n}\right)
\end{displaymath}
\begin{displaymath}
 \mu_{k} = \frac{Ek}{n}
\end{displaymath}

where $k$ is the iterator for the $n$ habitats.

\subsection{BBO Gene Selection Algorithm}
We present our BBO algorithm for performing gene selection.
For our problem, we treat a gene (identified by gene number) as an SIV for a habitat and each habitat has $m$ SIVs (arity $m$ as \textit{Habitat}  $H \in SIV^{m}$). For example, if a habitat is $Habitat_{1} = \{12, 345, 26, 7, 141\}$ then the SIVs 12, 345, \ldots, 141 are the selected gene numbers out of say a collection of 500 genes and the subset size is 5 genes. The tuned parameter values for our algorithm are given in the next section. The BBO gene selection algorithm is stated.

\begin{algorithm}
\caption{BBO for Gene Selection}
\label{bbo-gene}
\begin{algorithmic}[1]
 \STATE Initialize BBO parameters
 \STATE Initialize $ecosystem$ with randomly generated $n$ habitats
 \STATE Evaluate habitats : calculate the HSI of each habitat in the $ecosystem$ \textit{(Cross validation accuracies of each gene subset from a classifier)}
 \FOR{$G$ generations}
    \STATE Compute $\lambda_{i}$ and $\mu_{i}$ for each $Habitat_{i}$ based on its HSI
    \STATE Perform $migration$
    \STATE Perform $mutation$
    \STATE Re-evaluate $ecosystem$
    \STATE Perform $elitism$
 \ENDFOR
 \STATE Output the habitat with best HSI and its SIVs (selected genes).
\end{algorithmic}
\end{algorithm}

At each migration or mutation, we ensure that a gene is not duplicated within a single subset of genes, \textit{i.e.} within a single habitat. The subset sizes (the variable $m$) can be set during each run of the BBO algorithm. The ecosystem then has habitats all with same number of SIVs. This selection of subset size is predecided and can be tuned manually by running BBO for various subset sizes.

$Migration$: The migration procedure of the original BBO algorithm \cite{simonbbo} is retained. We produce the algorithm here for the purpose of completeness.
\begin{algorithm}
\caption{Migration}
\label{bbo-migration}
\begin{algorithmic}[1]
 \STATE Select $H_{i}$ with probability $\propto$ $\lambda_{i}$
 \IF{$H_{i}$ is selected}
    \FOR{$j=0$ \TO $n$}
    \STATE Select $H_{j}$ with probability $\propto$ $\mu_{i}$
    \IF{ $H_{j}$ is selected}
      \STATE Randomly select an SIV $\sigma$ from $H_{j}$
      \STATE Replace a random SIV in $H_{i}$ with $\sigma$
    \ENDIF
    \ENDFOR
 \ENDIF
\end{algorithmic}
\end{algorithm}

\textit{Mutation with Information Gain based Gene Ranking}: Given the vast search space formed by the possible genes, we keep the mutation rate to about 0.4 to 0.55 in order to graze other portions of this space with a good chance. We have used information gain heuristics as the additional information during mutation. The \textit{information gain} (IG) \cite{info-gain} of a gene is a measure of attribute selection. It stores the \textquoteleft information content\textquoteright \space of a gene with respect to the problem under consideration. The IG of a gene indicates the capability to separate instances for binary classification. We are specially interested in the non-zero IG values. Thus we partition the informative and non-informative genes into separate sets. For IG computation, we have used Weka \cite{weka} data mining software suite, which outputs the information gain based ranking of genes. This is fed to BBO for further computations. This informative gene set with non-zero infogain values are introduced in the population during the process of mutation.

While in mutation, our algorithm either randomly explores newer genes or exploits from the available set of genes with non-zero infogain values. We set a user defined exploitation probability $q_{0}$. This exploitation is done in a probabilistic manner as analogous with the exploration and exploitation in ant colony optimization \cite{aco-dorigo, aco-specialsection, aco-blum, vkj-aco, vkj-aco2}. To give an example, we have a total of  30 expressed genes and only the first 8 out of these 30 genes have a non-zero infogain. Assuming rand \textless \space $q_{0}$ is satisfied (step 5 of Algorithm 3), then we  select one out of these 8 genes, with a probability proportional to the information gain ranking scores, to be newly put in $Habitat_{i}$ in place of an existing one.  While on the other hand if rand \textless \space $q_{0}$ did not get satisfied, we execute the \textbf{\textit{else}} part (step 7 and 8) $i.e.$ randomly select a gene from all of the 30 genes expressed in the data.

Algorithm 3 gives the detailed mutation algorithm.

\begin{algorithm}[h]
\caption{Mutation}
\label{bbo-mutation}
\begin{algorithmic}[1]
 \FOR{$j=0$ \TO $m$}
    \STATE Use $\lambda_{i}$, $\mu_{i}$ of habitat $H_{i}$ to compute the probability $P_{i}$
    \STATE Select SIV(gene) $H_{i}(j)$ with probability $\propto$ $P_{i}$
    \IF{ $H_{i}(j)$ is selected}
      \IF{ rand \textless \space $q_{0}$  }
	  \STATE \textit{Exploit} : Replace $H_{i}(j)$ with a probabilistically selected a SIV (gene) from the rest (using their information gain)
      \ELSE
	  \STATE \textit{Explore} : Replace $H_{i}(j)$ with a random SIV (gene) out of the rest
      \ENDIF
    \ENDIF
    \ENDFOR
\end{algorithmic}
\end{algorithm}

$Elitism$: We implement elitism so that the best solutions obtained until a particular generation do not get corrupted.

\subsection{Support Vector Machines}
Support Vector Machines (SVMs) \cite{svm} were originally introduced by Vapnik and co-workers \cite{boser} and successively extended by a number of other researchers. SVM employs a maximum margin linear hyperplane for solving binary linear classification problems. For non-linearly separable problems, SVM first transforms the data into a higher dimensional feature and subsequently employs a linear hyperplane. To deal with computational intractability issues it further uses appropriate kernel functions facilitating all computations in the input space itself. Vapnik et al. in \cite{geneVapnik} have themselves used SVM with recursive feature elimination (RFE) for gene selection and achieved notably high accuracy levels. We discuss more about results in the subsequent section.

For our purposes we employ the libSVM \cite{libsvm} library for evaluation of our candidate solutions during each generation.

\subsection{Random Forests}

Random Forests (RF) were first introduced by Breimen and Cutler \cite{rfbreiman}. It is an ensemble of randomly constructed independent decision trees. It performs substantially better than single-tree classifiers such as CART \cite{CART} and C4.5 \cite{c45}. A random subset of attributes are used for node splitting while growing each decision tree. Normally, for each tree, a bootstrap set (with replacement) is drawn from the original training data, $i.e.$ an instance is picked from the training data and is replaced again before drawing the next instance. Likewise, $n$ such instances are taken to form \textquoteleft in bag\textquoteright \space set for a particular tree. For each of the bootstrap training sets, about one \textendash \space third of the samples, on an average, are unused for making the \textquoteleft in bag\textquoteright \space data and are called the \textquoteleft out of bag\textquoteright \space (OOB) data for that particular tree.
The classification tree is built with this \textquoteleft in bag\textquoteright \space data using the CART algorithm \cite{CART}. Separate test data is not required in RF for checking the overall accuracy of the forest. The OOB data is used for cross validation. When all the trees are grown, the $k^{th}$ tree classifies the samples that are OOB
for that tree (left out by the $k^{th}$ tree). In this manner, each instance is classified by about one third of the trees. A majority vote is then taken to decide on the class label for each case. The percentage of times that the voted class label is not equal to the original class of a sample, averaged over all the cases in the training data, is called as the OOB error rate \cite{bagging}.

We have used the randomForest package in $R$ for implementation purposes \cite{RFmanual}.

\section{Discussion and Results}
\subsection{Datasets}
The output of microarray experiments are the expression levels of different genes. Three such datasets were obtained from the Kent Ridge Biomedical datasets repository\cite{datasets} and libSVM repository \cite{libsvm} (made available from various other original sources).

The Colon Cancer dataset retrieved from the Kent Ridge Biomedical dataset repository consists of 62 instances representing cell samples taken from colon cancer patients. Among these, 40 are tumor samples while 22 otherwise \cite{colon-alon}. The breast cancer dataset is retrieved from the DUKE Breast Cancer SPORE frozen tissue bank \cite{breast-west}. Of the 44 samples we worked with, each sample with expressions for 7129 genes, 22 belong to class A (estrogen receptor-positive ER$+$) while 22 belong to class B (estrogen receptor-negative ER$-$). The Leukemia dataset \cite{leuk-golub} also retrieved from the Kent Ridge Biomedical dataset repository contains the expression of 7129 genes. These are total 72 samples taken from leukemia patients out of which 25 belong to the Acute Myeloid Leukemia (AML) class and 47 belong to the Acute Lymphoblastic Leukemia (ALL) class. These specifications are tabulated in Table I.

\begin{table}[h]
\caption{Dataset Specifications}
\label{dataset_info}
\centering
\begin{tabular}[h]{lccr}

\hline
 Cancer dataset &  \#genes &  \#classes& \#instances\\
 name (D)  &          &           &  (\#A \& \#B)\\
\hline
 Colon (C)  &  2000&  2&  62 (40 \& 22)\\
 Breast (B) &  7129&  2&  44 (22 \& 22)\\
 Leukemia (L)&  7129&  2& 72 (25 \& 47)\\
\hline
\end{tabular}
\end{table}

\subsection{Discussion and Results}

As discussed earlier while describing our algorithms, we have implemented BBO with and without heuristics. Very interestingly, both simple BBO and BBO with heuristics are successful in selecting a good set of features providing comparable classification results. As compared to very typical implementations of gene selection using genetic algorithms and other EAs like PSO \cite{pso-eberhart, ea-de}, which usually work over many generations with large population sizes \cite{compare-evo1}, our implementation of both versions of BBO (that with SVM and RF) for gene selection started showing better results very early with just 40-50 habitats in the ecosystem. We performed 50 simulations each with \#generations varying from 15 to 40. The algorithms almost always converged to comparable results by the end of 25 generations with very minute differences in further generations.

\begin{figure}[h]
\centering
\includegraphics[width=2.8in, height=2.5in]{graph.png}
\caption{Example convergence of population averages of simple vs. heuristic BBO over 25 generations; also with example non-monotonic behavior (boxed) of average suitability of habitats with simple BBO}
\label{fig_graph}
\end{figure}

For BBO with heuristics, we decided not to hinder the original BBO procedure a lot and incorporated the heuristics of information gain of each gene only during the mutation process; even this with only some degree and not for every mutation. This was done by roping in the analogy of exploration and exploitation from the ant colony optimization \cite{aco-dorigo, aco-specialsection, aco-blum} method to give the vast number of other genes a fair chance of inclusion; this degree of exploitation to be controlled by the user depending upon the problem at hand. As a result of this, we observed BBO-SVM and BBO-RF both to converge to an optimal or near optimal solution faster than their counterparts without heuristics. Also, the average suitability of habitats in the ecosystem almost always shows a monotonic improvement in this case unlike the earlier. In effect, it is more like the overall ecosystem (population) showing improvement. The inclusion of probabilistic selection of genes based on entropy, during mutation, adds to the improvement of the overall results. This behavior was consistently observed over 50 simulations. Figure \ref{fig_graph} shows an example of how the population average in BBO with information gain heuristics converged to higher accuracy in lesser generations as compared to simple BBO. The boxed portion demonstrates the typically observed non-monotonic behavior during some runs of simple BBO as against heuristic BBO.

With regards to the classifiers \textendash \space SVM and RF \textendash \space that we have used to evaluate BBO selected gene subsets, there has been work in literature that has reported to find SVM to outperform RF \cite{svm-rf-compare}, both on average and in majority of microarray datasets; and our results reiterate the same. We have verified this observation by training SVM and RF on the same subset of selected genes. The evaluation of each habitat to obtain its HSI from the classifier (the CV accuracies) can be run in parallel resulting in faster results by speeding up the whole process.

From literature, SVMRFE-RG \cite{geneVapnik} and Fisher-RG-SVMRFE \cite{azadeh} report high accuracies with SVM for classification. But the SVMRFE-RG is unable to tackle redundant genes \cite{azadeh}. While in \cite{azadeh}, the authors have used gene ontology to tackle redundant genes during gene selection, we have used the information gain based gene ranking. There are also other methods as in \cite{svm-furey} who have attempted classification of cancer tissue samples with SVMs without feature (gene) selection. They have reported approximately 85\% accuracy of classification (about 10-11 falsely classified samples out of nearly 70 in AML/ALL leukemia cancer case). In \cite{aco-rf}, the authors have used Ant Colony Optimization (ACO) for gene selection with Ant Miner (AM) and RF for classification.

The parameters and their corresponding tuned values used in our algorithms have been listed in Table II. These values were observed to give the most optimum results over extensive simulations.

\begin{table}[h]
\caption{Tuned algorithm parameters}
\label{table2_tuned_param}
\centering

\begin{tabular}[h]{l | r}
\multicolumn{2}{c}{$a$. For BBO with both SVM and RF} \\
\hline
Parameter&  Values  \\
\hline
 Population size (\#candidate solutions for each generation)&  50\\
 \#Generations						    &  25\\
 Mutation probability					    &  0.70\\
 Habitat modification probability		 	    &  1.00\\
 Exploitation probability during mutation (for heuristics)  &  0.55\\

\hline
\end{tabular}

\vspace{10pt}
\begin{tabular}[h]{l | r}
\multicolumn{2}{c}{$b$. For SVM} \\
\hline
 Parameter&  Values  \\
\hline
 Cost 	&  50\\
 Gamma (for Radial Basis Function as kernel)&  0.02\\
 Folds for cross-validation	& 10\\

\hline
\end{tabular}

\vspace{10pt}
\begin{tabular}[h]{l | r}
\multicolumn{2}{c}{$c$. For RF} \\
\hline
 Parameter&  Values  \\
\hline
 Trees in the forest&  500\\
 Features per tree&  $\sqrt{\scriptstyle features\_selected\_and\_fed\_to\_RF}$\\
\hline
\end{tabular}

\end{table}

\vspace{5pt}
$Results$ : Table III lists the sizes of gene subsets selected by BBO separately run with SVM and RF algorithms and the 10 \textendash \space fold cross validation accuracies obtained for the selected gene subsets.

\begin{table}[h]
\caption{Subset sizes and best 10-fold cross validation accuracies (CVAs) in \%}
\label{table3_size_accuracy}
\centering
\begin{tabular}[h]{l  c  c  c  c c}
\hline

 D&  Original& 	\#genes &	10-fold & 	 \#genes&	10-fold\\
  &   \#genes&	selected&	CVA for&	selected&	CVA for\\
  & 	   &     BBO-SVM& 	BBO-SVM&	  BBO-RF&	BBO-RF\\
\hline
 C&  2000&  	09&  		98.39&  	{11}&		92.34\\
 B&  7129& 	15&  		99.56&  	{20}&		94.38\\
 L&  7129&  	19&  		99.60&  	{20}&		93.20\\

\hline
\end{tabular}
\end{table}

With reference to literature, BBO-SVM and BBO-RF have fared well as compared to the previously best performing algorithms (for the same colon cancer dataset) namely SVMRFE-RG \cite{geneVapnik}, Fisher-RG-SVMRFE \cite{azadeh}, ACO-AM (Ant Colony Optimization\textendash Ant Miner) and ACO-RF \cite{aco-rf} which had demonstrated accuracies of 93.3, 94.7, 95.47 and 96.77\% respectively \cite{subha,liu}. Similarly, the best performing algorithms for leukemia cancer classification have shown accuracies in the range 91\textendash 97\%, with the best being 97.06\% \cite{aco-rf,geneVapnik,azadeh,roughsets,cong}. While \cite{aco-rf, geneVapnik} and \cite{azadeh} have worked with the same dataset for AML/ALL classification, \cite{roughsets} and \cite{cong} have worked with a different dataset (for Diffuse Large B-Cell Lymphoma (DLBCL)) but with similar properties, which makes us believe that our proposed algorithm with gene selection will also perform equally well as in AML/ALL classification. In case of breast cancer, the reported accuracies, for the dataset we have worked with, have been in the range of 91\textendash 94\%\cite{angela, leuk-golub}. Very clearly, our method of using BBO for gene selection in combination with SVM and RF has outperformed them with accuracies as shown in Table \ref{table3_size_accuracy}. It is worth to note that the methods in literature have almost always reported their best accuracies. In our work, we have reported the average accuracies for both, BBO\textendash SVM and BBO\textendash RF.

\section{Conclusion}

The hybrid BBO-SVM and BBO-RF techniques have shown consistently good results when compared against the highest accuracies for colon cancer, breast cancer and leukemia cancer datasets. Like other evolutionary algorithms, they are also simple to implement, robust and flexible since we can have various possible alternatives as suited to the problem and domain constraints. A significant speedup in the algorithm may be achieved by parallel implementations where the classification accuracies for individual candidate solutions may be computed in parallel.

\section{Future Work}

Like all EAs, our hybrid methods spur many possibilities of future work \textendash \space some problem dependent and some from the implementation perspective. From problem representation to specific migration and mutation strategies, we can have a variety of schemes. For example, with respect to problem representation, the other suitable scheme that one could explore for BBO here is: each habitat of the ecosystem could have an arbitrary number of attributes (SIVs) which remains fixed for itself across generations but may differ from other habitats in the ecosystem. This is similar to the variable population sizes proposed earlier for BBO  \cite{simonbbo} and other EAs \cite{var-pop-size-ga} but at a finer granularity, variable solution (habitats or chromosomes as the case may be) sizes. This could lead us to a compact implementation framework that can simultaneously output the better performing subset sizes along with the selected genes. Many such variations possible with other EAs could work here too.

\section*{Acknowledgment}

The authors acknowledge the Centre for Modeling and Simulation, University of Pune, India for its support. Also, VKJ gratefully acknowledges the Department of Science and Technology (DST), New Delhi, India for financial support.


\section{Introduction}\label{sec:Intro}
Strong stability preserving time discretization methods were originally
developed for the solution of nonlinear hyperbolic
partial differential equations (PDEs).  Solutions of such PDEs may
contain discontinuities even when the initial conditions are smooth.
Many numerical methods for their solution are based on a method-of-lines approach
in which the problem is first discretized in space to yield a system of ODEs.
The spatial discretization is often chosen to ensure the solution is total variation diminishing (TVD),
in order to avoid the appearance of spurious oscillations near discontinuities,
\emph{when coupled with first-order forward Euler time integration}.
Strong stability preserving (SSP) time discretizations (also known as TVD
discretizations \cite{Gottlieb/Shu:1998}) are high-order time
discretizations that guarantee the TVD property (or other convex functional bounds), with a
possibly different step-size restriction \cite{Gottlieb2011a}.
Section~\ref{sec:SSP} reviews Runge--Kutta methods and the concept of
strong stability preserving methods.

Explicit SSP Runge--Kutta methods cannot have order greater
than four \cite{Ruuth2002}.  However, a Runge--Kutta method may
achieve an effective order of accuracy higher than its classical order
by the use of special starting and stopping procedures.
The conditions for a method to have effective order $q$ are
in general less restrictive than the conditions for a method
to have classical order $q$.
Section~\ref{sec:Algebraic_RK} presents a brief overview of the
algebraic representation of Runge--Kutta methods, following Butcher
\cite{Butcher2008_book}.
This includes the concept of effective order and a list of effective
order conditions.

We examine the SSP properties of explicit Runge--Kutta methods whose
effective order is greater than their classical order.
Previous studies of SSP Runge--Kutta methods have considered only the
classical order of the methods.
Three natural questions are:
\begin{itemize}
    \item Can an SSP Runge--Kutta method have effective order of accuracy greater than four?
    \item If we only require methods to have {\em effective} order $q$, is it possible to achieve larger
            SSP coefficients than those obtained in methods with {\em classical} order $q$?
    \item SSP Runge--Kutta methods of order four require at least five stages.  Can SSP methods of
          effective order four have fewer stages?

\end{itemize}
We show in Section \ref{sec:ExRK_barrier} that the answer to the first question is negative.
We answer the second question by numerically solving the problem of optimizing
the SSP coefficient over the class of methods with effective order $q$;
see Section~\ref{sec:optimal_ESSPRK}.
Most of the methods we find are shown to be optimal, as they achieve a certain theoretical
upper bound on the SSP coefficient that is obtained by considering only
linear problems \cite{Kraaijevanger1986}.
We answer the last question affirmatively
by construction, also in Section~\ref{sec:optimal_ESSPRK}.
The paper concludes with numerical experiments in
Section~\ref{sec:numerics} and conclusions in
Section~\ref{sec:Conclusion}.

\section{Strong stability preserving Runge--Kutta methods}\label{sec:SSP}
Strong stability preserving (SSP) time-stepping methods were originally introduced
for time integration of systems of hyperbolic conservation laws
\cite{Shu/Osher:1988}
\begin{align}\label{eq:pde}
	\bm{U}_t + \nabla \cdot \bm{f}(\bm{U}) = 0,
\end{align}
with appropriate initial and boundary conditions.
A spatial discretization gives the system of ODEs
\begin{align}\label{eq:ode_system}
    \bm{u}'(t) = \bm{F}(\bm{u}(t)),
\end{align}
where $\bm{u}$ is a vector of continuous-in-time grid values approximating
the solution $\bm{U}$ at discrete grid points.
Of course, \eqref{eq:ode_system} can arise in many ways and $\bm{F}$
need not necessarily represent a spatial discretization.
Particularly, $\bm{F}$ may be time-dependent, but we can always make a
transformation to an autonomous form.
In any case, a time discretization then produces a sequence of
solutions $\bm{u}^{n} \approx \bm{u}(t_n)$.
This work studies explicit Runge--Kutta time discretizations.
An explicit $s$-stage Runge--Kutta method takes the form
\begin{align*}
	\bm{u}^{n+1} &= \bm{u}^{n} + \Delta t \sum_i^s b_i \bm{F}(\bm{Y}_i),
\end{align*}
where
\begin{align*}
	\bm{Y}_i &= \bm{u}^{n} + \Delta t \sum_j^{i-1} a_{ij} \bm{F}(\bm{Y}_j).
\end{align*}
Such methods are characterized by the coefficient matrix $A = (a_{ij}) \in
\mathbb{R}^{s \times s}$, the weight vector $\bm{b} = (b_i) \in \mathbb{R}^s$
and the abscissa $\bm{c} = (c_i) \in \mathbb{R}^s$, where
$c_i = \sum_{j=1}^{i-1}a_{ij}$.
The accuracy and stability of the method depend on the coefficients of the
Butcher tableau $(A,\bm{b},\bm{c})$ \cite{Butcher2008_book}.

In some cases, the solutions of hyperbolic conservation laws satisfy a
monotonicity property. For example, if \eqref{eq:pde} is scalar then solutions
are monotonic in the total variation semi-norm \cite{Ketcheson2008}.
For this reason, many popular spatial discretizations are designed such
that, for a suitable class of problems, the solution $\bm{u}$ in
\eqref{eq:ode_system} computed with the forward Euler scheme is
non-increasing (in time) in some norm, semi-norm, or convex functional; i.e.,
\begin{align}\label{eq:forwardEuler}
    \|\bm{u} + \Delta t\bm{F}(\bm{u})\| \le \|\bm{u}\|, \quad \text{for all } \bm{u} \text{ and for } 0 \le \Delta t \le \Delta t_{\textnormal{FE}}.
\end{align}
If this is the case, then an SSP method also generates a solution whose norm is
non-increasing in time, under a modified time-step restriction.
\begin{definition}[Strong Stability Preserving]
	A Runge--Kutta method is said to be \emph{strong stability preserving} with
	\emph{SSP coefficient} $\mathcal{C} > 0$ if, whenever the forward Euler condition
	\eqref{eq:forwardEuler} holds and
	\begin{align*}
		0 \leq \Delta t \leq \mathcal{C} \Delta t_{\textnormal{FE}},
	\end{align*}
	the Runge--Kutta method generates a monotonic sequence of solution values $\bm{u}^n$ satisfying
	\begin{align*}
  		\|\bm{u}^{n+1}\| \le \|\bm{u}^n\|.
	\end{align*}
\end{definition}

Note that $\Delta t_{\textnormal{FE}}$ is a property of the spatial discretization $\bm{F}$
and is independent of $\bm{u}$.
The SSP coefficient $\mathcal{C}$ is a property of the particular
time-stepping method and quantifies the allowable time step size relative
to that of the forward Euler method.
Generally we want the SSP coefficient to be as large as possible for efficiency.
To allow a fair comparison of explicit methods with different number of stages,
we consider the \emph{effective SSP coefficient}
\begin{align*}
	\sspcoef_{\textnormal{eff}} = \frac{\mathcal{C}}{s}.
\end{align*}
Note that the use of the word \emph{effective} here is unrelated to the
concept of \emph{effective order} introduced in Section~\ref{sec:Algebraic_RK}.

\subsection{Optimal SSP schemes}\label{subsec:Optimal_SSPRK}
We say that an SSP Runge--Kutta method is optimal if it has the largest
possible SSP coefficient for a given order and a given number of stages.
The search for these optimal methods was originally based on
expressing the Runge--Kutta method as combinations of forward Euler
steps (the Shu--Osher form) and solving a non-linear optimization
problem \cite{Gottlieb/Shu:1998, Gottlieb2001, Spiteri2003a, Spiteri2003b,
Ruuth2004, Ruuth:2006}.
However, the SSP coefficient is related to the
\emph{radius of absolute monotonicity} \cite{Kraaijevanger1991} and,
for irreducible Runge--Kutta methods, the two are equivalent
\cite{Ferracina2004, Higueras2004}.
This gives a simplified algebraic characterization of the SSP coefficient
\cite{Ferracina2005}; it is the maximum value of $r$ such that the following
conditions hold:
\begin{subequations} \label{eq:absmon}
\begin{align}
    K(I + rA)^{-1} \geq 0 \\
    \bm{e}_{s+1} - rK(I + rA)^{-1}\bm{e}_{s} \geq 0,
\end{align}
\end{subequations}
provided that $I + rA$ is invertible.
Here
\begin{equation*}
    K = \left(
            \begin{array}{c}
                     A              \\
                     \bm{b}^{\texttt{T}}
            \end{array}
         \right),
\end{equation*}
while $\bm{e}_s$ denotes the vector of ones of length $s$ and $I$ is the
$s \times s$ identity matrix.
The inequalities are understood component-wise.

The optimization problem of finding optimal SSP Runge--Kutta methods
can thus be written as follows:
\begin{equation}\label{eq:SSP_opt}
    \max_{A, \bm{b}, r} \; r \quad \text{subject to \eqref{eq:absmon} and } \Phi(K) = 0.
\end{equation}
Here \( \Phi(K) \) represents the  order conditions.

Following \cite{Ketcheson2008, Ketcheson/Macdonald/Gottlieb:2009},
we will numerically solve the optimization problem \eqref{eq:SSP_opt} to find
optimal explicit SSP Runge--Kutta methods for various effective orders of accuracy.
However, we first need to define the order conditions $\Phi(K)$ for these methods.
This is discussed in the next section.

\section{The effective order of Runge--Kutta methods}\label{sec:Algebraic_RK}

The definition, construction, and application of methods with an
effective order of accuracy relies on the use of starting and stopping
methods.
Specifically, we consider a \emph{starting method} $S$, a \emph{main
  method} $M$, and a \emph{stopping method} $S^{-1}$.
The successive use of these three methods results in a method $P =
S^{-1}MS$, which denotes the application of method $S$, followed by
method $M$, followed by method $S^{-1}$.
The method $S^{-1}$ is an ``inverse'' of method $S$.
We want $P$ to have order $q$, whereas $M$ might have lower classical
order $p < q$.
We then say $M$ has \emph{effective order} $q$.

When the method $P$ is used for $n$ steps,
$$P^n = (S^{-1}MS)^n = (S^{-1}MS) \cdots (S^{-1}MS) (S^{-1}MS),$$
it turns out that only $M$ need be used repeatedly, as in
$S^{-1} M^n S$,
because
$S S^{-1}$ leaves the solution unchanged up to order $q$.
The starting method introduces a perturbation to the solution,
followed by $n$ time steps of the main method $M$, and finally the
stopping method is used to correct the solution.
In Section~\ref{subsec:starting_stopping}, we propose alternative
starting and stopping procedures which allow the overall procedure to
be SSP.

The effective order of a Runge--Kutta method is defined in an abstract
algebraic context introduced by Butcher \cite{Butcher1969} and developed
further in \cite{Butcher1972, Hairer1974, Butcher1996, Butcher1998} and
others.
We follow the book \cite{Butcher2008_book} in our
derivation of the effective order conditions.

\subsection{The algebraic representation of Runge--Kutta methods}\label{subsec:Algebraic_representation}

\begin{table}
	\caption{Elementary weights $\alpha(t_i)$ of trees $t_i$ up to order five for a
  		Runge--Kutta method with Butcher tableau $(A,\bm{b},\bm{c})$.
  		Here $C$ is a diagonal matrix with components
  		$c_{i} = \sum_{j=1}^{i-1} a_{ij}$ and exponents of vectors
  		represent component exponentiation.
  		By convention $\alpha_0 = \alpha(t_{0}) = 1$, where $t_{0}$
  		denotes the empty tree.}
	\centering
	\begin{smalltrees}
		\begin{tabular}{ccccccccc}
    		\cline{1-4}\cline{6-9}
                \noalign{\vskip 2pt}
    		$i$ & tree $t_i$ & $\alpha(t_i)$ & $\gamma(t_i)$ & & $i$ & tree $t_i$ & $\alpha(t_i)$ & $\gamma(t_i)$ \\
                \noalign{\vskip 1pt}
    		\cline{1-4}\cline{6-9}
                \noalign{\vskip 3pt}
    		0 & $\emptyset$ \hspace{15pt}  & 1 & 0 & & 9 & \hspace{15pt} \tree{9} & $\bm{b}^T\bm{c}^4$ & 5 \\
    		1 & \hspace{15pt} \tree{1} & $\bm{b}^T\bm{e}$ & 1 & & 10 & \tree{10} \hspace{15pt} & $\bm{b}^TC^2A\bm{c}$ & 10 \\
    		2 & \tree{2} \hspace{15pt}  & $\bm{b}^T\bm{c}$ & 2 & & 11 & \hspace{15pt} \tree{11} & $\bm{b}^TCA\bm{c}^2$ & 15 \\
    		3 & \hspace{15pt} \tree{3} & $\bm{b}^T\bm{c}^2$ & 3 & & 12 & \tree{12} \hspace{15pt} & $\bm{b}^TCA^2\bm{c}$ & 30 \\
    		4 & \tree{4} \hspace{15pt}  & $\bm{b}^TA\bm{c}$ & 6 & & 13 & \hspace{15pt} \tree{13} & $\bm{b}^T(A\bm{c})^2$ & 20 \\
    		5 & \hspace{15pt} \tree{5} & $\bm{b}^T\bm{c}^3$ & 4 & & 14 & \tree{14} \hspace{15pt} & $\bm{b}^TA\bm{c}^3$ & 20 \\
    		6 & \tree{6} \hspace{15pt}  & $\bm{b}^TCA\bm{c}$ & 8 & & 15 & \hspace{15pt} \tree{15} & $\bm{b}^TACA\bm{c}$ & 40 \\
    		7 & \hspace{15pt} \tree{7} & $\bm{b}^TA\bm{c}^2$ & 12 & & 16 & \tree{16} \hspace{15pt} & $\bm{b}^TA^2\bm{c}^2$ & 60 \\
    		8 & \tree{8} \hspace{15pt}  & $\bm{b}^TA^2\bm{c}$ & 24 & & 17 & \hspace{15pt} \tree{17} & $\bm{b}^TA^3\bm{c}$ & 120 \tabularnewline\noalign{\vskip 2pt}
		\noalign{\vskip 3pt}
    		\cline{1-4}\cline{6-9}
  		\end{tabular}
  \end{smalltrees}
	\label{tab:elementary_weights}
\end{table}

According to Butcher's algebraic theory, irreducible Runge--Kutta methods
are placed in one-to-one correspondence with elements of a group
$G$, consisting of real-valued functions on the set of rooted trees \cite[Theorem~384A]{Butcher2008_book}.
A Runge--Kutta method corresponds to the map that takes each rooted tree $t$
to the corresponding elementary weight $\Phi(t)$ of that Runge--Kutta method.
Table~\ref{tab:elementary_weights} lists the elementary weights for trees of
up to degree five; a general recursive formula can be found in
\cite[Definition~312A]{Butcher2008_book}.
The ordering of trees given in Table~\ref{tab:elementary_weights} is used
throughout the remainder of this work; thus $t_9$ refers to the tree with
elementary weight $\bm{b}^T \bm{c}^4$.
For a function $\alpha \in G$ we write the values of the
elementary weights as $\alpha_{i} = \alpha(t_{i})$ for tree $t_{i}$.
A special element of the group $E \in G$ corresponds to the
(hypothetical) method that evolves the solution exactly.
The values of $E(t)$ are denoted $1/\gamma(t)$ \cite{Butcher2008_book}
and the values of $\gamma(t)$ are included in
Table~\ref{tab:elementary_weights}.
Classical order conditions are obtained by comparing the elementary weights
of a method with these values.

Let $\alpha, \beta \in G$ correspond to Runge--Kutta methods $M_1$ and $M_2$
respectively.
The application of method $M_1$ followed by method $M_2$ corresponds to
the multiplicative group operation $\alpha\beta$.\footnote{We write
	$M_2M_1$ to mean the application of $M_1$
	followed by the application of $M_2$
	(following matrix and operator ordering convention)
	but when referring to products of elements of $G$
        we use the reverse ordering ($\alpha\beta$)
	to match the convention in \cite{Butcher2008_book}.}
This product is defined by partitioning the input tree and computing
over the resulting forest \cite[\S~383]{Butcher2008_book}.

Two Runge--Kutta methods $M_1$ and $M_2$, are equivalent up to order
$p$ if their corresponding elements in $G$, $\alpha$ and $\beta$, satisfy
$\alpha(t) = \beta(t)$, for every tree $t$ with $r(t) \leq p$,
where $r(t)$ denotes the order of the tree (number of vertices).
We denote this equivalence relation by
$$M_1 \equival{p} M_2.$$
In this sense, methods have inverses: the product of $\alpha^{-1}$ and
$\alpha$ must match the identity method up to order $p$.
Note that inverse methods up to order $p$ are not unique and inverse
methods of explicit methods need not be implicit.
We can then define the effective order of accuracy of a method $M$
with starting method $S$ and stopping method $S^{-1}$.
\begin{definition}\cite[\S~389]{Butcher2008_book}\label{def:Effective_order}
  Suppose $M$ is a Runge--Kutta method with corresponding $\alpha \in G$.
  Then the method $M$ is of effective order $q$ if there exist methods
  $S,S^{-1}$ (with corresponding $\beta, \beta^{-1} \in G$) such that
	\begin{equation}\label{eq:Effective_order_1}
		(\beta\alpha\beta^{-1})(t) = E(t), \; \text{for every tree with } r(t) \leq q,
	\end{equation}
        where $\beta^{-1}$ is an inverse of $\beta$ up to order $q$; i.e.
        \begin{align*}
        		(\beta^{-1}\beta)(t) = 1(t), \; \text{for every tree with } r(t) \leq q.
		\end{align*}
        Here $1\in G$ is the identity element
        and $E\in G$ is the exact evolution operator.
\end{definition}

\begin{table}
    \caption{Effective order five conditions on $\alpha$ (main
      method $M$) in terms of order conditions on $\beta$
      (starting method $S$).
      See also \cite[\S~389]{Butcher2008_book}.
      Recall that $\alpha_i$ and $\beta_i$
      are the elementary weights associated with the index $i$ in
      Table~\ref{tab:elementary_weights}.
      We assume that $\beta_1=0$ (see
      Section~\ref{subsubsec:Main_starting_conditions}).}
  \small
  \setlength{\extrarowheight}{0.5pt}
  \centering

  \begin{tabular}{p{2mm}p{51mm}p{64mm}}
    \toprule
    $q$  &  Effective order conditions \\
    \midrule\noalign{\vskip 1.5pt}
    $1$  &
            $\alpha_1  = 1$. \\
    \mydashrule
    $2$  &
            $\alpha_2  = \tfrac{1}{2}$. \\
   \mydashrule
    $3$  &
            $\alpha_3  = \tfrac{1}{3} + 2\beta_2$,  \quad
            $\alpha_4  = \tfrac{1}{6}$.  \\
    \mydashrule
    $4$  &  \multicolumn{2}{l}{%
            $\alpha_5  = \tfrac{1}{4} + 3\beta_2 + 3\beta_3$, \hfill
            $\alpha_6  = \tfrac{1}{8} + \beta_2 + \beta_3 + \beta_4$, \hfill
            $\alpha_7  = \tfrac{1}{12} +\beta_2 - \beta_3 + 2\beta_4$, \hfill
            $\alpha_8  = \tfrac{1}{24}$.} \\
    \mydashrule
    $5$  &
            $\alpha_9  = \tfrac{1}{5} + 4\beta_2 + 6\beta_3 + 4\beta_5$,
         &  \hspace*{-3pt}$\alpha_{10} = \tfrac{1}{10} + \tfrac{5}{3}\beta_2 - 2\beta_2^{2} + \tfrac{5}{2}\beta_3 + \beta_4 + \beta_5 + 2\beta_6$, \tabularnewline\noalign{\vskip 2pt}
         &  $\alpha_{11} = \tfrac{1}{15} + \tfrac{4}{3}\beta_2 + \tfrac{1}{2}\beta_3 + 2\beta_4 + 2\beta_6 + \beta_7$,
         &  \hspace*{-3pt}$\alpha_{12} = \tfrac{1}{30} + \tfrac{1}{3}\beta_2 - 2\beta_2^{2} + \tfrac{1}{2}\beta_3 + \tfrac{1}{2}\beta_4 + \beta_6 + \beta_8$, \tabularnewline\noalign{\vskip 2pt}
         &  $\alpha_{13} = \tfrac{1}{20} + \tfrac{2}{3}\beta_2 - \beta_2^{2} + \beta_3 + \beta_4 + 2\beta_6$,
         &  \hspace*{-3pt}$\alpha_{14} = \tfrac{1}{20} + \beta_2 + 3\beta_4 - \beta_5 + 3\beta_7$, \tabularnewline\noalign{\vskip 2pt}
         &  $\alpha_{15} = \tfrac{1}{40} + \tfrac{1}{3}\beta_2 + \tfrac{3}{2}\beta_4 - \beta_6 + \beta_7 + \beta_8$,
         &  \hspace*{-3pt}$\alpha_{16} = \tfrac{1}{60} + \tfrac{1}{3}\beta_2 - \tfrac{1}{2}\beta_3 + \beta_4 - \beta_7 + 2\beta_8$, \,\! $\alpha_{17} = \tfrac{1}{120}$. \tabularnewline\noalign{\vskip 2pt}
            \bottomrule
    \end{tabular}
    \label{tab:effective_OCs_on_alpha}
\end{table}

\subsection{Effective order conditions}\label{sec:effOrderCond}
For the main method $M$ to have effective order $q$, its coefficients
and those of the starting and stopping methods must satisfy a set of
algebraic conditions.
These \emph{effective order conditions} can be found
by rewriting \eqref{eq:Effective_order_1} as
$(\beta\alpha)(t) = (E\beta)(t)$
and applying the group product operation.
For trees up to order five these are tabulated in Table~\ref{tab:effective_OCs_on_alpha} (and also in \cite[\S~389]{Butcher2008_book}).
In general, the effective order conditions allow more degrees of
freedom for method design than do the classical order conditions.
Note that the effective order conditions match the classical order conditions up to
second order.
\begin{remark}\label{rem:talltrees}
	The effective order conditions of the main method for the ``tall" trees
	$t_1, t_2, t_4, t_8, t_{17}, \dots$ match the classical order conditions
	and these are precisely the order conditions for linear problems.
	This follows from inductive application of the group product
	on the tall trees.
	Therefore, methods of effective order $q$ have classical order at least $q$
	for linear problems.
\end{remark}

\subsubsection{Order conditions of the main and starting methods}\label{subsubsec:Main_starting_conditions}
As recommended in \cite{Butcher2008_book},
we consider the elementary weights $\beta_{i}$ of the starting method as free
parameters when determining the elementary weights $\alpha_i$ of the main method.
The relationship in Table~\ref{tab:effective_OCs_on_alpha} between the
$\alpha_i$ and $\beta_i$ is mostly linear (although there are a few
$\beta_2^2$ terms).
It is thus straightforward to (mostly) isolate the equations for $\alpha_i$
and determine the $\beta_i$ as linear combination of the $\alpha_i$.
This separation provides maximal degrees of freedom and minimizes the number of
constraints when constructing the method $M$.
The resulting effective order conditions for the main method $M$ are given
in Table~\ref{tab:effective_OCs} (up to effective order five).
For a specified classical and effective order, these are the equality constraints
$\Phi(K)$ in the optimization problem \eqref{eq:SSP_opt} for method $M$.

Constructing the main method $M$ then determines the $\alpha$ values
and we obtain a set of order conditions on $\beta$ (for that
particular choice of $M$).  These are given in the right-half of
Table~\ref{tab:effective_OCs}.
We can also find the order conditions on $S^{-1}$ in terms of the
$\beta_i$ (see \cite[Table~386(III)]{Butcher2008_book}).
We note that increasing the classical order of the main method requires
$\alpha_i = 1/\gamma(t_i)$ and thus by Table~\ref{tab:effective_OCs_on_alpha} requires more of the $\beta_i$ to be zero.

\begin{table}
    \caption{Effective order $q$, classical order $p$ conditions on $ \alpha $ and $ \beta $ for the main and starting methods, $M$ and $S$ respectively.}
  \small
  \centering

    \begin{tabular}{M{2mm}M{2mm}M{68mm}M{67mm}}
      \toprule
        $q$ & $p$ & Order conditions for main method $M$ & Order conditions for starting method $S$ \tabularnewline\noalign{\vskip 2pt}
      \midrule
        \multirow{1}{*}{3} & \multirow{1}{*}{2} & {\small $\alpha_1 = 1$, $\alpha_2 = \frac{1}{2}$, $\alpha_4 = \frac{1}{6}$.} & {\small $\beta_1 = 0$, $\beta_2 = - \frac{1}{6} + \frac{1}{2}\alpha_3$.}\tabularnewline\noalign{\vskip 2pt}
      \mydashrule
        \multirow{3}{*}{4} & \multirow{3}{*}{2} & {\small $\alpha_1 = 1$, $\alpha_2 = \frac{1}{2}$, $\alpha_4 = \frac{1}{6}$,} & {\small $\beta_1 = 0$, $\beta_2 = - \frac{1}{6} + \frac{1}{2}\alpha_3$,}\tabularnewline\noalign{\vskip 2pt}
        & & {\small $\frac{1}{4} - \alpha_3 + \alpha_5 - 2\alpha_6 + \alpha_7 = 0$, $\alpha_8 = \frac{1}{24}$.} & {\small $\beta_3 = \frac{1}{12} - \frac{1}{2}\alpha_3 + \frac{1}{3}\alpha_5$, $\beta_4 = - \frac{1}{24} - \frac{1}{3}\alpha_5 + \alpha_6$.} \tabularnewline\noalign{\vskip 2pt}
      \mydashrule
        \multirow{3}{*}{4} & \multirow{3}{*}{3} & {\small $\alpha_1 = 1$, $\alpha_2 = \frac{1}{2}$, $\alpha_3 = \frac{1}{3}$, $\alpha_4 = \frac{1}{6}$,} & {\small $\beta_1 = 0$, $\beta_2 = 0$, $\beta_3 = - \frac{1}{12}  + \frac{1}{3}\alpha_5$,} \tabularnewline\noalign{\vskip 2pt}
        & & {\small $\frac{1}{12} - \alpha_5 + 2\alpha_6 - \alpha_7 = 0$, $\alpha_8 = \frac{1}{24}$.} & {\small $\beta_4 = - \frac{1}{24} - \frac{1}{3}\alpha_5 + \alpha_6$.} \tabularnewline\noalign{\vskip 2pt}
      \mydashrule
        \multirow{8}{*}{5} & \multirow{8}{*}{2} & {\small $\alpha_1 = 1$, $\alpha_2 = \frac{1}{2}$, $\alpha_4 = \frac{1}{6}$, $\alpha_8 = \frac{1}{24}$, $\alpha_{17} = \frac{1}{120}$,} & {\small $\beta_1 = 0$, $\beta_2 = - \frac{1}{6} + \frac{1}{2}\alpha_3$,} \tabularnewline\noalign{\vskip 2pt}
        & & {\small $\frac{1}{4} - \alpha_3 + \alpha_5 - 2\alpha_6 + \alpha_7 = 0$,} & {\small $\beta_3 = \frac{1}{12} - \frac{1}{2}\alpha_3 + \frac{1}{3}\alpha_5$, $\beta_4 = -\frac{1}{24} - \frac{1}{3}\alpha_5 + \alpha_6$} \tabularnewline\noalign{\vskip 2pt}
        & & {\small $\frac{1}{4}\alpha_9-\alpha_{10}+\alpha_{13}=\beta_2^{2}$, \: $\beta_2 = - \frac{1}{6} + \frac{1}{2}\alpha_3$,} & {\small $\beta_5 = -\frac{1}{120} + \frac{1}{4}\alpha_3 - \frac{1}{2}\alpha_5 + \frac{1}{4}\alpha_9$,} \tabularnewline\noalign{\vskip 2pt}
        & & {\small $\frac{3}{10} - \frac{3}{2}\alpha_3 + \alpha_5 + \frac{1}{2}\alpha_9 - 3\alpha_{10} + 3\alpha_{11} - \alpha_{14} = 6\beta_2^{2}$,} & {\small $\beta_6 = \frac{7}{720} + \beta_2^{2} + \frac{1}{12}\alpha_3 - \frac{1}{2}\alpha_6 - \frac{1}{8}\alpha_9 + \frac{1}{2}\alpha_{10}$,} \tabularnewline\noalign{\vskip 2pt}
        & & {\small $\frac{1}{15} - \frac{1}{2}\alpha_3 + \alpha_6 + \frac{1}{2}\alpha_9 - 2\alpha_{10} + \alpha_{11} + \alpha_{12} - \alpha_{15} = 2\beta_2^{2}$,} & {\small $\beta_7 = \frac{8}{45} - 2\beta_2^{2} - \frac{7}{12}\alpha_3 + \frac{1}{2}\alpha_5 - \alpha_6 + \frac{1}{4}\alpha_9 - \alpha_{10} + \alpha_{11}$,} \tabularnewline\noalign{\vskip 2pt}
        & & {\small $\frac{19}{60} - \alpha_3 + \alpha_5 - 2\alpha_6 + \alpha_{11} - 2\alpha_{12} + \alpha_{16} = 4\beta_2^{2}$.} & {\small $\beta_8 = -\frac{1}{120} + \beta_2^{2} + \frac{1}{8}\alpha_9 - \frac{1}{2}\alpha_{10} + \alpha_{12}$.} \tabularnewline\noalign{\vskip 2pt}
      \mydashrule
        \multirow{7}{*}{5} & \multirow{7}{*}{3} & {\small $\alpha_1 = 1$, $\alpha_2 = \frac{1}{2}$, $\alpha_3 = \frac{1}{3}$, $\alpha_4 = \frac{1}{6}$, $\alpha_8 = \frac{1}{24}$,} & {\small $\beta_1 = 0$, $\beta_2 = 0$, $\beta_3 = -\frac{1}{12} + \frac{1}{3}\alpha_5$} \tabularnewline\noalign{\vskip 2pt}
        & & {\small $\alpha_{17} = \frac{1}{120}$, $\frac{1}{12} - \alpha_5 + 2\alpha_6 - \alpha_7 = 0$,} & {\small $\beta_4 = -\frac{1}{24} - \frac{1}{3}\alpha_5 + \alpha_6$,} \tabularnewline\noalign{\vskip 2pt}
        & & {\small $\frac{1}{4}\alpha_9 - \alpha_{10} + \alpha_{13} = 0$,} & {\small $\beta_5 = \frac{3}{40} - \frac{1}{2}\alpha_5 + \frac{1}{4}\alpha_9$,} \tabularnewline\noalign{\vskip 2pt}
        & & {\small $\frac{1}{5} - \alpha_5 - \frac{1}{2}\alpha_9 + 3\alpha_{10} - 3\alpha_{11} + \alpha_{14} = 0$,} & {\small $\beta_6 = \frac{3}{80} - \frac{1}{2}\alpha_6 - \frac{1}{8}\alpha_9 + \frac{1}{2}\alpha_{10}$,} \tabularnewline\noalign{\vskip 2pt}
        & & {\small $\frac{1}{10} - \alpha_6 - \frac{1}{2}\alpha_9 + 2\alpha_{10} - \alpha_{11} - \alpha_{12} + \alpha_{15} = 0$,} & {\small $\beta_7 = -\frac{1}{60} + \frac{1}{2}\alpha_5 - \alpha_6 + \frac{1}{4}\alpha_9 - \alpha_{10} + \alpha_{11}$,} \tabularnewline\noalign{\vskip 2pt}
        & & {\small $\frac{1}{60} - \alpha_5 + 2\alpha_6 - \alpha_{11} + 2\alpha_{12} - \alpha_{16} = 0$.} & {\small $\beta_8 = -\frac{1}{120} + \frac{1}{8}\alpha_9 - \frac{1}{2}\alpha_{10} + \alpha_{12}$.} \tabularnewline\noalign{\vskip 2pt}
      \mydashrule
        \multirow{7}{*}{5} & \multirow{7}{*}{4} & {\small $\alpha_1 = 1$, $\alpha_2 = \frac{1}{2}$, $\alpha_3 = \frac{1}{3}$, $\alpha_4 = \frac{1}{6}$, $\alpha_5 = \frac{1}{4}$,} & {\small $\beta_1 = 0$, $\beta_2 = 0$,} \tabularnewline\noalign{\vskip 2pt}
        & & {\small $\alpha_6 = \frac{1}{8}$, $\alpha_7 = \frac{1}{12}$, $\alpha_8 = \frac{1}{24}$, $\alpha_{17} = \frac{1}{120}$,} & {\small $\beta_3 = 0$, $\beta_4 = 0$,} \tabularnewline\noalign{\vskip 2pt}
        & & {\small $\frac{1}{4}\alpha_9 - \alpha_{10} + \alpha_{13} = 0$,} & {\small $\beta_5 = -\frac{1}{20} + \frac{1}{4}\alpha_9$,} \tabularnewline\noalign{\vskip 2pt}
        & & {\small $\frac{1}{20} + \frac{1}{2}\alpha_9 - 3\alpha_{10} + 3\alpha_{11} - \alpha_{14} = 0$,} & {\small $\beta_6 = -\frac{1}{40} - \frac{1}{8}\alpha_9 + \frac{1}{2}\alpha_{10}$,} \tabularnewline\noalign{\vskip 2pt}
        & & {\small $\frac{1}{40} + \frac{1}{2}\alpha_9 - 2\alpha_{10} + \alpha_{11} + \alpha_{12} - \alpha_{15} = 0$,} & {\small $\beta_7 = -\frac{1}{60} + \frac{1}{4}\alpha_9 - \alpha_{10} + \alpha_{11}$,} \tabularnewline\noalign{\vskip 2pt}
        & & {\small $\frac{1}{60} - \alpha_{11} + 2\alpha_{12} - \alpha_{16} = 0$.} & {\small $\beta_8 = -\frac{1}{120} + \frac{1}{8}\alpha_9 - \frac{1}{2}\alpha_{10} + \alpha_{12}$.} \tabularnewline\noalign{\vskip 2pt}
        \bottomrule
    \end{tabular}
    \label{tab:effective_OCs}
\end{table}

Tables~\ref{tab:effective_OCs_on_alpha} and~\ref{tab:effective_OCs}
both assume that $\beta_1=0$ (i.e., the starting and stopping methods
perturb the solution but do not advance the solution in time).  This
assumption is without loss of generality following \cite[Lemma
389A]{Butcher2008_book}, the proof of which shows that we can always find
starting procedures with $\beta_1 = 0$ for which the main method has
effective order $q$, whenever this holds for a starting method with
$\beta_1 \neq 0$.

\section{Explicit SSP Runge--Kutta methods have effective order at most four}\label{sec:ExRK_barrier}
The classical order of any explicit SSP Runge--Kutta method cannot be greater
than four \cite{Ruuth2002}.
It turns out that the effective order of any explicit SSP Runge--Kutta method
also cannot be greater than four, although the proof of this result is more
involved.
We begin by recalling a well-known result.
\begin{lemma}\label{lem:positive_b}(see \cite[Theorem~4.2]{Kraaijevanger1991},\cite[Lemma 4.2]{Ruuth2002})
	Any irreducible Runge--Kutta method with positive SSP coefficient $\mathcal{C}>0$
	must have positive weights $\bm{b} > \bm{0}$.
\end{lemma}

Irreducibility \cite{dahlquist2006} is technically important in this
result and those that follow because a reducible SSP method might not
have positive weights (but it would be reducible to one that does, as
per the lemma).
The main result of this section is:
\begin{theorem}\label{thm:effective_barrier}
	Any explicit Runge--Kutta method with positive weights $\bm{b} > \bm{0}$
	has effective order at most four.
\end{theorem}

\begin{proof
	Any method of effective order five must have classical order at least two
	(see \cite{Butcher2008_book} or Table~\ref{tab:effective_OCs}).
    Thus it is sufficient to show that any method with all positive weights
    cannot satisfy the conditions of effective order five and classical order two.

    Let $(A,\bm{b},\bm{c})$ denote the coefficients of an explicit Runge--Kutta method with
    effective order at least five, classical order at least two, and positive
    weights $\bm{b} > \bm{0}$.
    The effective order five and classical order two conditions
    (see Table~\ref{tab:effective_OCs} with $q=5$ and $p=2$) include the following:
    \begin{subequations}\label{eq:theorem_cond}
    		\begin{align}
    			\bm{b}^T\bm{e} & = 1, \label{eq:theorem_cond_a} \\
             	\bm{b}^TA\bm{c} &= \frac{1}{6}, \label{eq:theorem_cond_b} \\
            	\frac{1}{2}\bm{b}^T\bm{c}^2 - \frac{1}{6} &= \beta_2, \label{eq:theorem_cond_c} \\
            	\frac{1}{4}\bm{b}^T\bm{c}^4 - \bm{b}^TC^2A\bm{c} + \bm{b}^T(A\bm{c})^2 &= \beta_2^2, \label{eq:theorem_cond_d}
        	\end{align}
	\end{subequations}
	where the powers on vectors are understood component-wise.
	Let

		$\bm{v} = \frac{1}{2}\bm{c}^2 - A\bm{c}$.

	Then substituting \eqref{eq:theorem_cond_b} into \eqref{eq:theorem_cond_c} and expressing \eqref{eq:theorem_cond_d} in terms of $\bm{v}$ gives

	\begin{align*}
		\bm{b}^T\bm{v} &= \beta_2, \\% \label{eq:btv} \\

		\bm{b}^T\bm{v}^2 &= \beta_2^2.
	\end{align*}

  Each of
  these is a strictly convex combination.

  Jensen's inequality (with a strictly convex function, as is the case with the
  square function here) then states $\bm{b}^T \bm{v}^2 \leq (\bm{b}^T
  \bm{v})^2$ with equality if and only if all components of $\bm{v}$ are equal

  \cite[Theorem 12, pg 31]{Bullen:inequalities}.

  Now $v_1 = 0$ for every explicit method so we deduce that $\bm{v}=\bm{0}$.
  That implies the method has stage order two, which is not possible for
  explicit methods \cite{Ruuth2002}.

  This contradiction completes the proof.
\end{proof}

\begin{corollary}\label{cor:no_SSP_5}
    Let $M$ denote an irreducible explicit Runge--Kutta method with $\mathcal{C}>0$.
    Then $M$ has effective order at most four.
\end{corollary}
\begin{proof}
	This follows immediately from Lemma~\ref{lem:positive_b} and
	Theorem~\ref{thm:effective_barrier}.
\end{proof}

\begin{remark}
    It is worth noting here an additional result that
    follows directly from what we have proved.
    Using Theorem~\ref{thm:effective_barrier} and \cite[Theorem~4.1]{dahlquist2006},
    it follows that any irreducible explicit Runge--Kutta method with positive radius of
    circle contractivity has effective order at most four.
\end{remark}

\section{Optimal explicit SSP Runge--Kutta schemes with maximal effective order}\label{sec:optimal_ESSPRK}
In this section, we use the SSP theory and Butcher's theory of effective
order (Sections \ref{sec:SSP} and \ref{sec:Algebraic_RK}) to find
optimal explicit SSP Runge--Kutta schemes with prescribed effective
order and classical order.
According to Corollary~\ref{cor:no_SSP_5}, there are no explicit SSPRK methods of
effective order five, and therefore we need only consider methods with
effective order up to four.

Recall from Section~\ref{sec:Algebraic_RK} that the methods with
an effective order of accuracy involve a main method $M$ as well as starting and
stopping methods $S$ and $S^{-1}$.
In Section~\ref{subsec:starting_stopping} we introduce a novel approach
to construction of starting and stopping methods in order to allow
them to be SSP.

We denote by ESSPRK($s,q,p$) an $s$-stage explicit SSP Runge--Kutta
method of effective order $q$ and classical order $p$.
Also we write SSPRK($s,q$) for an $s$-stage explicit SSP Runge--Kutta
method of order $q$.

\subsection{The main method}\label{subsec:main_method}

Our search is carried out in two
steps, first searching for optimal main methods $M$ and then for
possible corresponding methods $S$ and $S^{-1}$.
For a given number of stages, effective order, and classical order,
our aim is thus to find an optimal main method, meaning one with the
largest possible SSP coefficient $\mathcal{C}$.

To find a method ESSPRK($s,q,p$) with Butcher tableau $(A, \bm{b},
\bm{c})$, we consider the optimization problem~\eqref{eq:SSP_opt}
with $\Phi(K)$ representing the conditions for effective order
$q$ and classical order $p$ (as per Table~\ref{tab:effective_OCs}).
The methods are found through numerical search, using
\textsc{Matlab}'s optimization toolbox.
Specifically, we use \texttt{fmincon} with a sequential quadratic
programming approach \cite{Ketcheson2008, Ketcheson/Macdonald/Gottlieb:2009}.
This process does not guarantee a global minimizer, so many searches
from random initial guesses are performed to help
find methods with the largest possible SSP coefficients.

\subsubsection{Optimal SSP coefficients}\label{subsubsec:optimal_SSP_coeff}
Useful bounds on the optimal SSP coefficient can be obtained
by considering an important relaxation.
In the relaxed problem, the method is required to be accurate and strong
stability preserving only for linear, constant-coefficient initial value problems.
This leads to a reduced set of order conditions and a relaxed absolute
monotonicity condition \cite{Kraaijevanger1986,Ketcheson2008,ketcheson2009a}.
We denote the maximal SSP coefficient for linear problems
(maximized over all methods with order $q$ and $s$ stages) by $\sspcoef^{\textnormal{lin}}_{s,q}$.

Let $\mathcal{C}_{s,q}$ denote the maximal SSP coefficient (relevant to
non-linear problems) over all methods of $s$ stages with order $q$.  Let
$\mathcal{C}_{s,q,p}$ denote the object of our study, i.e. the maximal SSP
coefficient (relevant to non-linear problems) over all methods of $s$ stages
with effective order $q$ and classical order $p$.
From Remark~\ref{rem:talltrees} and the fact that the ESSPRK($s,q,p$) methods
form a super class of the SSPRK($s,q$) methods, we have
\begin{align} \label{ineq:clin}
        \mathcal{C}_{s,q} \le \mathcal{C}_{s,q,p} \le\sspcoef^{\textnormal{lin}}_{s,q}.
\end{align}

The effective SSP coefficients for methods with up to eleven stages are shown in
Table~\ref{tab:eff_SSP_coeff}.
Recall from Section~\ref{sec:ExRK_barrier} that $q=5$ implies a zero
SSP coefficient and from Section~\ref{sec:Algebraic_RK} that for
$q=1,2$, the class of explicit Runge--Kutta methods with effective order $q$
is the simply the class of explicit Runge--Kutta methods with order $q$.
Therefore we consider only methods of effective order $q=3$ and $q=4$.
Exact optimal values of $\sspcoef^{\textnormal{lin}}_{s,q}$ are known for many classes of methods; for
example see \cite{Kraaijevanger1986,Ketcheson2008,ketcheson2009a}.
Those results and \eqref{ineq:clin} allow us to determine the optimal value
of $\mathcal{C}_{s,q,p}$ {\em a priori} for the cases $q=3$ (for any $s$) and
for $q=4,s=10$, since in those cases we have $\mathcal{C}_{s,q}=\sspcoef^{\textnormal{lin}}_{s,q}$.

\begin{table}
\caption{Effective SSP coefficients $ \sspcoef_{\textnormal{eff}} = \mathcal{C}/s$ of the best known

      ESSPRK($s,q,p$) methods.
    		Entries in bold achieve the bound $\sspcoef^{\textnormal{lin}}_{s,q}$ given by the linear SSP coefficient and are therefore optimal.
    		If no positive $\mathcal{C}$ can be found, we use ``$-$'' to indicate non-existence.
    		The optimal fourth-order linear SSP coefficients are $\mathcal{C}^{\textnormal{lin}}_{4,4}=0.25$,
    		$\mathcal{C}^{\textnormal{lin}}_{5,4}=0.40$ and $\mathcal{C}^{\textnormal{lin}}_{6,4}=0.44$.}
    \centering
    \begin{tabular}{ccccccccccccc}
        \toprule

        \multirow{2}{*}{$q$} &
        \multirow{2}{*}{$p$}
            & \multicolumn{11}{c}{stages $s$} \\
            \cmidrule{3-13}
& & $1$ & $2$ & $3$ & $4$ & $5$ & $6$ & $7$ & $8$ & $9$ & $10$ & $11$ \\
        \midrule
        $3$ & $2$ & $-$ & $-$ & $\bf 0.33$ & $\bf 0.50$ & $\bf 0.53$ & $\bf 0.59$ & $\bf 0.61$ & $\bf 0.64$ & $\bf 0.67$ & $\bf 0.68$ & $\bf 0.69$\\

        $4$ & $2$ & $-$ & $-$ & $-$ & $0.22$ & $0.39$ & $\bf 0.44$ & $\bf 0.50$ & $\bf 0.54$ & $\bf 0.57$ & $\bf 0.60$ & $\bf 0.62$ \\

        $4$ & $3$ & $-$ & $-$ & $-$ & $0.19$ & $0.37$ & $0.43$ & $\bf 0.50$ & $\bf 0.54$ & $\bf 0.57$ & $\bf 0.60$ & $\bf 0.62$ \\
        \bottomrule
    \end{tabular}
    \label{tab:eff_SSP_coeff}
\end{table}

\subsubsection{Effective order three methods}\label{subsubsec:3rd_ESSPRK}
Since $\mathcal{C}_{s,q}=\sspcoef^{\textnormal{lin}}_{s,q}$ for $q=3$, the optimal effective order three methods
have SSP coefficients equal to the corresponding optimal classical order three methods.
In the cases of three and four stages, we are able to determine exact coefficients for
families of optimal methods of effective order three.
\begin{theorem}\label{thm:ESSPRK(3,3,2)}
	A family of optimal three-stage, effective order three SSP Runge--Kutta
	methods of classical order two, with SSP coefficient $\mathcal{C}_{3,3,2} = 1$, is given by
    \begin{displaymath}
    		\begin{split}
    			\bm{Y}_1 &= \bm{u}^n, \\
    			\bm{Y}_2 &= \bm{u}^n + \Delta t\bm{F}(\bm{Y}_1), \\
    			\bm{Y}_3 &= \bm{u}^n + \gamma\Delta t\bm{F}(\bm{Y}_1) + \gamma\Delta t\bm{F}(\bm{Y}_2), \\
    			\bm{u}^{n+1} &= \bm{u}^n + \frac{5\gamma-1}{6\gamma}\Delta t\bm{F}(\bm{Y}_1) + \frac{1}{6}\Delta t\bm{F}(\bm{Y}_2) + \frac{1}{6\gamma}\Delta t\bm{F}(\bm{Y}_3),
        \end{split}
    \end{displaymath}
    where $1/4 \leq \gamma \leq 1$ is a free parameter.

\end{theorem}
\begin{theorem}\label{thm:ESSPRK(4,3,2)}
	A family of optimal four-stage, effective order three SSP Runge--Kutta
	methods of classical order two, with SSP coefficient $\mathcal{C}_{4,3,2} = 2$ is given by
    \begin{displaymath}
    		\begin{split}
    			\bm{Y}_1 &= \bm{u}^n, \\
    			\bm{Y}_2 &= \bm{u}^n + \frac{1}{2}\Delta t\bm{F}(\bm{Y}_1), \\
    			\bm{Y}_3 &= \bm{u}^n + \frac{1}{2}\Delta t\bm{F}(\bm{Y}_1) + \frac{1}{2}\Delta t\bm{F}(\bm{Y}_2), \\
    			\bm{Y}_4 &= \bm{u}^n + \gamma\Delta t\bm{F}(\bm{Y}_1) + \gamma\Delta t\bm{F}(\bm{Y}_2) + + \gamma\Delta t\bm{F}(\bm{Y}_3), \\
    			\bm{u}^{n+1} &= \bm{u}^n + \frac{8\gamma-1}{12\gamma}\Delta t\bm{F}(\bm{Y}_1) + \frac{1}{6}\Delta t\bm{F}(\bm{Y}_2) + \frac{1}{6}\Delta t\bm{F}(\bm{Y}_3) + \frac{1}{12\gamma}\Delta t\bm{F}(\bm{Y}_4),
        \end{split}
    \end{displaymath}
    where $ 1/6 \leq \gamma \leq 1/2 $ is a free parameter.

\end{theorem}
\begin{proof}
	In either theorem, feasibility can be verified by direct calculation of the
	conditions in problem~\eqref{eq:SSP_opt}. Optimality follows because
	$\mathcal{C}_{s,3,2} = \mathcal{C}^{\textnormal{lin}}_{s,3}$.
\end{proof}

Theorem~\ref{thm:ESSPRK(3,3,2)} gives a \emph{family} of three-stage
methods.
The particular value of $\gamma = 1/4$ corresponds to the classical
Shu--Osher SSPRK($3,3$) method \cite{Gottlieb/Shu:1998}.
Similarly, in Theorem~\ref{thm:ESSPRK(4,3,2)} the particular value of
$\gamma = 1/6$ corresponds to the usual SSPRK($4,3$) method.
It seems possible that for each number of stages, the
ESSPRK($s, 3, 2$) methods may form a family in which an optimal
SSPRK($s$, $3$) method is a particular member.

\subsubsection{Effective order four methods}\label{subsubsec:4th_ESSPRK}
The ESSPRK($s,4,p$) methods can have classical order $p=2$ or $3$.
In either case, for stages $7 \le s \le 11$ the methods found are
optimal because the SSP coefficient attains the upper bound of
$\sspcoef^{\textnormal{lin}}_{s,q}$.
For fewer stages, the new methods still have SSP coefficients up to
30\% larger than that of explicit SSPRK($s,q$) methods.
In the particular case of four-stage methods we have the following
\begin{remark}
	In contrast with the non-existence of an SSPRK(4,\,4) method
	\cite{Gottlieb/Shu:1998,Ruuth2002},
	we are able to find ESSPRK(4,\,4,\,2) and ESSPRK(4,\,4,\,3) methods.
	The coefficients of these methods are found in
	Tables~\ref{tab:ESSPRK(4,4,2)_scheme}
	and~\ref{tab:ESSPRK(4,4,3)_scheme}.
\end{remark}

Additionally, we find two families of methods with effective order four,
for which $\sspcoef_{\textnormal{eff}}$ asymptotically approaches unity.
The families consist of second order methods with $s = n^2+1$ stages and
SSP coefficient $\mathcal{C}_{s,4,2} = n^2-n$.
They are optimal since $\mathcal{C}_{s,4,2}  = \mathcal{C}^{\textnormal{lin}}_{s,4}$
\cite[Theorem~5.2(c)]{Kraaijevanger1986}.
It is convenient to express the coefficients in the modified Shu--Osher form
\cite{Gottlieb2011a}
\begin{align*}
	\bm{Y}_i &= v_i\bm{u}^n + \sum_{j=1}^{i-1}\bigl(\alpha_{ij}\bm{Y}_j + \Delta t\beta_{ij}\bm{F}(\bm{Y}_j)\bigr), \; 1 \leq i \leq s+1 \\
	\bm{u}^{n+1} &= \bm{Y}_{s+1},
\end{align*}
because of the sparsity of the matrices $\alpha, \beta \in \mathbb{R}^{(s+1)\times s}$
and vector $\bm{v} \in \mathbb{R}^s$.
For $n \geq 3$ the non-zero elements are given by
\begin{align*}
	v_1 &= 1, \quad\quad v_{n^2+2} = \frac{2}{(n^2+1)\bigl((n-1)^2+1\bigr)}, \\
	\alpha_{n^2-2n+4,(n-2)^2} &= \frac{n^2-1 \pm \sqrt{n^3-3n^2+n+1}}{4n^2-6n+2}, \\
	\alpha_{n^2+2,n^2+1} &= \frac{n(n-1)^2}{(2n-1)(n^2+1)(1-\alpha_{n^2-2n+4,(n-2)^2})}, \\
	\alpha_{n^2+2,n^2-2n+2} &= 1 - v_{n^2+2} - \alpha_{n^2+2,n^2+1}, \\
	\alpha_{i+1,i} &= \begin{cases}
								1 - \alpha_{i+1,(n-2)^2}, & i = n^2-2n+3 \\
								1,  &\mbox{otherwise,}
							\end{cases}
\end{align*}
where $ 1 \leq i \leq n^2$ and
\begin{align*}
	\beta_{i,j} & = \frac{\alpha_{i,j}}{n^2-n}, \quad 1 \leq i \leq n^2+2, \;\; 1 \leq j \leq n^2+1.
\end{align*}
In \cite[\S~6.2.2]{Gottlieb2011a}, a similar pattern was found
for SSPRK($s,3$) methods.

\subsection{Starting and stopping methods}\label{subsec:starting_stopping}
Provided an ESSPRK($s,q,p$) scheme that can be used as the main
method $M$, we want to find perturbation methods $S$ and $S^{-1}$ such that the
Runge--Kutta scheme $S^{-1}MS$ attains classical order $q$, equal to the
effective order of method $M$.
We also want the resulting overall process to be SSP.
However at least one of the $S$ and $S^{-1}$ methods is not SSP:
if $\beta_1 = 0$ then $\sum_i b_i = 0$ implies the presence of at
least one negative weight and thus neither scheme can be SSP.
Even if we consider methods with $\beta_1 \neq 0$, one of $S$ or
$S^{-1}$ must step backwards and thus that method cannot be SSP
(unless we consider the downwind operator
\cite{Ruuth2004,Gottlieb/Ruuth:SSPfastdownwind,Ketcheson:2011:downwind}).

In order to overcome this problem and achieve ``bona fide''
SSPRK methods with an effective order of accuracy, we need to choose different starting and stopping methods.
We consider methods $R$ and $T$ which each take a positive step such that
$R \equival{q} MS$ and $T \equival{q} S^{-1}M$.
That is, the order conditions of $R$ and $T$ must match those of
$MS$ and $S^{-1}M$, respectively, up to order $q$.
This gives a new $TM^{n-2}R$ scheme which is equivalent up to order $q$
to the $S^{-1}M^nS$ scheme and attains classical order $q$.
Each starting and stopping procedure now takes a positive step forward
in time.

To derive order conditions for the $R$ and $T$ methods, consider their
corresponding functions in group $G$ to be $\rho$ and $\tau$
respectively.
Then the equivalence is expressed as
\begin{equation} \label{eq:R_T_OCs}
    \rho(t) = (\beta\alpha)(t) \text{ and } \tau(t) = (\alpha\beta^{-1})(t), \quad \text{for all
    trees $t$ with $r(t) \leq q$.}
\end{equation}
Rewriting the second condition in \eqref{eq:R_T_OCs} as
$(\tau\beta)(t) = \alpha(t)$, the order conditions for the starting and stopping
methods can be determined by the usual product formula and are given in Table~\ref{tab:rho_tau_OCs}.
These conditions could be constructed more generally but here we have
assumed $\beta_1=0$ (see Section~\ref{subsubsec:Main_starting_conditions}); this
will be sufficient for constructing SSP starting and stopping
conditions.

\begin{table}
  	\caption{Order conditions on $\rho$ and $\tau$ up to effective order four for starting
  		and stopping methods $R$ and $T$, respectively.
  		The upper block represents the effective order three conditions.
         As in Table~\ref{tab:effective_OCs_on_alpha} and Table~\ref{tab:effective_OCs} we assume
         $\beta_1 = 0$.}
	\centering
	\begin{tabular}{lcl}
		\toprule

    		$\rho(t) = (\beta\alpha)(t)$ & & $\tau(t) = (\alpha\beta^{-1})(t)$ \\
    		\midrule
    		 $\rho_1 = \alpha_1$ & & $\tau_1 = \alpha_1$ \\
    		$\rho_2 = \alpha_2 + \beta_2$ & & $\tau_2 = \alpha_2 - \beta_2$ \\
    		$\rho_3 = \alpha_3 + \beta_3$ & & $\tau_3 = \alpha_3 - 2\alpha_1\beta_2 - \beta_3$ \\
    		$\rho_4 = \alpha_4 + \alpha_1\beta_2 + \beta_4$ & & $\tau_4 = \alpha_4 - \alpha_1\beta_2 - \beta_4$ \\
                \mydashrule
		$\rho_5 = \alpha_5 + \beta_5$ & & $\tau_5 = \alpha_5 - 3\alpha_1^2\beta_2 - 3\alpha_1\beta_3 - \beta_5$ \\
		$\rho_6 = \alpha_6 + \alpha_2\beta_2 + \beta_6$ & & $\tau_6 = \alpha_6 - (\alpha_1^2 + \alpha_2 -\beta_2)\beta_2 -\alpha_1\beta_3 - \alpha_1\beta_4 - \beta_6$ \\
		$\rho_7 = \alpha_7 + \alpha_1\beta_3 + \beta_7$ & & $\tau_7 = \alpha_7 - 2\alpha_1\beta_4 - \alpha_1^2\beta_2 - \beta_7$ \\
		$\rho_8 = \alpha_8 + \alpha_1\beta_4 + \alpha_2\beta_2 + \beta_8$ & & $\tau_8 = \alpha_8 - \alpha_1\beta_4 - \alpha_2\beta_2 + \beta_2^2 -  \beta_8$
                \\
                \bottomrule
  	\end{tabular}
  	\label{tab:rho_tau_OCs}
\end{table}

\subsubsection{Optimizing the starting and stopping methods}\label{subsubsec:opt_methods}
It turns out that the order conditions from \eqref{eq:R_T_OCs} do not
contradict the SSP requirements.
We can thus find methods $R$ and $T$ using the optimization procedure
described in Section~\ref{subsec:Optimal_SSPRK} with the order conditions
given by Table~\ref{tab:rho_tau_OCs} for $\Phi(K)$ in \eqref{eq:SSP_opt}.

The values of $\alpha_i$ are determined by the main method $M$.
Also note that for effective order $q$, the algebraic expressions on
$\beta$ up to order $q-1$ are already found by the optimization procedure of
the main method (see Table~\ref{tab:effective_OCs}).
However, the values of the order $q$ elementary weights on $\beta$ are not
known; these are $\beta_3$ and $\beta_4$ for effective order three and
$\beta_5$, $\beta_6$, $\beta_7$ and $\beta_8$ for effective order four.
From Table~\ref{tab:rho_tau_OCs}, we see that both the $R$ and $T$
methods depend on these parameters.
Our approach is to optimize for both methods at once: we solve a
modified version of the optimization problem \eqref{eq:SSP_opt} where
we simultaneously maximize both SSP coefficients subject to the
constraints given in \eqref{eq:R_T_OCs} and conditions on $\beta$ given by
Table~\ref{tab:effective_OCs}.
The unknown elementary weights on $\beta$ are used as free parameters.
In practice, we maximize the objective function $\min(r_1,r_2)$, where $r_1$
and $r_2$ are the radii of absolute monotonicity of the methods $R$ and $T$.

We were able to construct starting and stopping schemes for each main
method, with an SSP coefficient at least as large as that of the main method.
This allows the usage of a uniform time-step $\Delta t \leq \mathcal{C}\Delta t_{\textnormal{FE}}$,
where $\mathcal{C}$ is the SSP coefficient of the main method.
The additional computational cost
of the starting and stopping methods is minimal:
for methods $R$ and $T$ associated with an $s$-stage main method,
at most $s + 1$ and $s$ stages, respectively, appear to be required.
Tables \ref{tab:ESSPRK(4,4,2)_scheme} and \ref{tab:ESSPRK(4,4,3)_scheme}
show the coefficients of the schemes where the main
method is ESSPRK($4,4,2$) and ESSPRK($4,4,3$), respectively.

It is important to note that in practice, if accurate values are needed at
any time other than the final time, the computation must invoke the
stopping method to obtain them.  Furthermore, changing step-size would require first applying the stopping method with the old step-size and then applying the starting method with the new step-size.

\begin{table}
    \caption{ESSPRK(4,4,2): an effective order four SSPRK method with
      four stages and classical order two with its associated starting
      and stopping methods.}
    \setlength{\tabcolsep}{2pt}

    \footnotesize
    \subfloat[Main method $M$, ESSPRK($4,4,2$) \label{ESSPRK(4,4,2)_scheme_a}]{
        \begin{tabular}{c | c c c c}
             $0$ & & & & \\
             $0.730429885783319$ & $0.730429885783319$ & & & \\
             $0.644964638145795$ & $0.251830917810810$ & $0.393133720334985$ & & \\
             $1.000000000000000$ & $0.141062771617064$ & $0.220213358584678$ & $0.638723869798257$ & \\
             \hline
             & $0.384422161080494$ & $0.261154113377550$ & $0.127250689937518$ & $0.227173035604438$
        \end{tabular}
    }\\
    \subfloat[Starting method $R$ \label{ESSPRK(4,4,2)_scheme_b}]{
        \begin{tabular}{c | c c c c c}
			$0$ & & & & & \\
			$0.545722177514735$ & $0.545722177514735$ & & & & \\
			$0.842931687441527$ & $0.366499989048164$ & $0.476431698393363$ & & & \\
			$0.574760809487828$ & $0.135697968350722$ & $0.176400587890242$ & $0.262662253246864$ & & \\
			$0.980872743236632$ & $0.103648417776838$ & $0.134737771331049$ & $0.200625899485633$ & $0.541860654643112$ & \\
             \hline
             & $0.233699169638954$ & $0.294263351266422$ & $0.065226988215286$ & $0.176168374199685$ & $0.230642116679654$
        \end{tabular}
    }\\
    \subfloat[Stopping method $T$ \label{ESSPRK(4,4,2)_scheme_c}]{
        \begin{tabular}{c | c c c c}
			$0$ & & & & \\
			$0.509877496215340$ & $0.509877496215340$ & & & \\
			$0.435774135529007$ & $0.182230305923759$ & $0.253543829605247$ & & \\
			$0.933203341300203$ & $0.148498121305090$ & $0.206610981494095$ & $0.578094238501017$ & \\
             \hline
             & $0.307865440399752$ & $0.171863794704750$ & $0.233603236964822$ & $0.286667527930676$
        \end{tabular}
    }
    \label{tab:ESSPRK(4,4,2)_scheme}
\end{table}

\begin{table}
    \caption{ESSPRK(4,4,3): an effective order four SSPRK method with
      four stages and classical order three with its associated starting
      and stopping methods.}
    \setlength{\tabcolsep}{2pt}

    \footnotesize
    \subfloat[Main method $M$, ESSPRK($4,4,3$) \label{ESSPRK(4,4,3)_scheme_a}]{
        \begin{tabular}{c | c c c c}
             $0$ & & & & \\
             $0.601245068769724$ & $0.601245068769724$ & & & \\
             $0.436888719886063$ & $0.139346829159954$ & $0.297541890726109$ & & \\
             $0.747760163757110$ & $0.060555450075478$ & $0.129301708677891$ & $0.557903005003740$ & \\
             \hline
             & $0.220532078662434$ & $0.180572397883936$ & $0.181420582644840$ & $0.417474940808790$
        \end{tabular}
    }\\
    \subfloat[Starting method $R$ \label{ESSPRK(4,4,3)_scheme_b}]{
        \begin{tabular}{c | c c c c c}
			$0$ &  & & & & \\
			$0.438463764036947$ & $0.438463764036947$ & & & & \\
			$0.639336395725557$ & $0.213665532574654$ & $0.425670863150903$ & & & \\
			$0.434353425654020$ &$0.061345094040860$ & $0.122213530726218$ & $0.250794800886942$ & & \\
			$0.843416464962307$ & $0.039559973266996$ & $0.078812561688700$ & $0.161731525131914$ & $0.563312404874697$ & \\
             \hline
             & $0.154373542967849$ & $0.307547588471376$ & $0.054439037790856$ & $0.189611674483496$ & $0.294028156286422$
        \end{tabular}
    }\\
    \subfloat[Stopping method $T$ \label{ESSPRK(4,4,3)_scheme_c}]{
        \begin{tabular}{c | c c c c}
			$0$ & & & & \\
			$0.556337718891090$ & $0.556337718891090$ & & & \\
			$0.428870688216872$ & $0.166867537553458$ & $0.262003150663414$ & & \\
			$0.815008947642716$ & $0.104422177204659$ & $0.163956032598547$ & $0.546630737839510$ & \\
             \hline
             & $0.203508169408374$ & $0.096469758967330$ & $0.321630956102914$ & $0.378391115521382$
        \end{tabular}
    }
    \label{tab:ESSPRK(4,4,3)_scheme}
\end{table}

\section{Numerical experiments}\label{sec:numerics}
Having constructed strong stability preserving $TM^{n-2}R$ schemes in the previous
section, we now numerically verify their properties.
Specifically, we use a convergence study to show that the procedure
attains order of accuracy $q$, the effective order of $M$.
We also demonstrate on Burgers' equation that the
SSP coefficient accurately measures the maximal time-step for which the
methods are strong stability preserving.

\subsection{Convergence study}\label{subsec:convergence}
We consider the van der Pol system \cite{Hairer1987_book}
\begin{equation}\label{eq:conv_eq}
	\begin{split}
    		u_1'(t) &= u_2(t), \\
                u_2'(t) &= \mu \bigl(1 - u_1^2(t)\bigr)u_2(t) - u_1(t),
    \end{split}
\end{equation}
over the time interval $t \in [0, 50]$ with $\mu = 2$ and initial
values $u_1(0) = 2$ and $u_2(0) = 1$.
The reference solution for the convergence study is calculated by \textsc{Matlab}'s
\texttt{ode45} solver with relative and absolute tolerances set to $10^{-13}$.

We solve the initial value problem \eqref{eq:conv_eq}
using SSP $TM^{n-2}R$ schemes.
The solution is computed using $n = 100 \cdot 2^{k}$ time steps for
$k = 2, \dots, 7$.
The error at $t = 50$ with respect to time-step is shown in
Figure~\ref{fig:conv_study} on a logarithmic scale.
\begin{figure}
	\centering
     \subfloat[ESSPRK($s,3,2$)]{\label{fig:conv_study_a}
     \includegraphics[width=0.45\textwidth]{convergence_3_order}}
   \qquad
     \subfloat[ESSPRK($s,4,2$)]{\label{fig:conv_study_b}
    \includegraphics[width=0.45\textwidth]{convergence_4_order}}
    \caption{Convergence study for the van der Pol system using $TM^{n-2}R$
    Runge--Kutta schemes when (a) $M$
    is an ESSPRK($s,3,2$) method and (b) $ M $ is an ESSPRK($s,4,2$) method.}
    \label{fig:conv_study}
\end{figure}
The convergence study is performed for $TM^{n-2}R$ schemes with
various number of stages $s$ and the results show that the schemes
attain an order of accuracy equal to the effective order of their main
method $M$.
It is important in doing this sort of convergence study that the
effective order of accuracy can only be obtained after the stopping method is
applied.
Intermediate steps will typically only be order $p$ accurate (the classical
order of the main method).
Finally, we note that the methods with more stages generally exhibit smaller errors
(for a given step size).

\subsection{Burgers' equation}\label{subsubsec:burgers}
The inviscid Burgers' equation consists of the scalar hyperbolic conservation law
\begin{align}\label{eq:HCL}
    U_{t} + f(U)_{x} = 0,
\end{align}
with flux function $f(U) = \frac{1}{2}U^{2}$.
We consider initial data
$U(0,x)  = \frac{1}{2} - \frac{1}{4}\sin{\pi x}$,
on a periodic domain $x \in [0,2)$.
The solution advances to the right where it eventually exhibits a shock.
We perform a semi-discretization
using an upwind approximation to obtain the system of ODEs
\begin{align*}\label{eq:burgers_flux}
	\frac{\textrm{d}}{\textrm{d} t} u_i = -\frac{f(u_{i}) - f(u_{i-1})}{\Delta x}.
\end{align*}
This spatial discretization is total-variation-diminishing (TVD) when
coupled with the forward Euler method under the restriction
\cite{Laney:1998}
$$\Delta t \leq {\Delta t}_{\text{FE}} = \Delta x / \|U(0,x)\|_{\infty}.$$
Recall that a time discretization with SSP
coefficient $\mathcal{C}$ will give a TVD solution for $\Delta t \leq
\mathcal{C} {\Delta t}_{\text{FE}}$.

Burgers' equation was solved using an SSP $TM^{n-2}R$ scheme with time-step
restriction $\Delta t = \sigma{\Delta t}_{\text{FE}}$, where $\sigma$ indicates the size
of the time step.
We integrate to roughly time $t_\text{f} = 1.62$ with $200$ points in space.
Figure~\ref{fig:burgers_cont} shows that if $\sigma$ is chosen less than the SSP
coefficient of the main method, then no oscillations are observed.
If this stability limit is violated, then oscillations may appear,
as shown in Figure~\ref{fig:burgers_cont_b}.
We measure these oscillations by computing the total variation of the
numerical solution.

\begin{figure}
    \centering
    \subfloat[$\sigma = 0.88$]{\label{fig:burgers_cont_a}
      \includegraphics[width=0.45\textwidth]{burgers_cont}}
    \qquad
    \subfloat[$\sigma = 1.60$]{\label{fig:burgers_cont_b}
      \includegraphics[width=0.45\textwidth]{burgers_cont_no_tvd}}
    \caption{Solution of Burgers' equation at the final time with continuous initial data, using a
    $TM^{n-2}R$ scheme, where $M$ is the optimal ESSPRK($4,4,2$).
    The time-step used is $\Delta t = \sigma{\Delta t}_{\text{FE}}$:
    in (a) we take $\sigma$ as the SSP coefficient
    $\sigma = \mathcal{C} = s \times 0.22 = 0.88$
    (see Table~\ref{tab:eff_SSP_coeff}) and no oscillations are observed.
    However, in (b) we take $\sigma > \mathcal{C}$ and we observe oscillations
    (note (b) is magnified to show these).

    Here $TV$ denotes the discrete total variation of the solution at the final time:
    a value greater than 1 (the total variation of the initial condition)
    indicates a violation of the TVD condition.

    }
    \label{fig:burgers_cont}
\end{figure}

We also consider Burgers' equation with a discontinuous
square wave initial condition
\begin{equation}\label{eq:burgers_discont_IC}
    U(0,x)  = \left\{
                \begin{array}{ll}
                  1, & \hbox{$0.5 \leq x \leq 1.5$} \\
                  0, & \hbox{otherwise.}
                \end{array}
              \right.
\end{equation}
The solution consists of a rarefaction (i.e., an expansion fan) and a
moving shock.
Again we use $200$ points in space and we compute the solution until
roughly time $t_\text{f} = 0.6$, using a time-step $\Delta t = \sigma{\Delta t}_{\text{FE}}$.
Figure~\ref{fig:burgers_discont} shows the result of solving the
discontinuous problem using an SSP $TM^{n-2}R$ scheme, where $M$ is an
ESSPRK($5,4,2$) method with SSP coefficient $\mathcal{C} = 1.95$.
In this case, $\sigma = 1.98$ appears to be the largest value
for which the total variation is monotonically decreasing during the
calculation.
This is only $2\%$ larger than the value of the SSP coefficient.
Figure~\ref{fig:burgers_discont_b} shows part of the solution exhibiting
oscillations when $\sigma$ is larger than the SSP coefficient.

\begin{figure}
    \centering
    \subfloat[$\sigma = 1.95$]{\label{fig:burgers_discont_a}
      \includegraphics[width=0.45\textwidth]{burgers_discont}}
    \qquad
    \subfloat[$\sigma = 2.15$]{\label{fig:burgers_discont_b}
      \includegraphics[width=0.45\textwidth]{burgers_discont_no_tvd}}
    \caption{Solution of Burgers' equation at the final time with discontinuous initial data, using a
    $TM^{n-2}R$ scheme, where $M$ is ESSPRK($5,4,2$) method.
    The time-step used is $\Delta t = \sigma{\Delta t}_{\text{FE}}$:
    in (a) we take $\sigma$ as the SSP coefficient
    $\sigma = \mathcal{C} = s \times 0.39 = 1.95$

    and no oscillations are observed.
    However, in (b) we take $\sigma > \mathcal{C}$ and, when magnified,
    we observe oscillations.

    Here $TV$ denotes the total variation semi-norm of the solution at the final time:
    a value greater than 2 indicates an increase in total variation.

    }
    \label{fig:burgers_discont}
\end{figure}

For various schemes, Table~\ref{tab:observed_SSP_coeff} shows the
maximum observed values of $\sigma$ for which the numerical solution
is total variation decreasing for the entire computation.  With the
exception of the four-stage effective order four methods, we note good
agreement between these experimental values and the SSP coefficients
predicted by the theory.

\begin{table}
    \caption{Maximum observed coefficients exhibiting the TVD property on
    the Burgers' equation example with discontinuous data \eqref{eq:burgers_discont_IC}.
    The numbers in parenthesis indicate the increase relative to the corresponding SSP
    coefficients.}
    \centering
    \begin{tabular}{>{\hspace*{-4pt}}c@{\hspace{3pt}}c@{\hspace{3pt}}c@{\hspace{2pt}} c@{\hspace{2pt}} c@{\hspace{2pt}} c@{\hspace{2pt}} c@{\hspace{2pt}} c@{\hspace{2pt}} c@{\hspace{2pt}} c@{\hspace{2pt}} c@{\hspace{2pt}}}
        \toprule

        \multirow{2}{*}{$q$} &
        \multirow{2}{*}{\;$p$\;}
               &   \multicolumn{9}{c}{stages $s$} \\
            \cmidrule{3-11}
        &      &   $3$ & $4$ & $5$ & $6$ & $7$ & $8$ & $9$ & $10$ & $11$ \\
        \midrule
        $3$ & $2$ & \small$1.04(4\%)$ & \small$2.00(0\%)$ & \small$2.65(0\%)$ & \small$3.52(0\%)$ & \small$4.29 (0\%)$ & \small$5.11(0\%)$ & \small$6.00(0\%)$ & \small$6.79(0\%)$ & \small$7.63(0\%)$\\
        $4$ & $2$ & \small$-$ & \small$1.07(22\%)$ & \small$1.98(2\%)$ & \small$2.69(2\%)$ & \small$3.56(1\%)$ & \small$4.33(1\%)$ & \small$5.16(1\%)$ & \small$6.05(1\%)$ & \small$6.84(1\%)$ \\
        $4$ & $3$ & \small$-$ & \small$1.05(35\%)$ & \small$1.89(3\%)$ & \small$2.63(2\%)$ & \small$3.53(1\%)$ & \small$4.31(1\%)$ & \small$5.16(1\%)$ & \small$6.04(1\%)$ & \small$6.85(1\%)$ \\
        \bottomrule
    \end{tabular}
    \label{tab:observed_SSP_coeff}
\end{table}

We also note the necessity of our modified starting and stopping methods in the $RM^{n-2}T$ approach: in this
example if we use the original approach of $S$ and $S^{-1}$,
the solution exhibits oscillations immediately following the
application of the starting perturbation method $S$.

\section{Conclusions}\label{sec:Conclusion}
We use the theory of strong stability preserving time discretizations
with Butcher's algebraic interpretation of order to construct
explicit SSP Runge--Kutta schemes with an effective order of accuracy.
These methods, when accompanied by starting and stopping
methods, attain an order of accuracy higher than their (classical) order.
We propose a new choice of starting and stopping methods to allow the
overall procedure to be SSP.
We prove that explicit Runge--Kutta methods with strictly positive
weights have at most effective order four.
This extends the barrier already known in the case of classical order
explicit SSPRK methods.

SSP Runge--Kutta methods of effective order three and four
are constructed by numerical optimization.
Most of the methods found are optimal because they achieve
the upper bound on the SSP coefficient known from linear
problems.
Also, despite the non-existence of four-stage, order four explicit SSPRK methods,
we find effective order four methods with four stages (of classical
order two and three).
We perform numerical tests which confirm the accuracy and
SSP properties of the new methods.

The ideas here are applied to explicit Runge--Kutta methods, but they
could also be applied to other classes of methods including implicit
Runge--Kutta methods, general linear methods, and Rosenbrock methods.

\section*{Acknowledgments}{
                The authors would like to thank the anonymous referees for their helpful
                and valuable suggestions on the paper.
        }

        \bibliographystyle{acm}