, 3 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 \in\mathbb{C}, \quad t_0, t\in\mathbb{R},$$

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}, \quad h\in\mathbb{R}, \quad k,m \in\mathbb{N}, $$

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 the 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$, where $\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.

Added 09-May-2025: To get the ball rolling, a first simple approach is to implement a command line program, which takes a list of formulas and stepsize changes and produces JavaScript output for Apache ECharts. That is then displayed in the browser.

Added 08-Jul-2025: For the program stabregion.c we use $m=1$. Therefore

$$ \det Q(\mu,H) = \det\left[ (A_0-HB_0) + (A_1-HB_1)\mu \right] = 0, \quad H=\lambda h\in\mathbb{C} $$

When $\mu=e^{i\varphi}$ traverses the unit circle, then framing above equation as generalized eigenvalue problem

$$ Ax=\lambda Bx, $$

we get

$$ A = A_0 + \mu A_1, \quad B = B_0 + \mu B1, $$

and solve for $H=H(\mu)$. This gives the stability regions.

If, however, we want to depict the stability mountains, we solve for $\mu=\mu(H)$, and have

$$ A = A_0 - H B_0, \quad B = -(A_1 - H B_1). $$