Special-case Runge-Kutta methods to exploit structure in linear ODE?

Question

I am interested in numerical solutions of a linear, time-dependent ODE of the form

$$ \dot y = A(t)y - Ry, $$

A good model is the following problem in $\mathbb R^2$:

$$ A(t) = \begin{bmatrix}0 & -\omega(t)\\ \omega(t) & 0\end{bmatrix} \\ R = \begin{bmatrix}R_1 & 0\\ 0 & R_2\end{bmatrix}, \\ $$

The full problem is in $\mathbb R^3$ and a bit more complicated, but I think the above example captures most of it.

The solution has two main features: an "oscillating" part and an exponential decay. The function $\omega(t)$ can take on values as high as $10^6$, while the decay factor $R$ is on the order of $10^0$ to $10^2$. Therefore, the oscillation is very rapid in comparison to the exponential decay. However, $\omega(t)$ has a quite slow rate of change, on the order of $10^1$ to $10^2$ per unit time.

I tried the classical 4th-order Runge-Kutta method, and it works quite well. One advantage is that it only requires samples of $A(t)$ at the beginning, middle and end of each time step, which simplifies my implementation.

Now to my question(s): Is it possible to exploit the structure of this system in some way to obtain a higher-performance method? Specifically, can a "better" Runge-Kutta method (higher order, fewer stages, lower principal error term constant) be derived using the additional structure? Or would it be better to use some other kind of solution strategy?

Edit: Here's the Butcher tableau for the method I'm using. Constant timesteps for simplicity since it will run in parallel for many different instances of $\omega(t)$ on a GPU in the end. I wrote my own implementation and I have verified that a log-log plot of error vs. timestep has slope 4.

\begin{array}{c|cccc} 0 & & & & \\ 1/2 & 1/2 & & \\ 1/2& & 1/2 & \\ 1 & & & 1 \\ \hline & 1/6 & 1/3 & 1/3 & 1/6 \end{array}

Answer 1

There are many kinds of RK methods that have extensions to exploit linearity. They all use some form of exponential or Lie Group idea (again exponential) to do so. Thus they generally do some form of integrating factor method and then apply the Runge-Kutta method to the IF transformed equation. If you want to see a listing of some of these, you can check out the methods found in DifferentialEquations.jl.

There are standard RK methods that achieve superconvergence on linear problems, i.e. hit an order higher. DiffEqDevTools.jl implements almost every tableau I could find in the literature, and from that, we found the following superconvergence cases:

• Dormand and Prince's 6th order method
• Tsitouras and Papakostas's 6th order method
• Verner's old 6th order method
• Enright and Verner's 7th order method

You can dig up the papers for these methods by finding them in the 8,000 lines of coefficients tableau page. Looking at the papers, this superconvergence behavior seems to be by accident. But it does show it's possible to achieve an order higher than expected on linear equations.

Now onto the more sophisticated methods. In the Split ODE Solvers section you will find the methods compatible with the form:

$$u' = Au + f(u,p,t)$$

also known as a semilinear ODE. These are known as the exponential Runge-Kutta methods (exponential RK). This is probably the most common use, but there are others.

If you just have $u' = Au$ then, of course, $u(t) = \exp(TA)u(0)$ but you wouldn't want to just use the scaling and squaring methods built into most languages so you'd want to use something that exploits the Krylov subspace for large exponential times vector calculations, hence the LinearExponential method.

Another common case are methods that specialize in $u' = A(t)u$, which is what most Magnus methods specialize in. You can find the list of methods for this category as the state-independent solvers. The problem you write downfalls into the Magnus methods since you can make $\hat{A}(t) = A(t) - R$.

Then there are the Lie Group methods, like the Runge–Kutta–Munthe-Kaas method, which specializes in the form $u' = A(u)u$. Then of course there are a few methods that allow for time and state dependence, i.e. $u' = A(u,t)u$, specifically some Crouch-Grossman Runge Kutta methods, but as you specialize less it tends to get less and less benefit.

As with any of the DifferentialEquations.jl methods, you can from the Julia REPL do ?alg on the algorithm to dig up the original source if you want more information.

I want to leave with one last detail though that people tend to ask about this kind of thing. Let's say you have an affine equation:

$$u' = A(t)u + b(t)$$

could you still exploit the linear structure with a Magnus method? The answer is yes. What you need to do is exploit the fact that an affine transformation is just a linear transformation in a larger space. I.e. you define $\hat{u}$ such that $\hat{u} = [u;1]$, so you append a 1 to your original space, and then $\hat{A}(t)$ has a block that is $A$, with a row of 0's appended on the bottom and the column $b$ appended to the right. If you work it out you'll see that the extended portion of the equation has an initial condition 1 and 0 derivative, and so $\hat{A}\hat{u} = A(t)u + b(t)$ on the non-constant portion. This means you can also encode affine transformations into a form to use these specialized RK methods.

So Magnus, RK of Crouch–Grossman type, Runge–Kutta–Munthe-Kaas, exponential Runge-Kutta, Cayley methods, and more. There are tons of methods that exploit this property, and there are libraries that implement them. In particular, DifferentialEquations.jl uses specialized matrix exponential methods, is GPU-compatible, automatically constructs adjoints, etc. That said, this is all far from complete and there are many more things that can be done. But that's the most complete listing that I know of. For more on the backgrounds of these methods, you might want to check out Hairer's Geometric Numerical Integration book which has a very good discussion on a subset of the methods described here.