Hello. New to Julia, so my question is tentative (and might be being posted in the wrong place). I am looking for singularity detection features in the solvers in DifferentialEquations.jl, and I am not seeing any. The usual heuristics in non-Julia solvers (eg ode45 in Matlab will report one when t+dt==t) are based on the step-size going to zero, and they’re pretty reliable, and remarkably accurate. But the messages I get from the Julia solvers are like the first one in this post: “There is either an error in your model specification or the true solution is unstable”. While that might be technically true (a solution that is singular won’t be stable in the usual sense) that word isn’t used for singular solutions, and indeed what I want to know is the location of the singularity.
Simple examples: u’ = t^2 + u^2 has a singularity between Pi/4 and 1. Matlab’s ode45 locates it nicely.
I’m very impressed with the Julia codes overall (and by now I have cited your work, Chris, several times because of the breadth of coverage you have) and I am now trying to learn how to use them effectively. So, if there is a known method for locating singularities, I would like to learn the syntax. While I am new to Julia, one of the problems is that I am an old-fashioned programmer and quite a lot of the words being used are unfamiliar to me (e.g. “callback” for event location, but I have started to get that to work); I suspect that the unfamiliarity of the words is because they are modern usage and I’m just old. So even just learning the words you use to describe the features you have for singularity detection and handling would be helpful for me in reading the documentation and locating what I want.
julia> using OrdinaryDiffEq
julia> prob = ODEProblem((u, p, t) -> t^2 + u^2, 1.0, (0.0, 1.0))
ODEProblem with uType Float64 and tType Float64. In-place: false
timespan: (0.0, 1.0)
u0: 1.0
julia> sol = solve(prob) ;
┌ Warning: dt(1.1102230246251565e-16) <= eps(t)(0.9697906948802055) , and step error estimate = 2.5757996612379013. Aborting. There is either an error in your model specification or the true solution is unstable (or the true solution can not be represented in the precision of Float64).
└ @ SciMLBase C:\Users\____\.julia\packages\SciMLBase\roUH9\src\integrator_interface.jl:633
julia> sol.t[end]
0.9697906948802055
Expanding a bit on @rveltz’s answer, one way of finding the singularity might be to just solve, ignore the warning and take the last time of the returned solution (the solution is still saved, even though it reports an instability).
With an even simpler example (where the analytic function is known) it could look like this:
using OrdinaryDiffEq
using Plots
function f!(du, u, p, t)
du[1] = u[1]^2
end
ode = ODEProblem(f!, [1.0], (0.0,10.0))
sol = solve(ode, Tsit5())
singularity = last(sol.t)
scatter(sol.t, only.(sol.u); label="solve", ylim=(1.0, 1e3), yscale=:log10)
plot!(sol.t, (t -> 1 / (1 - t)).(sol.t), label="analytic")
vline!([singularity], linestyle=:dash, label="singularity")
In my understanding it’s just not possible to know a priori if the user just made a typo or if the ODE is correct and singular. I would guess that the former case is a lot more common, hence the way the output is presented feels more like a failure. But if you already expect the solution to have a singularity, all the output is there.
Since I have no idea about how Matlab does it, could you share an example of how the output looks like there? I am also curious to learn if there is an option in OrdinaryDiffEq.jl that allows one to return the singularities (but I’m not aware of any).
(Larry Shampine and I wrote a paper about the Maple codes in 2000: Redirecting . Larry wrote the solvers (I think Chris had some questions about the choices, in his “Comparison” paper) and I wrote the graphical interface—which later got re-done by others. Larry’s code has been extended and maintained in-house by Allan Wittkopf.)
I originally typed that as one “reply” but the system would not let me (a new user) put two screenshots and a DOI link in the same box, or reply too quickly after I had split them (within 13 seconds, which I find amusing). I expect that over time, as I get used to the community here, my posts will become more (shall we say) culturally attuned!
Hey sorry was out for a few days. It seems like the issue here is just the wording of the message? One thing I can tell you is that users will never believe they have done anything incorrectly unless you’re explicit about it. If the warning message said “your model is likely singular”, I can guarantee you that the standard response would then be to write in the paper “and Julia’s differential equation solver reported that the equation is singular”. I’m almost at 10 years of daily debugging people’s code online and that’s pretty much a given .
So the “There is either an error in your model specification” is very purposefully put in there, and it was added later. Probably 2018 or so? And if you search the old vs the new error message, you see a really dramatic drop in the number of people who report the error. And if you look through the trends you will see why. Before that was added, there were many people who said things like “my model went to dtmin, how do I force it to stay above dtmin?”, and the thread would eventually find they used x instead of -x and then it magically didn’t need to change dtmin anymore because now the protein level does not go to infinity. After it was added, well people now treat dtmin issues as an error and presumably check their model before positing.
So I understand that dt->dtmin can be used on valid models to do things like detect a singularity, but I think we have to be really cautious about what we tell most users. Most users will take for granted whatever the warning says, and if in most cases if you’re integrating over a singularity, it probably wasn’t there on purpose . But yes for people who really know dynamics and are trying to estimate the location of a singularity, your model specification is fine and the true solution is unstable. Since that or is in there, “or the true solution is unstable (or the true solution can not be represented in the precision of Float64)”, I think the error message is okay, but it very purposefully suggests you to check your model first.
