Playing with Runge-Kutta methods and SageMath.

To use the code, you must install SageMath.

The goal of this code is to verify different properties of a given Runge--Kutta method, defined by it's Butcher array (all the notions used here are defined in the Bible, see below).

• With SageMath, it is possible to make exact computations in some useful sets of numbers, and thus we make the hypothesis that the coefficients of the methods live in the set of (real) algebraic numbers (called AA in SageMath); this is not a really restrictive hypothesis. But then, the result obtained are proofs (if both my codes and SageMath are correct!).

• To compute the order of a Runge-Kutta method, one use the so called rooted trees for which we have a SageMath implementation, coded by Florent Hivert.

• We also provide a function to compute the Butcher array of a method defined by collocation. A classical application is the set of Gaussian Runge-Kutta methods.

Some sage/jupyter notebooks are provided: they probably provide the best way to learn and to test this code (even, you can run the notebooks on binder, see below).

Note that the cost of the computation grows very fast with tne number of steps of the method; also, the rooted_tree machinery in sage is quite slow.

References:

• To learn about Ordinary Differential Equations solvers, you should read the Bible:

1- Solving Ordinary Differential Equations I, by Hairer, Nørsett,, Wanner (HNW),

2- Solving Ordinary Differential Equations II Stiff and Differential-Algebraic Problems by Hairer and Wanner (HW),

3- Geometric Numerical Integration by Hairer, Lubich and Wanner (HLW).

If you want to learn and to understand rooted trees and B-series, I recommend to start first by reading reference 3. Being the latest book, things have become much easier to understand than in reference 1.

• If you want to learn SageMath, you can read the book Mathematical Computation with Sage (which now is available in French, English and German), and for which freely available pdf files can be downloaded here and there.

Using the code

• First, you must define a Runge--Kutta method. To do this you must write a (simple) python class: have a look at formulas.py in methods/ to undestand what to do.

  The class must derive from the "RungeKutta" class and include a constructor.

The constructor must:

1- Define the arrays A and B of the Butcher array (the C part is generally not necessary).

2- Give a title.

3- Call the base RungeKutta class constructor.

Remember that the "formula.py" file in methods/ give examples of such classes.

Be careful: The coefficients of A and B must be algebraic numbers (AA or QQbar) or rational numbers (QQ). This is not absolutely trivial since, in sage:

2/3, for example, is directly evaluated as a float (0.66666...). This is a Python intrinsic, and cannot be avoided. So as you need to enter an algebraic (real) number or a rational, you must write QQ(2/3) or AA(2/3) or 2/QQ(3) or 2/AA(3) instead of 2/3. Again, have a look at methods/formula.py.

Have a look at "GoodAndBad.ipynb": this notebook shows this more in details.

• Then, the best is to look at the notebook Example1.ipynb and adapt it to your needs. For this, launch sage like this:

>sage -n jupyter

or

>sage -n jupyterlab

and then, launch the notebook Exemple1.ipynb.

Gaussian formulas (and other methods obtained by collocation)

Gausian formulae with n steps are obtained by collocation at the roots of the Legendre P polynomials of degree n, shifted from [-1,1] to [0,1].

RKcolloc.colloc computes the Butcher arrays, given collocation points in [0,1] and returns a Runge-Kutta method class (note that collocation points are not necessary Gaussian points).

See the notebook Gaussian.ipynb.

For the fun, formulas obtained by collocation at roots of Tchebychev polynomials are computed in "Tchebychev.ipynb" (they certainly have no interesting application).

Notebooks

All are Sage/Jupyter notebooks. Launch Sage by typing:

sage -n jupyter

or

sage -n jupyterlab

1. Example.ipynb : a tour of the system.

2. Gaussian.ipynb: construct and test some of the Gauss methods.

3. AllProperties.ipynb: show how to compute all possible properties of a method.

4. Test.ipynb : just testing that everything works fine on a set of formulas.

5. RadauByCollocation.ipynb : Radau methods are collocation methods. Play with this.

6. GoodAndBad.ipynb : is supposed to show what to do and what not to do when coding a Runge-Kutta formula.

Running the notebooks on Binder

Just Click me.

otherwise, if it does not work, have a look at "Dockerfile".

Implementation:

• The code uses a decorator @_persistance to avoid recomputing known properties (which is often expensive).

A precedent version was using Sage's @lazy_attribute decorator, which could be disturbing.