Blog
A 3-stage 3rd order L-stable DIRK scheme with "nice" coefficients
I often solve stiff PDEs. Stiff PDEs have large stable exponents and require implicit methods to take reasonable timesteps. There has been extensive research into operator splitting, in which one handles the stiff part of dynamics implicitly and the nonlinear parts explicitly. However, I am a simple man. I like to understand my numerical integrators exactly and have them be as general as possible. I would rather treat everything implicitly, even the nonlinear bits.
Assuming the reader is familiar with Runge-Kutta methods, here is a Butcher table for a 3-stage 3rd order method with nice stability properties. It is a Diagonally Implicit Runge Kutta (DIRK) scheme.
This scheme has several features.
It is 3rd order; it solves all four order conditions.
It is L-stable.
It is a Fehlberg-inspired method. A general 3-stage DIRK has 9 degrees of freedom but the order conditions and L-stability are only 5 constraints. Extra Fehlberg conditions b1 = 0 and sum_j aij cj = 1/2 ci^2 for i=2,3 make one of the order conditions redundant.
The coefficients are rational and simple.
Here is a proof of 3rd order accuracy applied to the Kuramoto-Sivashinsky equation and a visualization of the stability function.
Periodic orbits of the 3-body problem
Last Christmas break, I tried to finalize some code for converging periodic orbits of the 3-body problem. I investigated the equal mass case. Candidates were identified from symplectic forward shooting and converged with Newton-lsqr. The whole code for converging 3D orbits is in C and it is blazing fast.
The C code for 3D solutions can be found here. I have spectral MATLAB code that computes these things globally and to machine precision. The MATLAB code can solve for the monodromy matrix and the Lyapunov exponents to machine precision with Newton-GMRES as well.
The end goal of this is to study Covariant Lyapunov Vectors (CLV) in this setting, since CLVs are exactly computable on periodic and relative periodic orbits.