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}

  • $\begingroup$ Did you implemented RK4 by yourself, or did you use some built-in function ? (odeint,ode45 ,...)? You have different scales in your solution, so classical explicit RK methods may not perform so good $\endgroup$
    – VoB
    Commented Dec 30, 2020 at 14:05
  • $\begingroup$ Yes, I implemented it myself. I'll add this with the Butcher tableau in the question. I need a quite simple implementation, since it will be done in parallel on a GPU in the end. For that reason, I'd like to avoid any fancy stuff like variable timesteps etc. $\endgroup$ Commented Dec 30, 2020 at 14:54
  • $\begingroup$ @VoB, also as far as I understand, the problem is not stiff and I don't have any stability problems. So I'd rather not go for implicit methods unless there is some other advantage besides stability. $\endgroup$ Commented Dec 30, 2020 at 15:06
  • $\begingroup$ My bad, reading your post I assumed you had different scales. You can certainly save same evaluations of $A(t_n + \frac{h}{2})$ in your code. Btw, in which sense are you doing this in parallel? The stages are sequential @johannestoger $\endgroup$
    – VoB
    Commented Dec 30, 2020 at 15:26
  • $\begingroup$ @VoB I'm solving for many different $\omega(t)$ in parallel. $\endgroup$ Commented Dec 30, 2020 at 15:28

1 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:

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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.