RK-Opt: A package for the design of numerical ODE solvers

Ordinary and partial differential equations (ODEs and PDEs) are used to model many important phenomena. In most cases, solutions of these models must be approximated by numerical methods. Most of the relevant algorithms fall within a few classes of methods, with the properties of individual methods determined by their coefficients. The choice of appropriate coefficients in the design of methods for specific applications is an important area of research. RK-Opt is a software package for designing numerical ODE solvers with coefficients optimally chosen to provide desired properties. It is available from https://github.com/ketch/RK-Opt, with documentation at http://numerics.kaust.edu.sa/RK-Opt/. The primary focus of the package is on the design of Runge-Kutta methods, but some routines for designing other classes of methods such as multistep Runge-Kutta and general linear methods are also included.


Summary
Ordinary and partial differential equations (ODEs and PDEs) are used to model many important phenomena. In most cases, solutions of these models must be approximated by numerical methods. Most of the relevant algorithms fall within a few classes of methods, with the properties of individual methods determined by their coefficients. The choice of appropriate coefficients in the design of methods for specific applications is an important area of research. RK-Opt is a software package for designing numerical ODE solvers with coefficients optimally chosen to provide desired properties. It is available from https://github.com/ketch/RK-Opt, with documentation at http://numerics.kaust.edu.sa/RK-Opt/. The primary focus of the package is on the design of Runge-Kutta methods, but some routines for designing other classes of methods such as multistep Runge-Kutta and general linear methods are also included.

Statement of need
Over the last several decades, a great deal of work has gone into the design of numerical ODE solvers. Initially this work was aimed at developing general purpose solvers, but over time the emphasis shifted increasingly toward development of optimized methods for specific applications. Different accuracy, stability, performance, and other properties may be relevant or essential depending on the nature of the equations to be solved.
An s-stage Runge-Kutta method has roughly s^2 coefficients (roughly s^2/2 for explicit methods), which can be chosen so as to provide high accuracy, stability, or other properties. Historically, most interest in Runge-Kutta methods has focused on methods using the minimum number of stages for a given order of accuracy. However, in the past few decades there has been increasing recognition that using extra stages can be worthwhile in order to improve other method properties. Some areas where this is particularly useful are in the enhancement of linear and nonlinear stability properties, the reduction of storage requirements, and the design of embedded pairs. Methods with dozens or even hundreds of stages are not unheard of.
At the same time, most existing Runge-Kutta methods have been designed by hand, by researchers laboriously solving the order conditions. When using extra stages, the number of available parameters makes the selection of a near-optimal choice by hand impossible, and one resorts to computational optimization. This leads to a different paradigm of numerical method design, in which we use sophisticated numerical (optimization) algorithms to design sophisticated numerical (integration) algorithms. It can be expected that this trend will accelerate in the future, and perhaps one day simple manually-constructed algorithms will be the exception.
RK-Opt contains a set of tools for designing Runge-Kutta methods in this paradigm. It provides code that can enforce desired properties and/or objective functions. The constraints and objective are then used within an optimization framework, to determine coefficients of methods that best achieve the desired goal. Thus, RK-Opt is a sort of meta-software, consisting of algorithms whose purpose is to create other algorithms.
Typically, the most obvious formulation of the corresponding optimization problem is intractable. Therefore, these problems are reformulated in ways that make them amenable to available techniques. These reformulations include, for instance, turning a nonconvex problem into a sequence of convex problems or even linear programs. The resulting algorithms can often guarantee optimality of their output. However, for the general problem of determining Runge-Kutta coefficients, the nonconvex problem must be attacked directly and optimality cannot be guaranteed.
RK-Opt is written entirely in MATLAB, and leverages the MATLAB Optimization Toolbox as well as the Global Optimization Toolbox. Its development has been motivated largely by research needs and it has been used in a number of papers (see below).

RK-Opt includes the following subpackages.
polyopt This package computes optimal stability functions for Runge-Kutta methods. Here optimal means that the stable step size is maximized for a given ODE spectrum. The corresponding optimization problem is intractable under a direct implementation. The package uses the algorithm developed in (Ketcheson & Ahmadia, 2012), which relaxes the global optimization problem by solving a sequence of convex subproblems. Under certain technical assumptions, the result is guaranteed to be the optimal solution of the original problem. polyopt relies on CVX (Grant & Boyd, 2008 to solve the convex subproblems. This package is usually used as the first step in designing a Runge-Kutta method.

RK-Coeff-Opt
This package computes optimal Runge-Kutta coefficients based on a desired set of constraints and an objective. Available constraints include: • The number of stages and order of accuracy • The class of method (explicit, implicit, diagonally implicit, low-storage) • The coefficients of the stability polynomial (usually determined using polyopt) Two objective functions are provided; methods can be optimized for the strong stability preserving (SSP) coefficient or the principal error norm (a measure of the leading-order truncation error coefficients). In addition to standard Runge-Kutta methods, various classes of multistep Runge-Kutta methods can also be optimized.
The optimization problem in question is highly nonconvex and the available solvers may fail to find a solution, or may converge to a non-optimal solution. For this reason, the implementation is based on solving many local optimization problems in parallel from different random initial points, using MATLAB's Global Optimization Toolbox.
The packages dwrk-opt and low-storage are specialized but less full-featured versions of RK-Coeff-Opt that were developed for specific research projects involving downwind Runge-Kutta methods and low-storage Runge-Kutta methods, respectively.

am_radius-opt
Whereas the previous two subpackages are fairly general-purpose tools, this package solves a very specific and discrete set of problems described in (Ketcheson, 2009). Specifically, the provided routines determine the coefficients of multistep methods (including classes of general linear methods) with the largest possible SSP coefficient (also known as radius of absolute monotonicity). The corresponding optimization problem had previously been attacked using brute force search, but this limited its solvability to methods with very few steps. In this package the problem is reformulated as a sequence of linear programming problems, enabling its efficient solution for methods with many steps.

Related research and software
RK-Opt development has proceeded in close connection to the NodePy package (https:// github.com/ketch/NodePy). Whereas RK-Opt is focused on the design of numerical methods, NodePy is focused more on their analysis. A common workflow involves generating new methods with RK-Opt and then studying their properties in more detail using NodePy.
Because of the nature of RK-Opt, applications often involve writing some additional code to impose special constraints, or simply using the existing code as a template. A number of related optimization routines written for similar purposes in this vein can be found at https://github.com/SSPmethods.