In DifferentialEquations.jl, how can I find out the default options for a specific solver? In particular, I’d like to know the default options for the ODE solver DP5. I see a description of options in the docs, but I don’t see what are the defaults.
I’m asking because I’m trying to recapitulate results that I’ve previously calculated in Matlab with the ode45 function, which I believe is the same solver as DP5. The default options in Matlab are documented here.
I find the results in DifferentialEquations.jl do not exactly match Matlab’s ode45 because different time steps are chosen. I’ve attempted to make the time steps the same by setting the same abstol and reltol parameters as in Matlab, but this does not seem sufficient as DifferentialEquation.jl still tends to use less time steps. In some problem instances, I have to set abstol and reltol to be much more sensitive at 1e-8 to more closely match the final solutions generated using only 1e-4 in Matlab, but then the running time becomes slower. I suspect that changing other options like the maximum time step (qmax) in DifferentialEquations.jl may help in equalizing the results.
There are two things that are going on. For one, what MATLAB returns to you are not the actual steps. It lies. Well, it doesn’t lie, it is explained in the MATLAB ODE Suite paper https://www.mathworks.com/help/pdf_doc/otherdocs/ode_suite.pdf that the steps of a 5th order method are generally too large, so it would make bad plots. Therefore, a 4th order interpolation is used to fill in 2 intermediate points at 1/3 and 2/3 of the step. If you’re trying to match the stepping behavior of DiffEq to the output behavior of MATLAB, you can see why this is doomed to fail. Nonetheless, DiffEq just returns a continuous solution sol(t) instead, so it’s just an API difference here.
The second thing is that MATLAB sets dtmax to be the total time divided by 100, so you’re guaranteed 100 steps. We don’t care if you solve it in 5 steps as long as its to tolerance. So you’d need to change that to match MATLAB more closely.
As the other commenter mentioned, you could just do adaptive=false,tstops=my_times in order to exactly specify all of the dts for the integration. That said, you’ll want to only step to every 3rd timestep of MATLAB’s output of course .
Moral of the story, trying to make the time stepping behavior exactly the same between a closed source code and DiffEq is going to be hard because there are a lot of variables going on. We haven’t even mentioned the PI-adaptivity going on yet which could be different, but we cannot know without breaking the license… Instead of benchmarking like this, do it fairly with a work-precision plot:
Generate a reference solution with abstol=reltol=1e-14
Solve both with a range of tolerances, and calculate the error according to some norm against the reference.