, 2 min read

Design Notes on System for the Analysis of Order- and Stepsize Changes for Cyclic Composite Multistep Methods

For the numerical solution of the ordinary differential equation, initial value problem

$$ y'=f(t,y), \quad y(t_0)=y_0,$$

one can use a cyclic composite multistep method of the form

$$ \sum\limits_{i=0}^{m} A_i Y_{k+i} = h \sum\limits_{i=0}^{m} B_i Y'_{k+i}, \qquad A_i,B_i\in\mathbb{R}^{\ell\times\ell},$$

and

$$ Y_k = \left( \begin{array}{c} y_{k\ell+1} \\ \vdots \\ y_{k\ell+\ell} \end{array} \right), \quad Y'_k = \left( \begin{array}{c} y'_{k\ell+1} \\ \vdots \\ y'_{k\ell+\ell} \end{array} \right), \quad y_i=y(t_i), \quad y'_i=f(t_i,y_i).$$

$y'$ denotes the derivative of $y$, $h$ is the stepsize, $\ell$ is the cycle length. The stepsize does not need to be kept constant. If the stepsize changes one can still use above formula but adds a sub-cycle of the form

$$ \tilde{y_i}=\hbox{linear combination / interpolation polynomial with previous grid elements}$$

Some notes:

  1. Individual stages of the cyclic composite multistep method do not have to be of same order.
  2. Some stages can be explicit, some implicit.

The method can still be stiffly stable even if some stages are explicit, assuming that the explicit stages are not used for advancing the solution.

For the test equation

$$ y'=\lambda y, \quad \lambda\in\mathbb{C},$$

one has to analyze the roots of the determinant of the matrix polynomial

$$ \det Q(\mu,H)=\det \left[ \sum\limits_{i=0}^{m}\left(A_i-H B_i\right)\mu^i \right], \quad H=\lambda h.$$

The analysis system or workbench we are after will now

  1. use the companion matrix to solve $\det Q(\mu,H)=0$

  2. use $\tilde{y_i}$ for the interpolation step, i.e., stepsize change

  3. order change is trivial

  4. calculate Widlund wedge

  5. picture of 3D images of $\det Q(\mu,H)$ and conic section of $\det Q(e^{i\varphi},H)=0, >\varphi\in\left[0,2\pi\right]$

  6. propose a list of stepsize and order changes, which can be edited easily, possibly with mouse pointer or direct edit

    1. linear increase / descrease
    2. sawtooth / zig-zag
    3. Sierpinski shape

  7. ideally be able to change some of the $A_i, B_i$ or automatically tune them with a minimization procedure according some criteria, like

    1. Widlund wedge
    2. maximal stepsize possible not exceeding minimum Widlund wedge
    3. sum of absolute values of stage coefficients

There are still a number of open issues and several decisions have to be taken:

  1. use Octave, Maxima, or C
  2. use OpenGL, or Gnuplot
  3. incorporate CUDA in C or Octave (see: Drop-in Acceleration of GNU Octave), if dimension of companion matrix becomes large and calculations therefore slow

Added 02-Jun-2020: At the time of posting mid 2013 the Julia programming language was just almost freshly born. Now it fits in nicely here: It offers speed, graphics, and support for CUDA.