[New algorithm] Runge-Kutta in the interaction picture (RKIP) #2736

Azercoco · 2025-06-10T15:33:28Z

The RKIP is a method to solve SplitODEProblem where f1 is a AbstractSciMLOperator which has an efficient way of computing/storing exp(A*t). It performs a change of variable to turn the (stiff) problem into a non-stiff one which can be solved with any Runge-Kutta methods.

The RKIP method is well suited for pseudospectral method for solving elliptic PDE (such as the ones from FourierFlows.jl).

This method is similar to the exponential Runge-Kutta (like ETDRK4) but allows for adaptative time stepping and arbitrary order. In counterpart, it requires a constant linear operator (so no Jacobian can be used).

oscardssmith · 2025-06-10T16:21:09Z

lib/OrdinaryDiffEqRKIP/src/algorithms.jl

+
+KEYWORD_DESCRIPTION = """
+- `nb_of_cache_step::Integer`: the number of steps. Default is `100`.
+- `tableau::ExplicitRKTableau`: the Runge-Kutta Tableau to use. Default is `constructDormandPrince6()`.


why is this the default? The Dormand Prince integrators tend to not be very good.

I need to recheck but from my benchmark on NLSE-like equations, it seems to perform the best. What others methods would you suggest ?

Tsit5 or Vern6/7 tend to be good see https://docs.sciml.ai/SciMLBenchmarksOutput/stable/NonStiffODE/FitzhughNagumo_wpd/ for details.

the tableau methods do not have the same performance and that isn't showing accuracy, just tolerance.

Comparison with lower order method. The slopes seem all correct. The last plot has possibly a too long integration time compared to Lyapunov time. The only thing strange is that the order 4 and 5 method seems to be have the same slope. But from this, I would suggest having either CashKarp (for medium accuracy) or Verner7 (for high accuracy) as the default tableau.

It's very odd for Cash Karp to be good, I wonder if it has some odd super convergence on this problem 😅

That's entirely possible that it is a specific of this problem. Also, the error control is made on the transformed equation so the change of coordinate can increase (or decrease) the error. So you would suggest to set Verner6 as default ?

I think that would be a generally good choice.

Verner6 is now the default tableau.

src/OrdinaryDiffEq.jl

ChrisRackauckas · 2025-06-17T19:37:51Z

lib/OrdinaryDiffEqRKIP/test/type_test.jl

+# fails as calling a MatrixOperator on  a StaticVector change its type to Vector, causing a type instability
+
+# @testset "Out-of-place Diagonal 1x1 Operator SVector Adaptative" begin
+#     test((μ) -> DiagonalOperator([μ]), u0 -> SVector(u0), false)
+# end
+
+# @testset "Out-of-place Matrix 1x1 Operator SVector Adaptative" begin
+#     test((μ) -> MatrixOperator([μ;;]), u0 -> SVector(u0), false)
+# end
+
+# # Fails for a strange reason as calling ldiv!(a, B, c) on a MatrixOperator dispatch toward an inexisting method
+
+# @testset "In-Place MatrixOperator Vector Float Adaptative Ldiv" begin
+#     test(μ -> MatrixOperator([μ;;]), u0 -> [u0], true; use_ldiv=true)
+# end


Set as @test_broken

ChrisRackauckas · 2025-06-17T19:38:15Z

lib/OrdinaryDiffEqRKIP/test/runtests.jl

+@safetestset "Type Safety Tests" include("type_test.jl")
+@safetestset "Cache Test" include("cache_recycling_test.jl")
+@safetestset "Fourier Semilinear PDE Tests" include("semilinear_pde_test_cpu.jl")
+@safetestset "CUDA Test" include("semilinear_pde_test_cuda.jl")


This needs to move to the CUDA queue

This mean moving the test to OrdinaryDiffEq/test/gpu ?

lib/OrdinaryDiffEqRKIP/src/algorithms.jl

ChrisRackauckas · 2025-06-17T19:43:26Z

lib/OrdinaryDiffEqRKIP/src/algorithms.jl

+solve(ode_prob_1, rkip, t1)
+solve(ode_prob_2, rkip, t2)
+````
+
+will reuse the precomputed exponential cached during the first `solve` call.
+This can be useful for solving several times successively problems with a common `A`.


t1 and t2 aren't defined here. And this is interface breaking: it should just use init and reinit! / solve! and it'll do this without the extra work anyways.

ChrisRackauckas · 2025-06-17T19:44:06Z

lib/OrdinaryDiffEqRKIP/src/algorithms.jl

+This solver is only implemented for semilinear problem: `SplitODEProblem` when the first function `f1` is a `AbstractSciMLOperator` A implementing:
+
+```julia
+LinearAlgebra.exp(A, t) # = exp(A*t)


This should instead be using ExponentialUtilities.jl algorithm choice for exponential!

As the implementation of the exponential is very operator dependent, are you sure that is should not use the AbstractSciMLOperator interface instead ?

cf : https://github.com/SciML/SciMLOperators.jl/blob/7ba386430a229776b41f481ff352eacd7c9f09d4/src/interface.jl#L382C1-L382C60

I checked the code of ExponetialUtilities and it seems to be useful when the exponential matrix output is dense. But in that case, caching the matrix as made here does not make a lot of sense and one would be better to used an ETD method with Krylov based expmv.

As the implementation of the exponential is very operator dependent, are you sure that is should not use the AbstractSciMLOperator interface instead ?

ExponentialUtilities.jl will specialize on claimed properties like symmetric.

I checked the code of ExponetialUtilities and it seems to be useful when the exponential matrix output is dense. But in that case, caching the matrix as made here does not make a lot of sense and one would be better to used an ETD method with Krylov based expmv.

You mean sparse? expv is for the sparse case. But exponential! is for the dense case, and it will be faster for standard dense matrices.

Co-authored-by: Christopher Rackauckas <[email protected]>

Azercoco added 5 commits June 10, 2025 15:41

First commit, draft implementation of RKIP

cf38dcd

Reformat Style

60016de

Use generic_solver_docstring for the documentation

0cb3e46

fix typos and improve docs

a42e3b4

Add aligned in math mode

b2902e9

oscardssmith reviewed Jun 10, 2025

View reviewed changes

Azercoco added 3 commits June 12, 2025 11:28

Fix cache recycling with test and change default tableau to Verner6

d22b707

Dispatch p and t to the operator

5e40016

Explicitely use Verner6 as the defaut tableau

6c7924a

ChrisRackauckas reviewed Jun 18, 2025

View reviewed changes

Azercoco and others added 2 commits June 18, 2025 11:14

Fix typo

d9a8e3f

Co-authored-by: Christopher Rackauckas <[email protected]>

Remove direct export

8343936

Co-authored-by: Christopher Rackauckas <[email protected]>

Uh oh!

[New algorithm] Runge-Kutta in the interaction picture (RKIP) #2736

Are you sure you want to change the base?

[New algorithm] Runge-Kutta in the interaction picture (RKIP) #2736

Uh oh!

Conversation

Azercoco commented Jun 10, 2025

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Azercoco Jun 18, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Azercoco Jun 18, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Azercoco Jun 18, 2025 •

edited

Loading

Azercoco Jun 18, 2025 •

edited

Loading