A parallel-in-time approach for wave-type PDEs
Abe C. Ellison Bengt Fornberg
math.NA cs.NA 65L06, 65L05, 65M20, 65Y05







    Department of Applied Mathematics 

    University of Colorado 

    Boulder, CO 80309 

    abe.ellison@colorado.edu






A parallel-in-time approach for wave-type PDEs
    Abe C. Ellison Bengt Fornberg
    Received: date / Accepted: date
==============================================




Numerical solutions to wave-type PDEs utilizing method-of-lines require
the ODE solver's stability domain to include a large stretch of the
imaginary axis surrounding the origin.  We show here that extrapolation 
based solvers of Gragg-Bulirsch-Stoer (GBS) type can meet this requirement.
Extrapolation methods utilize several independent time stepping sequences,
making them highly suited for parallel execution.  Traditional extrapolation schemes 
use all time stepping sequences to maximize the method's order of accuracy.
The present method instead maintains a desired order of accuracy while employing
additional time stepping sequences to shape the resulting stability domain.
We optimize the extrapolation coefficients to maximize the stability domain's imaginary axis coverage.
This yields a family of explicit schemes that approaches maximal time step size for wave propagation problems.
On a computer with several cores we achieve both high order and fast time to solution compared
with traditional ODE integrators.


65L06 
65L05 
65M20 
65Y05 






§ INTRODUCTION


Time integration of ODEs is an inherently sequential process, since
each forward step ought to be based on the most recent information
available. Three conceivable options for achieving some level of parallel-in-time
are (i) to have correction calculations follow the explicit forward
steps as closely behind as possible, letting them catch up frequently,
(ii) to carry out `preparatory' calculations that are based on trying
to anticipate later solution states, and (iii) to exploit extrapolation
ideas. While all of these concepts have been pursued for systems of
ODEs, as summarized in <cit.>, their performance is unclear
for ODE systems arisen from method-of-lines (MOL) discretization of
wave-type PDEs. The additional requirement that arises
then is that the ODE solver's stability domain must include
a quite large stretch of the imaginary axis surrounding the origin.
We show here that extrapolation-based ODE solvers of Gragg-Bulirsch-Stoer
(GBS) type can meet this requirement. In particular, one such scheme
that we will focus on steps forward explicitly using six cores as
fast as Forward Euler (FE) does on one core, but combines eighth order
of accuracy with a generously sized stability domain. In contrast
to linear multistep methods, it needs no back levels in time to get
started. The present approach is compared against explicit Runge-Kutta (RK) methods
for a PDE test problem. 

Standard Richardson extrapolation schemes utilize a square Vandermonde-type system
to compute the extrapolation weights.  This system is constructed
to cancel successive terms in the asymptotic error expansion of the time stepper.
By allowing more columns than rows in the system - that is, more extrapolation
components than order constraints - we create an underdetermined system
that grants degrees of freedom to optimize the extrapolated stability domain.
For wave-type PDEs stepped with method-of-lines we optimize the stability
domain along the imaginary axis.  We achieve stability domains far larger
than those of both the square extrapolation systems and other standard ODE integrators,
thereby enabling large time step sizes and thus faster time to solution.



§ GBS-TYPE ODE SOLVERS




 §.§ GBS concept


We consider first the problem of advancing forward in time an ODE of
the form y'(t)=f(t,y), where the unknown function y(t) is either
scalar or vector valued. The complete time interval of interest is
split into N sections. For each of these sections, the
basic (unextrapolated) GBS scheme consists of the steps

    {[       y_1-y_0/h=f(t_0,y_0)         Forward Euler (FE); y_n+1-y_n-1/2h_=f(t_n,y_n)  Leap-frog (LF), n=1,2,…,N; y_N^*=14(y_N-1+2y_N+y_N+1)                 Averaging, ].

after which y_N^* is accepted as the new value at time t_N.
The initial FE step is accurate to first order, while the subsequent
LF steps are second order accurate. One would therefore expect y_N^*
to be accurate to at most second order, and have an error expansion
in which all further powers of h would be present. Remarkably,
for any smooth (linear or nonlinear) function f(t,y), it transpires
that, if N is even, all odd powers in the expansion will vanish
<cit.>:

    Error=y_N^*-y(t_N)=a_2h^2+a_4h^4+a_6h^6+…

The form of the expansion makes Richardson extrapolation
particularly efficient since, each time this is applied, the result
will gain two orders of accuracy. For example, the results from four
completely independent calculations over the same section in time,
using different h-values, can be combined to give an 𝒪(h^8)-accurate
result. These four calculations require no communications between
each other, and can therefore be run simultaneously on separate cores.

Even when not counting the cost of the work on the extra cores, the GBS approach does not
offer any striking benefits for standard ODE systems, unless possibly
if extrapolated to very high orders. However, in the present context
of wave-type PDEs, the situation becomes different, since GBS methods
can be designed to feature particularly favorable stability domains. 



 §.§ Stability domains for GBS-type methods


<ref> briefly summarizes the definition of an ODE solver's
stability domain, explains its significance in the context of
MOL time stepping, and provides stability domain information for some
well-known explicit ODE solvers. These domains should be contrasted
to the corresponding ones for GBS methods described below. The imaginary
stability boundary (ISB) of an ODE solver is defined as the largest
value such that the imaginary axis is included from -i·ISB
to +i·ISB. For solvers that lack any imaginary axis coverage,
we define their ISB to be zero. 

In order to provide a fair comparison between different methods, we
will from now on further normalize all ISB values by the number of
function evaluations that each step requires and denote this ISB_n.
For example, we divide
RK_4's ISB, stated in <ref> as 2.8284,
by four to compensate for its four stages (with one function evaluation
in each), i.e. we list its ISB_n
as 0.7071 (=1/√(2)). Similarly, the ISB_n
for the 13-stage RK_8 method becomes 0.2848.
With this normalization, the largest feasible ISB_n
for any explicit method becomes one <cit.>, which is realized
by the LF scheme. Since the longest distance a solution can be advanced
per function evaluation is proportional to the time stepping method's
ISB_n, a key goal will be to design a method that has both
high order and a large ISB_n. 

Stability domains and ISBs for GBS-type schemes do not appear to have been
studied until in <cit.>. One key observation made there was that GBS schemes of
orders 4, 8, 12, ... will feature positive ISBs, whereas schemes of order 6,
10, 14, ... will not.  Hence, in what follows we study only schemes with order
divisible by four.



§ OPTIMIZING THE STABILITY DOMAIN




 §.§ Introduction to ISB Optimization


Stability domain optimization has been well studied in the literature. 
The class of steppers that maximizes ISB given stability polynomial order N+1 was found
independently by Kinnmark and Gray <cit.> and  Sonneveld and van Leer <cit.>.
The methods divide a time interval into N evenly spaced steps.  A Forward Euler predictor
and Backward Euler corrector pair initiates the time step, then N-1 leap frogs
bring us to the end of the time interval.  This class of methods has order of accuracy at most
two, and achieves an ISB_n=N/(N+1).

Kinnmark and Gray demonstrate third
and fourth order accurate stability polynomials in <cit.> that converge
to the optimal ISB as number of subintervals increases.  Interestingly, the first
two methods of this class are the third order and fourth order explicit Runge-Kutta methods
with three and four stages, respectively.  Thus RK_4 is optimal in the sense that it
fully utilizes its four function evaluations to maximize time step for wave-type problems.
It is therefore an excellent candidate for comparison with the optimized schemes that follow.



 §.§ GBS Stability Domain Optimization


In Richardson extrapolation schemes one sets up a square Vandermonde
system to compute the weights guaranteeing a specified order of accuracy.
If we allow the number of components in the extrapolation scheme to
increase beyond those necessary for maintaining order of accuracy
we obtain an underdetermined system with degrees of freedom. We
utilize these degrees of freedom to optimize the stability domain
along a contour in the complex plane. 

Extrapolation allows us to eliminate successively higher order terms
in the asymptotic error expansion of our solution. To do so, for each of m
integrators we divide the time interval H into n_i steps of size h_i=H/n_i, i=1,2,,m.  We then
construct a linear system to eliminate terms through order p-1 in the error expansion,
yielding a p-order accurate solution.
In the case of GBS integrators, the odd coefficients in the asymptotic expansion
are zero.  Thus we may drop the constraint equations for odd powers
of h_i, obtaining a system of p/2 equations:

    [16pt][       1       1               1;   h_1^2   h_2^2           h_m^2;       ⋮       ⋮       ⋱       ⋮; h_1^p-2 h_2^p-2         h_m^p-2 ]_p/2× m_V[5pt][ [1]
        c_1;         c_2;            ;           ⋮;            ;         c_m ]_m
    _c
    =
    [16pt][ [1]
        1;         0;         ⋮;         0 ]_p/2.
    _b

When m=p/2 we have a square matrix which corresponds to
the usual Richardson extrapolation schemes. The matrix is invertible
when n_i n_j,i j, and so we solve for the weight vector c,
which we apply to the individual integrated solutions to form the combined
solution at the end of the time interval.

By allowing m>p/2 the system becomes underdetermined and we may
enforce order constraints while optimizing selected features of the stability
domain.  The optimization algorithm was adapted from the polynomial optimization formulation
in <cit.>; details are provided in <ref>.




 §.§ Fully-Determined Optimization Results


We first investigate optimal step count selection for fully determined
extrapolation schemes.  For schemes of order p we test each combination of p/2 step counts
up to a set maximum, here chosen to be 24.  Each combination yields a set of uniquely determined 
extrapolation weights.  We then select the combination of step counts
that maximizes ISB_n of the extrapolated stability domain.  Tab. <ref> contains the tabulated
results for orders eight, twelve and sixteen.  The schemes all have generous imaginary
axis coverage and can be implemented efficiently on three, four and five cores respectively.




  §.§.§ Eighth Order

The eighth order, three-core scheme has ISB_n = 0.5799 with the following step counts and 
uniquely determined weights:

    [1.5] Step Counts 
    Weights    [1.5] : 
     : [1.5]
        2,     16,     18,     20 
    
        -1/498960,    65536/9639,     -531441/25840,    250000/16929.





  §.§.§ Twelfth Order

The twelfth order, four-core scheme has ISB_n = 0.4515 with the following step counts and weights:

    [1.5] Step Counts 
    Weights    [1.5] : 
     : [1.5]
        2,     8,     12,     14,     16,     20 
    
        -1/157172400,     4096/155925,      -59049/15925,     282475249/15752880,      -4194304/178605,     9765625/954261.





  §.§.§ Sixteenth Order

The sixteenth order, five-core scheme has ISB_n=0.4162 and utilizes the step count
sequence {2,8,10,12,14,16,18,22}.  Extrapolation weights can be computed by solving
the corresponding Vandermonde system (<ref>).  We omit them here since the weights are ratios of
large integers in both the numerators and denominators.




 §.§ Underdetermined Optimization Results


Using the optimization methodology described in <ref> we optimize the ISB of GBS-type methods up to order
sixteen.  Increasing the number of extrapolation components leads to an
increase in ISB_n.  Since for explicit schemes
the maximum ISB_n is one we expect the relative gains of adding more components to eventually saturate.
By evenly distributing work across CPU cores we can demonstrate the relationship
between available processors and maximal time step size.
This correspondence between core count and ISB_n is shown in Fig. <ref>.
Here we observe that efficiency saturation occurs around ten cores for all orders of accuracy;
the saturation value itself is strongly dependent on the order of accuracy.


We achieve the best optimization results by utilizing all (even) step counts up to a maximum
dependent on the number of available CPU cores.  The m extrapolation
components then have subinterval counts {2, 4, , 2m}, with m set by the available processing resources.
We denote the number of subintervals at which ISB_n can no longer be increased N_opt and
aggregate the results for each order of accuracy in Tab. <ref>.
The stability domains for each core count are plotted
in Fig. <ref>.
Results demonstrate a tradeoff between optimal ISB and order of accuracy as is typical
of explicit time integrators.


The capping of N_opt is an artifact of the optimizer.  Our convex solver
fails to produce methods with larger ISB if we increase the number of free variables
beyond those presented in Tab. <ref>.  We believe that by addressing the conditioning
of the Vandermonde system as in <cit.> one can continue further along the curves presented in
Fig. <ref>.  Extrapolating these curves shows the schemes do not converge
to the optimal ISB_n=1; the exact tradeoff between order of accuracy and optimal
ISB is a topic of future research.





 §.§ Leading Order Error

Let the time integrator's stability polynomial, as defined in <ref>, be denoted R(ξ).
For a method of order p, the stability domain's boundary follows the imaginary axis
surrounding the origin linearly up to deviation on the order p+1.  To compute
the leading error coefficient we set the stability polynomial R(ξ) = e^iθ.  Taking the
complex logarithm of the polynomial and Taylor expanding yields a power series for θ(ξ).
We then compute the inverse series to find
ξ(θ) = iθ+a_p+1(iθ)^p+1+a_p+2(iθ)^p+2+𝒪(θ^p+3).
Since we consider only methods with order divisible by four we simplify as follows:

    ξ(θ) = iθ + i a_p+1θ^p+1 - a_p+2θ^p+2 + 𝒪(θ^p+3).

Departure from the imaginary axis is governed by the a_p+2 coefficient.  We then require a_p+2 to be
negative to ensure the stability domain has a positive ISB.  Accuracy is determined by
the a_p+1 coefficient; Fig. <ref> presents this coefficient
for each method as a function of number of cores.

The leading order error coefficient decays toward zero as number of cores
increases while holding order fixed.  This implies the underdetermined extrapolation
schemes do not trade away numerical precision for achieving large ISBs - they instead
gain significantly in accuracy.




 §.§ Core Partitioning

In order to achieve the theoretical efficiencies presented in Sect. <ref> we require 
a work partitioning scheme that distributes the individual time steppers amongst the cores.
For a specified maximum subinterval count N_max we achieve the largest ISB_n by utilizing all even step counts up
to N_max.  For a fixed number of cores, here denoted n_cores, we can always evenly distribute
the work when we choose N_max=4· n_cores-2.  This corresponds to folding together onto a single core
pairs of integrators with step counts adding to N_max.  For example, with N_max=10,
we evenly load three cores with step counts {10, {8, 2}, {6, 4}}.

The GBS scheme requires N+1 function evaluations for an integrator with N subintervals.
When stacking multiple integrators on a single core we share the first function evaluation
since it takes identical arguments for all time steppers.  We can in principle share this evaluation
among all cores but communication overhead may make this approach less efficient.




 §.§ Method Specification

Let the m×p/2-sized Vandermonde matrix in (<ref>) be denoted V.
When we have m>p/2 the system is underdetermined.  We
then choose p/2 dependent step counts to meet the order constraints
and collect these into the step count sequence {n_dep}.  The associated extrapolation
weight vector is denoted c_dep and the corresponding columns of V denoted V_dep.
Likewise collect the remaining step counts into the (m-p/2)-length step count sequence {n_free}
and label corresponding extrapolation weights c_free and Vandermonde columns V_free.
Then the constraint equations may be written

    V_depc_dep+V_freec_free = b.

The Vandermonde system order constraints must be met exactly.
We thus omit floating-point coefficients for c_dep in the text; they are best computed
symbolically then converted to the desired floating point format.  We may readily compute c_dep
with the relation 

    c_dep = V_dep^-1(b-V_freec_free),

so one only needs the step count sequences {n_dep} and {n_free} and  weight vector
c_free to fully specify a scheme.




 §.§ Methods of Choice

In this section we present two eighth order methods and one twelfth order method
with rational coefficients for convenient use.  In order to provide
robustness to spurious discretized eigenvalues sitting slightly in
the right-half plane we push the optimization curve into the positive
reals.  This disturbs the ISB_n very little and makes the method 
suitable for local differentiation stencils generated for example by
RBF-FD (radial basis function-generated finite difference) approximations
<cit.>.  

To generate the following schemes we first optimize the free coefficients using
the optimization methodology described in <ref>.  The optimization contour is
chosen to trade off a small amount of imaginary axis coverage for an area containing
the positive reals away from the origin.  We then perform a search
over a set of rational numbers that closely approximate the floating point free
coefficients.  We select a set with small integers in the numerator and denominator
which disturbs the scheme's stability domain very little.  These coefficients are
reported below.




  §.§.§ The Eighth Order, Six Core Method: GBS_8,6

The following eighth order method achieves a robust stability domain for wave-type PDEs
on six cores and is therefore dubbed GBS_8,6.
The scheme's ISB_n is 0.7675, a 0.26% reduction from the optimal six core value of 0.7695.
The scheme can be implemented using the following step counts and extrapolation weights:

    {n_dep}   = { 2,     4,     6,     10 }
    [1.5] {n_free}
    c_free   [1.5] = 
    = [1.5]
        {    8,     12,     14,     16,     18,     20,     22    }
    
        [    2165/767488,    13805/611712,    4553/72080,    14503/66520,    27058/7627,     -86504/5761,    40916/3367   ]^T.

The dependent weights can be computed
exactly by inverting the corresponding Vandermonde system.
The scheme's stability domain is plotted in Fig. <ref>, with a
zoom-in around the imaginary axis on the right-hand side.

This method can be run efficiently on six cores with time steps 6.25 times larger
than those of RK_4.  After normalizing for number of function evaluations, the scheme
achieves time-to-solution 0.7675/0.7071 = 8.5% faster than RK_4 but with eighth order of accuracy.  
Compared to RK_8, though, we achieve time-to-solution 269% faster.
As will be seen in Sect. <ref>, speed-up to achieve a specified accuracy
is far improved over RK_4 due to the size of the leading order error term combined
with eighth order convergence to the true solution.






  §.§.§ The Eighth Order, Eight Core Method: GBS_8,8

Likewise for the six core method, we produce a robust stability domain for an eighth order
scheme that runs efficiently on eight cores, called GBS_8,8. 
The scheme's ISB_n is 0.8176, a 0.24% reduction from the optimal eight core value, 0.8196.
The scheme utilizes the following step counts and extrapolation weights:

    {n_dep}   = { 2,     26,     28,     30 }
    {n_free}   = { 4,     6,     8,     10,     12,     14,     16,     18,     20,     22 }
    
    c_free   = 
    
    
    
    
    
    
    
    [1.5]
        [    6833/476577792,    10847/91078656,    15235/34643968,    383/321152,    543/198784,    
       [1.5]
           9947/1741056,    6243/543104,    6875/296192,    1401/28496,    17713/152688,    6375/19264   ]^T.

The scheme's stability domain is plotted in Fig. <ref>, with a
zoom-in around the imaginary axis on the right-hand side.

This method can be run efficiently on eight cores with time steps 8.96 times larger
than those of RK_4.  After normalizing for number of function evaluations the scheme
achieves time-to-solution 15.6% faster than RK_4, and 287% faster
than RK_8.  The two additional cores grant us a 6.5% increase in efficiency over GBS_8,6.






  §.§.§ The Twelfth Order, Eight Core Method: GBS_12,8

The following twelfth order method runs on eight cores and is therefore dubbed GBS_12,8.
The scheme's ISB_n is 0.7116, a 0.17% reduction from the optimal eight core value of 0.7128.
The scheme can be implemented using the following step counts and extrapolation weights:

    {n_dep}   = { 2,     8,     10,     16,     24,     26 }
    {n_free}   = { 4,     6,     12,     14,     18,     20,     22,     28,     30 }
    
    c_free   = 
    
    
    
    
    
    
    
    [1.5]
        [    235/21030240256,    4147/1612709888,    11521/39731200,    2375/3528704,    6435/708736,    
       [1.5]
           1291/15780,    11311/4672,     -180864/751,    222080/2079   ]^T.

Stable time steps with GBS_12,8 are 7.79 times larger than those of RK_4
and, after normalization, time-to-solution is improved by 0.6%.
This (very) modest efficiency improvement is drastically offset by the twelfth order of convergence
of the method - wall time to achieve a desired accuracy is far shorter than that of RK_4.






§ NUMERICAL RESULTS


The following results demonstrate the optimized time
steppers on a test problem with known analytic solutions.



 §.§ One-Way Wave Equation


To demonstrate the performance benefit over standard time steppers from the
literature we run the periodic one-way wave equation, 

    ∂∂ tu(x,t) + ∂∂ xu(x,t) = 0,        0 ≤ x ≤ 1, t>0, 
        u(x,0) = 12(1-cos2π x),        0 ≤ x ≤ 1,

utilizing the rational-coefficient GBS_8,6 and GBS_12,8 methods, with RK_4 as reference.
Spatial derivatives are spectral to ensure errors are due to the time stepping algorithm alone.  We run
all time steppers near their respective limits of stability, at λ=Δ t/Δ x=.99/π×ISB,
where the factor of π arises from the spectral spatial derivatives.
After convecting the wave once around the periodic interval we compute the absolute error with
respect to the analytic solution, then refine in both time and space.  

Convergence to the analytic solution for the various methods
is demonstrated in Fig. <ref>.  For fair comparison across methods,
the horizontal axis is time step normalized by the number of function evaluations per step.
Thus vertical slices correspond to equal time-to-solution, neglecting the overhead of sharing
data across cores.  We use a high precision floating point library <cit.> for computation since
machine precision is achieved in the high order methods before we can establish a 
trendline.  Coefficient truncation to double precision causes error to
stagnate at 10^-14 and 10^-12 for the eighth and twelfth order methods, respectively.  
To obtain full floating point precision to 10^-16 the extrapolation coefficients
must be precise to twenty significant digits.






§ CONCLUSIONS


We have presented a scheme for maximizing the time step size for extrapolation
based ODE solvers.  To do so we construct an underdetermined Vandermonde system,
then optimize the weights to maximize the stability domain along a given curve in
the complex plane.  For wave-type PDEs we utilize GBS integrators and optimize the
methods for imaginary axis coverage.  We achieve large ISB
values for methods through order sixteen which, when implemented on a computer with
several cores, yield faster time to solution than standard Runge-Kutta integrators.

The optimization method leaves both the time integrator and desired contour of stability
as user parameters.  Changing the ODE integrator in turn changes the stability polynomial
basis, which immediately affects the resulting extrapolated stability domains.  The GBS
integrator maintains large ISB through extrapolation; other integrators may be better suited
for different desired stability domains.  Future work therefore involves identifying other
suitable integrators for stability domain optimization in different contexts.  

The underdetermined extrapolation scheme saturates in parallelism
around ten cores.  We can improve scalability by incorporating the optimized schemes as local
building blocks in time-parallel solvers like Parareal <cit.>.  These solvers are known to be
less efficient with wave-type PDEs <cit.>.  Stable algorithms may be achieved by optimizing the
global time-parallel integrators rather than optimizing the coarse and fine grid propagators individually.
These degrees of freedom provide more flexibility than optimizing only the local schemes and is
a promising research direction for improving time to solution for wave-type equations.






§ STABILITY DOMAINS AND IMAGINARY AXIS COVERAGE FOR SOME STANDARD CLASSES OF ODE SOLVERS






 §.§ Stability domains and their significance for MOL time stepping


Each numerical ODE integration technique has an associated stability
domain, defined as the region in a complex ξ-plane, with ξ=h λ,
for which the ODE method does not have any growing solutions when
it is applied to the constant coefficient ODE 

    y^'=λ y .

For a one-step method the stability polynomial, here denoted
R(ξ), is the numerical solution after one step for Dahlquist's
test equation (<ref>) <cit.>.  The method's stability
domain is then 

    S={ξ∈ℂ : |R(ξ)| ≤ 1}.

When solving ODEs, the stability domain can provide a guide to the
largest time step h that is possible without a decaying solution
being misrepresented as a growing one. In the context of MOL-based
approximations to a PDE of the form ∂ u/∂ t=L(x,t,u),
the role of the stability domain becomes quite different, providing
necessary and sufficient conditions for numerical stability under
spatial and temporal refinement: all eigenvalues to the discretization
of the PDE's spatial operator L must fall within the solver's stability
domain. For wave-type PDEs, the eigenvalues of L will predominantly
fall up and down the imaginary axis. As long as the time step h is
small enough, this condition can be met for solvers that feature a
positive ISB, but never for solvers with ISB = 0. 



 §.§ Runge-Kutta methods


All p-stage RK methods of order p feature the same stability
domains when p=1,2,3,4. For higher orders of accuracy, more than p stages
(function evaluations) are required to obtain order p. The RK_4
scheme used here is the classical one, and the RK_8
scheme is the one with 13 stages, introduced by Prince and Dormand
<cit.>, also given in <cit.> Table 6.4. Their normalized stability
domains are shown in Fig. <ref>.
Their ISBs are 2.8284 and 3.7023, respectively. 







§ ISB OPTIMIZATION ALGORITHM






 §.§ Optimization Formulation

Let the extrapolated GBS stability polynomial be denoted R(ξ), and the individual
stability polynomials from each of m extrapolation components be denoted P_i(ξ).
Then we have

    R(ξ) = ∑_i=1^mc_iP_i(ξ)

for extrapolation weights c_i, 1 ≤ i ≤ m.
We collect the monomial coefficients of each P_i(ξ) into the rows of a matrix, denoted P(ξ).
Then R(ξ) can be more compactly expressed as R(ξ)=c^TP(ξ).  Now let the
left-hand-side Vandermonde matrix from (<ref>) be denoted V, and the right-hand-side
constraint vector be denoted b.  Then our order constraint equation (<ref>) can be
rewritten as Vc=b.

Denoting the time step size h, and given a curve  Λ⊂ℂ,
we specify the optimization problem as follows:

    c_1,c_2,…,c_mmaximize       h 
       subject to       |R(hλ)|-1 ≤ 0,  ∀λ∈Λ. 
              Vc = b


Following the work of Ketcheson and Ahmadia in <cit.>, we reformulate the
optimization problem in terms of an iteration over a convex subproblem.  Minimizing
the maximum value of |R(hλ)|-1 over the weights c_i is a convex problem 
(see <cit.>).  We therefore define the subproblem as follows:

    c_1,c_2,…,c_mminimize      max_λ∈Λ(|R(hλ)|-1)
       subject to        Vc = b

Calling the minimax solution to (<ref>) r(h,Λ), we can now
reformulate the optimization problem as:

    c_1,c_2,…,c_mmaximize       h 
       subject to       r(h,Λ) ≤ 0

The optimization routine was implemented with the CVX toolbox for MATLAB <cit.>
using a bisection over time step h.  Results presented in this paper use the software
OPTISB <cit.> to optimize the stability domains.



 §.§ Comparison to Optimizing Monomial Coefficients


The main theoretical difference between our current algorithm and the algorithm
presented in <cit.> is the basis over which coefficients are
optimized.  In <cit.> the authors optimize directly the coefficients
to the stability polynomial in the monomial basis.
This yields an optimal stability polynomial that must be approximated
with a Runge-Kutta integrator.  The polynomial is therefore fed into a second
optimization routine to compute the Runge-Kutta coefficients.

In the extrapolation coefficient optimization we operate directly on linear
combinations of the time stepper stability polynomials.  The true optimal stability
polynomial therefore may not be in the space of extrapolated GBS time stepper
stability polynomials.  However, the resulting stability polynomial is immediately
realizable and we require no further optimization stage to generate our time
stepping algorithm.  



 §.§ Implementing Order Constraints


To guarantee accuracy we require order constraints to be satisfied to machine
precision.  Most optimization routines accept equality constraints that will hold within a
certain tolerance.  Due to ill-conditioning of the Vandermonde systems we
prefer to explicitly enforce the order constraints in the convex optimization.
As in Sect. <ref> we split the stability polynomials into two groups which take on the “dep" and
“free" subscripts, denoting dependent and optimized quantities, respectively.  The dependent
weights guarantee the extrapolation scheme achieves the specified order of accuracy.
The remaining weights are our optimization variables.   Thus the stability
polynomial is computed as follows:

    R(ξ) = c_dep^TP_dep(ξ) + c_free^TP_free(ξ).

Order constraints take the form (<ref>) which yields the
dependent weight computation (<ref>).
Splitting the weights apart reduces the number of design variables and, in practice, leads to
better solutions than when utilizing equality constraints.





spmpsci


