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},


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}=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

Leave a Reply

Please log in using one of these methods to post your comment: Logo

You are commenting using your account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s