SplitODEProblem efficiency against classical ODEProblem

General Usage Performance
,
1

Hello,

I read the documentation of the DifferentialEquations.jl package.

I had an ODE problem where the RHS is f(u, p, t) = ( L + L(t) ) * u, where L and L(t) are two sparse matrices (usually 600x600 matrix with a lot of zeros) and of course L(t) depends on time.
I’m yet able to solve this problem as a classical one, using the Tsit5() method it tooks about 2 seconds with 300x300 sparse matrices.

Reading the doc I noticed that this problem can be solved by the SplitODEProblem with several methods (for example ETD2()) which require a fixed time step.
This second method tooks about 7 seconds to solve the same problem with the same resolution. My question is why we need to use SplitODEProblems, when the classical one is much faster?

Here a minimal example of both methods:

function drho_dt(rho, p, t)
    A_l, w_l = p
    return (L + A_l * sin(w_l * t) * L_t) * rho
end

prob = ODEProblem(drho_dt, rho0_vec, tspan, p)
@time sol0 = solve(prob, Tsit5(), rtol = 1e-8, atol = 1e-8)

and

function f2(drho, rho, p, t)
    A_l, w_l = p
    mul!(drho, A_l * sin(w_l * t) * L_t, rho)
    #return (A_l * sin(w_l * t) * L_t) * rho
end

function update_func(A,u,p,t)
    A .= A_l * sin(w_l * t) * L_t
end

A1 = DiffEqArrayOperator(L)
A2 = DiffEqArrayOperator(0 * L_t, update_func = update_func)
prob = SplitODEProblem(A1, A2, rho0_vec, tspan, p)
@time sol = solve(prob, ETDRK2(), rtol = 1e-8, atol = 1e-8, dt = (tspan[2] - tspan[1]) / 4000)
2

One is adaptive and one is not.

1 Like
3

ETDRK2 is a low order and fairly inefficient method. It’s probably not recommended here. Did you try something like KenCarp4, specializing the derivative on the IMEX sparsity?

1 Like
4

How can i specialize the derivative on the IMEX sparsity?

Is this case a semilinear problem, expressed here?

I tried to simply replace ETDRK2() with KenCarp4(). It tooks the same time with totally wrong results.

5

Yeah sorry, just caught back up from JuliaCon. Your functions allocate a lot BTW so that would be the first thing to improve. Also using Octavian.jl, etc.

But yes, what you’d do is

prob = SplitODEProblem(f1, f2, rho0_vec, tspan, p)
@time sol = solve(prob, KenCarp47())

and it would be semi-implicit. Then just like how you’d specify sparsity in the ODEProblem:

ODEFunction(f,jac_prototype=sparsity_pattern)

you would do

SplitODEFunction(f1,f2,jac_prototype=sparsity_pattern)

and it will specialize on the sparsity, where here it’s just the sparsity of the f1 part.

I need to work on a doc tutorial around the IMEX methods showcasing this, but let me know if that’s not enough to go off of.

BTW, rtol,atol → reltol,abstol.

2 Likes