Finding Roots with Roots.jl on ODE solutions

Specific Domains Numerics
,
1

I can find roots of an ODE solution nicely in this working, simplified example

using DifferentialEquations
using Plots
using Roots
prob = ODEProblem(
    function(u::Float64,p,t)
	return -u
    end,
    1.0,
    (0.0,1.0))
sol = solve(prob)
root find_zero(t -> sol(t) - 0.5, (0.0, 1.0))

But when I try this more complex Julia 0.6 example (they seemingly changed function names for 1.0 so it does not run there) I get an error:

using OrdinaryDiffEq
using Roots

function dgl(dr,r,p,t)
    a=1
    dr[1] = r[4]
    dr[2] = r[5]
    dr[3] = r[6]
    dr[4] = -a*sqrt(r[4]^2+r[5]^2+r[6]^2)*r[4]
    dr[5] = -a*sqrt(r[4]^2+r[5]^2+r[6]^2)*r[5]
    dr[6] = -a*sqrt(r[4]^2+r[5]^2+r[6]^2)*r[6] - 1
end

r0 = [0.0; 0.0; 0.0; 100; 0.0; 100]
tspan = (0.0,10.0)
prob = ODEProblem(dgl, r0, tspan)
sol = solve(prob, Tsit5())

t = sol.t
z = [el[3] for el in sol.u]
i = find(z .> 0)[end]
println(t[i])   # ≈ 5.0
println(t[i+1]) # ≈ 5.9
find_zero(t -> sol(t), (5.0,5.9), Bisection())

Error Message:

WARNING: sign(A::AbstractArray) is deprecated, use sign.(A) instead.
Stacktrace:
 [1] depwarn(::String, ::Symbol) at ./deprecated.jl:70
 [2] sign(::Array{Float64,1}) at ./deprecated.jl:57
 [3] init_state(::Roots.Bisection, ::Roots.DerivativeFree, ::Tuple{Float64,Float64}) at /home/christian/.julia/v0.6/Roots/src/bracketing.jl:181
 [4] #find_zero#12(::Roots.NullTracks, ::Bool, ::Array{Any,1}, ::Function, ::Function, ::Tuple{Float64,Float64}, ::Roots.Bisection) at /home/christian/.julia/v0.6/Roots/src/bracketing.jl:367
 [5] find_zero(::Function, ::Tuple{Float64,Float64}, ::Roots.Bisection) at /home/christian/.julia/v0.6/Roots/src/bracketing.jl:364
 [6] execute_request(::ZMQ.Socket, ::IJulia.Msg) at /home/christian/.julia/v0.6/IJulia/src/execute_request.jl:180
 [7] (::Compat.#inner#14{Array{Any,1},IJulia.#execute_request,Tuple{ZMQ.Socket,IJulia.Msg}})() at /home/christian/.julia/v0.6/Compat/src/Compat.jl:332
 [8] eventloop(::ZMQ.Socket) at /home/christian/.julia/v0.6/IJulia/src/eventloop.jl:8
 [9] (::IJulia.##15#18)() at ./task.jl:335
while loading In[3], in expression starting on line 25

DimensionMismatch("Cannot multiply two vectors")

Stacktrace:
 [1] *(::Array{Float64,1}, ::Array{Float64,1}) at ./linalg/rowvector.jl:184
 [2] init_state(::Roots.Bisection, ::Roots.DerivativeFree, ::Tuple{Float64,Float64}) at /home/christian/.julia/v0.6/Roots/src/bracketing.jl:182
 [3] #find_zero#12(::Roots.NullTracks, ::Bool, ::Array{Any,1}, ::Function, ::Function, ::Tuple{Float64,Float64}, ::Roots.Bisection) at /home/christian/.julia/v0.6/Roots/src/bracketing.jl:367
 [4] find_zero(::Function, ::Tuple{Float64,Float64}, ::Roots.Bisection) at /home/christian/.julia/v0.6/Roots/src/bracketing.jl:364
 [5] execute_request(::ZMQ.Socket, ::IJulia.Msg) at /home/christian/.julia/v0.6/IJulia/src/execute_request.jl:180
 [6] (::Compat.#inner#14{Array{Any,1},IJulia.#execute_request,Tuple{ZMQ.Socket,IJulia.Msg}})() at /home/christian/.julia/v0.6/Compat/src/Compat.jl:332
 [7] eventloop(::ZMQ.Socket) at /home/christian/.julia/v0.6/IJulia/src/eventloop.jl:8
 [8] (::IJulia.##15#18)() at ./task.jl:335
2

Looking at your code you actually want

find_zero(t -> sol(t)[3], (5.0,5.9), Bisection())

otherwise Roots doesn’t know which element of the vector you want to be zero.

1 Like
3

BTW: to make this work in Julia 1.0 the only thing that you need to do is replace the line

i = find(z .> 0)[end]

with

i = findlast(x -> x > 0, z)
4

Thank you, you’re right
find was renamed to findall actually

5

If you need the roots of the solution I think that using callbacks may be more accurate. (With a continuous callback you can force the solution to pass through the root)
See http://docs.juliadiffeq.org/latest/features/callback_functions.html

6

Depends on what you’re doing. Discontinuous changes to the ODE at specific event points? Yes, use a callback. Checking the points where the function crosses zero afterwards to generate a Poincure diagram? Interpolate afterwards.

1 Like
7

@SebastianM-C thanks for the tip – I’ll read
@ChrisRackauckas also thanks. I’m definitely only generating diagrams from the roots of the solution:)