Problem understanding results of DiffEqUncertainty

[MrBellamy]

I have issues using the librairie DiffEqUncertainty.jl or at least understanding the results. I have a DAE made of 4 differential equations and 1 algebric equations. I defined it using a mass matrix and I am solving it with an ODE solver like this

solve(prob_MM, Rodas5(), reltol = 1e-8, abstol=1e-12, callback=cbs, maxiters = 1e7)

The solver manage to solve it without issue :
9.159245 seconds (21.20 M allocations: 1.964 GiB, 9.29% gc time, 0.10% compilation time)

and I can plot the solution. It looks like that

which is what I was expecting. So far so good. Now I tried to see how sensitive my solutions was, so I used the librairie DiffEqUncertainty.jl with the lines

cb_uncert = AdaptiveProbIntsUncertainty(5)
ensemble_prob = EnsembleProblem(prob_MM)
sim = solve(ensemble_prob, Rodas5(), reltol = 1e-10, abstol=1e-12, trajectories=5, callback=cbs)
(cb_uncert is added to cbs using cbs = CallbackSet(others_cb, cb_uncert ))

Then the simulation takes a really long time with A LOT of allocation

6588.658787 seconds (21.72 G allocations: 3.086 TiB, 52.83% gc time, 0.00% compilation time)

and the solution is really not smooth, like this :

I was expecting results similar in smoothness as those ones https://diffeq.sciml.ai/stable/analysis/uncertainty_quantification/
What could be the cause of this behaviour ? Is my system extremely unstable or is it something else ? Why does it take so many allocations ? I would say the specificity of my system is that I interpolate many thermodynamical properties, so the code is quite heavy and not very practical to post here.

[ChrisRackauckas]

@MrBellamy

I defined it using a mass matrix and I am solving it with an ODE solver like this

solve(prob_MM, Rodas5(), reltol = 1e-8, abstol=1e-12, callback=cbs, maxiters = 1e7)

Don’t use Rodas5 for this. The documentation says why. Try Rodas4 here.

[MrBellamy]

What I found in the doc about Rodas5 is that : “Currently has a Hermite interpolant because its stiff-aware 3rd order interpolant is not yet implemented”. Is this the reason why I should not use it ? I don’t understand it. Anywway, I tried Rodas4 as suggested, it took much longer :

16113.129405 seconds (23.94 G allocations: 3.447 TiB, 74.66% gc time, 0.00% compilation time)

but the results seem a bit better though really not smooth

How come it takes so long for the uncertainty calculations when individually the simulation is relatively fast (9.1 s) and why does it take so much allocations ? I guess it is linked with the non-smooth nature of the solution but I cannot understand why ?

[ChrisRackauckas]

@MrBellamy

What I found in the doc about Rodas5 is that : “Currently has a Hermite interpolant because its stiff-aware 3rd order interpolant is not yet implemented”. Is this the reason why I should not use it ?

It’s much more explicit than that in the docs.

https://diffeq.sciml.ai/stable/solvers/dae_solve/#Rosenbrock-Methods

• Rodas5 - A 5th order A-stable stiffly stable Rosenbrock method. Currently has a Hermite interpolant because its stiff-aware 3rd order interpolant is not yet implemented. This means the interpolation is inaccurate on algebraic variables, meaning this algorithm should not be used with saveat or post-solution interpolation on DAEs.

So if the non-smooth behavior was in the algebraic variables, then that would be the explanation. If it’s in the differential variables, and that explosion… I think you might be doing something wrong in your f function. See:

PSA: How to help yourself debug differential equation solving issues

Debugging differential equation solver issues almost always boils down to doing the same thing, so this is a summary to help you out. For more information on specific issues, check out the FAQ section of the DifferentialEquations.jl documentation which highlights common issues and questions: https://diffeq.sciml.ai/dev/basics/faq/ How do I debug why the differential equation solver is diverging? dt <= dtmin. Aborting. There is either an error in your model specification or the true solution i…

One guess, do you happen to modify u?

[MrBellamy]

Indeed it is more explicit than that, I only looked at this page :
https://diffeq.sciml.ai/stable/solvers/ode_solve/
where the explanation was more succinct.

I had looked at your PSA post which is very useful, thanks for that, unfortunately I checked already those points. The problem is that my code works perfectly when I solve it “alone” and most of the tips on this page are for ODE system that cannot get solved.

I am not modifying u or f either in my code. I had a “max” function in my f function that I transformed into a smooth max using 0.5* x * (1 + erf(x/sqrt(2.0))) - Wolfram|Alpha
and I relaunched the uncertainty code but it took even longer and displayed similar behaviour :

Could it be a normal behaviour linked with the stiffness of my model ?

[ChrisRackauckas]

Nah, that’s not a stiffness thing. That looks like a caching problem. If you can open an issue with your code it could be looked into more.

[MrBellamy]

Ok, I cleaned the code a bit and posted an issue on the github :
https://github.com/SciML/DiffEqUncertainty.jl/issues/57

I hope it is the proper place, I am quite new to this.

[ChrisRackauckas]

That’s the right place, thanks.