Algorithm comparisons: Brusselator model

[mforets]

This might be an interesting problem to study the trade-off between algorithm choices. Consider the Brusselator model. Below I compared validated_integ and validated_integ2.

const A = 1.0
const B = 1.5
const B1 = B + 1

@taylorize function brusselator!(du, u, p, t)
    x, y = u
    x² = x * x
    aux = x² * y
    du[1] = A + aux - B1*x
    du[2] = B*x - aux
    return du
end

# set of initial states
U0(r) = Singleton([1.0, 1.0]) ⊕ BallInf(zeros(2), r)

# initial-value problem
bruss(r) = @ivp(u' = brusselator!(u), u(0) ∈ U0(r), dim: 2)

prob = bruss(0.1)

# validated_integ2
@time sol_3 = solve(prob, T=25.0, alg=TMJets(orderT=6, orderQ=3, adaptive=false, absorb=false));
  2.283377 seconds (16.94 M allocations: 1.497 GiB, 44.55% gc time)

# validated_integ
@time sol_2 = solve(prob, T=30.0, alg=TMJets21a(orderT=6, orderQ=2));
  1.345497 seconds (4.99 M allocations: 418.477 MiB, 68.60% gc time)

However, the first algorithm fails if the time window is slightly larger:

@time sol_3 = solve(prob, T=28.0, alg=TMJets(orderT=6, orderQ=3, adaptive=false, absorb=false));

I'm confused as to why this happens, since you can check by plotting with plot(sol_3, vars=(1, 3)) that the flowpipe is effectively shrinking, so my gut feeling is that it would be easier to validate those final steps.

[lbenet]

Thanks for reporting. I can reproduce that validated_integ finishes the integration (up to T=30), but validated_integ2 stops just short before T=25; for simplicity I am using adaptive=true. As far as I have checked, this occurs for orderQ=3 as well as for orderQ=2. I'm not yet sure what's happening, but suddenly the remainder produced by validated_integ2 becomes very large. Give me some time to dig into this.