Magnus solver slow compared with Tsitouras for linear state independent ODE

Specific Domains Numerics
,
1

Hi! I was playing a bit with the Magnus solvers included in the examples for linear ODEs in the docs: Non-autonomous Linear ODE / Lie Group ODE Solvers.

I was particularly interesting in benchmarking the solver I am currently using for a state independent solver and I noticed that the gold standard Tsitouras solver still performs better in these cases, even when Magnus solvers should be in principle better for ODEs of the form u' = A(t) u.

Here is a simple MWE with the example included in the documentation:

using OrdinaryDiffEqCore, OrdinaryDiffEqTsit5
using OrdinaryDiffEqLinear
using BenchmarkTools

function update_func(A, u, p, t)
    A[1, 1] = cos(t)
    A[2, 1] = sin(t)
    A[1, 2] = -sin(t)
    A[2, 2] = cos(t)
end

A = DiffEqArrayOperator(ones(2, 2), update_func = update_func)
prob = ODEProblem(A, ones(2), (1.0, 6.0))

@btime solve(prob, Tsit5(), dt=0.1, abstol=1e-6, reltol=1e-6);
# 26.245 μs (280 allocations: 24.20 KiB)

@btime solve(prob, MagnusGL6(), dt=0.1, abstol=1e-6, reltol=1e-6);
# 208.365 μs (3481 allocations: 335.88 KiB)
@btime solve(prob, MagnusGL4(), dt=0.1, abstol=1e-6, reltol=1e-6);
# 93.735 μs (981 allocations: 94.47 KiB)
@btime solve(prob, MagnusLeapfrog(), dt=0.1, abstol=1e-6, reltol=1e-6);
# 67.176 μs (532 allocations: 45.33 KiB)

Also, the Magnus solvers don’t allow don’t passing dt, reason why the comparison is with fixed timestep.

Does this result makes sense? I was expecting the Magnus solvers to be the best solution for my problem u' = A(t)u, but it looks like Tsitouras also is better here!

Any feedback is more than welcome, I am trying to understand if I am missing something here.

1 Like
2

for small problems like this, a good out of place solver is really hard to beat.

3

First of all, you need adaptive=false here to force fixed dt on Tsit5. You only set the initial dt.

The Magnus methods are non-adaptive so the tolerance does not effect them. Their error is determined by the dt.

The Magnus methods will could use some more optimization. But I think for very small cases it will be hard for them to outperform the highly optimized RK methods. Where they shine is not in terms of performance but instead enforcing that certain constraints are satisfied, like positivity.

1 Like
4

They should be norm preserving, for example, if the Lie group is SO(N)

5

Thank you all for the responses.

I made those changes and, as expected after this discussion, the results had not really change:

@btime solve(prob, Tsit5(), dt=0.1, adaptive=false);
# 38.446 μs (487 allocations: 41.09 KiB)

@btime solve(prob, MagnusGL6(), dt=0.1);
# 208.365 μs (3481 allocations: 335.88 KiB)
@btime solve(prob, MagnusGL4(), dt=0.1);
# 93.735 μs (981 allocations: 94.47 KiB)
@btime solve(prob, MagnusLeapfrog(), dt=0.1);
# 67.176 μs (532 allocations: 45.33 KiB)

I guess I was expecting to gain some performance with the Magnus methods but that was not the case :sob:

Thank you all!

6

We can probably optimize them some more, but I wouldn’t expect them to generally be faster.

1 Like