Strange result of multiplication

[MartinOtter]

When evaluating Measurements.jl with DifferentialEquations.jl on a simple bouncing ball example, a very strange effect occurs (performing a multiplication outside of this example gives a different result as in this example):

Here is the bouncing ball example with Measurements.jl

module Test_bouncingBall_with_measurement

using DifferentialEquations, Measurements

function f(du,u,p,t)
    du[1] = u[2]
    du[2] = -9.81
end

function condition(u,t,integrator)
    z = u[1] + 1e-12
    return z
end

const e = 0.7 ± 0.1   # coefficient of restitution

function affect_neg(integrator)
    println("Event at time = ", integrator.t, ", h = ", integrator.u[1])

    v_before = integrator.u[2]
    println("   v_before = ", v_before, ", e = ", e)
    v_after = -e*v_before
    println("   v_after  = ", v_after)

    integrator.u[2] = v_after
    auto_dt_reset!(integrator)
    set_proposed_dt!(integrator, integrator.dt)
end

cb = ContinuousCallback(condition,nothing,affect_neg! = affect_neg)

u0 = [1.0 ± 0.2, 0.0]
tspan = (0.0,2.0)
prob = ODEProblem(f,u0,tspan)
sol = solve(prob,Tsit5(),saveat=0.01,callback=cb)

end

Results in the following output:

Event at time = 0.4515236409859562, h = -0.0 ± 0.2
   v_before = -4.42944691807223 ± 0.0, e = 0.7 ± 0.1
   v_after  = 3.1 ± 0.44
Event at time = 1.0836567383662918, h = -0.0 ± 0.34
   v_before = -3.1 ± 0.44, e = 0.7 ± 0.1
   v_after  = 2.170428989855373 ± 2.8e-15

Note that v_after = -e*v_before and it seems surprising that the uncertainty after this operation is nearly zero (2.8e-15). However, when perform this operation separately (outside of this example), the result is the correct one:

using Measurements

v_before = -3.1 ± 0.44
e        =  0.7 ± 0.1
v_after  = -e*v_before
@show v_after

Resulting output:

v_after = 2.17 ± 0.44

It is strange that when applying the same operation in the affect_neg callback, a wrong result is calculated. Is this a bug in Measurements.jl or did I made a mistake?

[giordano]

Your example doesn't run for me, I get the following error:

ERROR: LoadError: MethodError: no method matching Float64(::Measurements.Measurement{Float64})
Closest candidates are:
  Float64(::Real, ::RoundingMode) where T<:AbstractFloat at rounding.jl:200
  Float64(::T) where T<:Number at boot.jl:718
  Float64(::Int8) at float.jl:60
  ...

with a very long stacktrace pointing to line

sol = solve(prob,Tsit5(),saveat=0.01,callback=cb)

Edit:

Ok, I can run the code if I change

tspan = (0.0,2.0)

to

tspan = (0.0 ± 0.0, 2.0 ± 0.0)

[giordano]

I think that the differential equation example is working correctly. What you're seeing is elimination of uncertainty because Measurements.jl takes care of functional correlation, so for example:

julia> x = 0.7 ± 0.1
0.7 ± 0.1

julia> x - x
0.0 ± 0.0

julia> (0.7 ± 0.1) - (0.7 ± 0.1)
0.0 ± 0.14

julia> x / x
1.0 ± 0.0

julia> (0.7 ± 0.1) / (0.7 ± 0.1)
1.0 ± 0.2

For example, you can use Measurements.derivative to see how the derivative of v_before with respect to e changes over time:

module Test_bouncingBall_with_measurement

using DifferentialEquations, Measurements

function f(du,u,p,t)
    du[1] = u[2]
    du[2] = -9.81
end

function condition(u,t,integrator)
    z = u[1] + 1e-12
    return z
end

const restitution = 0.7 ± 0.1   # coefficient of restitution

function affect_neg(integrator)
    println("Event at time = ", integrator.t, ", h = ", integrator.u[1])

    v_before = integrator.u[2]
    @show v_before
    @show Measurements.derivative(v_before, restitution)
    v_after = -restitution*v_before
    @show v_after
    @show Measurements.derivative(v_after, restitution)

    integrator.u[2] = v_after
    auto_dt_reset!(integrator)
    set_proposed_dt!(integrator, integrator.dt)
    println()
end

cb = ContinuousCallback(condition,nothing,affect_neg! = affect_neg)

u0 = [1.0 ± 0.2, 0.0]
tspan = (0.0 ± 0.0, 2.0 ± 0.0)
prob = ODEProblem(f,u0,tspan)
sol = solve(prob,Tsit5(),saveat=0.01,callback=cb)

end

Output:

Event at time = 0.45152364098595665 ± 0.0, h = -0.0 ± 0.2
v_before = -4.42944691807223 ± 0.0
Measurements.derivative(v_before, restitution) = 0.0
v_after = 3.1 ± 0.44
Measurements.derivative(v_after, restitution) = 4.42944691807223

Event at time = 1.0836567383662938 ± 0.0, h = -0.0 ± 0.34
v_before = -3.1 ± 0.44
Measurements.derivative(v_before, restitution) = 4.42944691807223
v_after = 2.17042898985538 ± 1.8e-15
Measurements.derivative(v_after, restitution) = -1.7763568394002505e-14

Event at time = 1.5261499065325295 ± 0.0, h = -0.0 ± 0.34
v_before = -2.170428989855389 ± 1.8e-15
Measurements.derivative(v_before, restitution) = -1.7763568394002505e-14
v_after = 1.52 ± 0.22
Measurements.derivative(v_after, restitution) = 2.1704289898554014

Event at time = 1.8358951242488952 ± 0.0, h = -0.0 ± 0.4
v_before = -1.52 ± 0.22
Measurements.derivative(v_before, restitution) = 2.1704289898554014
v_after = 1.0635102050291443 ± 3.3e-16
Measurements.derivative(v_after, restitution) = -3.3306690738754696e-15

Your second example doesn't give you the same result because this doesn't take into account that v_before actually has a functional relation with e, i.e. a non-null derivative.

Does this make any sense?

[MartinOtter]

O.k. Thanks for your explanation. Yes this makes sense (so everything is corret, just my intuitive understanding was wrong).

[giordano]

Anyway, do you have an analytic solution for this problem to compare with?

[MartinOtter]

Anyway, do you have an analytic solution for this problem to compare with?

You are right, this should be checked. Since the equations are so simple, I derived the analytic solution manually now:

module Test_bouncingBall_with_measurement_analytic

using Measurements

g = 9.81
e = 0.7 ± 0.1
h0 = 1.0 ± 0.2

#=
der(h) = v  ; h(0) = h0
der(v) = -g ; v(0) = 0

when h <= 0 then
	v = -e*v
end when

Analytic solution of differential equation with
h(t_ev) = h_ev, v(t_ev) = v_ev:

	v = -g*(t-t_ev) + v_ev
	h = -g/2*(t-t_ev)^2 + v_ev*(t-t_ev) + h_ev

First phase (t_ev = 0, v_ev=0, h_ev=h0)
	v = -g*t
	h = -g/2*t^2 + h0

First event:
	hev1 = 0 -> t_ev1 = sqrt(2*h0/g)
	            v_ev1_before = -g*t_ev1
				v_ev1_after  = -e*v_ev1_before

Second phase (t_ev = t_ev1, v_ev = v_ev1_after, h_ev = 0)
    v = -g*(t-t_ev1) + v_ev1_after
    h = -g/2*(t-t_ev1)^2 + v_ev1_after*(t-t_ev1)

Second event:
	h_ev2 = 0 = -g/2*(t-t_ev1)^2 + vev1_after*(t-tev_1)
	          = (t-t_ev1)*(-g/2*(t-t_ev1) + v_ev1_after)

	-> t_ev2-t_ev1 = 2*v_ev1_after/g

	v_ev2_before = -g*(t_ev2-t_ev1) + v_ev1_after
	             = -g*2*v_ev1_after/g + v_ev1_after
	   			 = -v_ev1_after
	v_ev2_after  = -e*v_ev2_before
=#

t_ev1        = sqrt(2*h0/g)
v_ev1_before = -g*t_ev1
v_ev1_after  = -e*v_ev1_before
t_ev2        = t_ev1 +2*v_ev1_after/g
v_ev2_before = -v_ev1_after
v_ev2_after  = -e*v_ev2_before

@show t_ev1
@show v_ev1_before
@show v_ev1_after
@show t_ev2
@show v_ev2_before
@show v_ev2_after

end

Resulting in

                               # Numeric solution:
t_ev1        = 0.452 ± 0.045   # 0.45152364098595665 ± 0.0
v_ev1_before = -4.43 ± 0.44    # -4.429446918072235 ± 0.0
v_ev1_after  = 3.1 ± 0.54      # 3.1 ± 0.44
t_ev2        = 1.08 ± 0.14     # 1.083656738366296 ± 0.0
v_ev2_before = -3.1 ± 0.54     # -3.1 ± 0.44
v_ev2_after  = 2.17 ± 0.66     # 2.1704289898553957 ± 1.3e-16

Hm. So, this is not the result computed by DifferentialEquations.jl and Measurements.jl and even the uncertainty at the first event is 20 % different to the analytic solution. I simulated with reltol=1e-8 and saveat=0.0001, and also tried another solver (BS5) but the numeric solution did not change much.

So, there is an issue.

The problem seems to be that the analytic solution computes an uncertainty in the event times (t_ev1, t_ev2), whereas the solution with DifferentialEquations.jl/Measurements.jl has no uncertainty in the event times.

[giordano]

A couple of comments:

• my understanding is that the time doesn't propagate the uncertainty (maybe it's just sampled at regular times?), so your times should be

  t_ev1        = Measurements.value(sqrt(2*h0/g))
  t_ev2        = Measurements.value(t_ev1 +2*v_ev1_after/g)

  With this change I get as output (extending to 4 bounces)

  t_ev1 = 0.4515236409857309
  v_ev1_before = -4.4294469180700204
  v_ev1_after = 3.1 ± 0.44
  t_ev2 = 1.083656738365754
  v_ev2_before = -3.1 ± 0.44
  v_ev2_after = 2.17 ± 0.62
  t_ev3 = 1.5261499065317703
  v_ev3_before = -2.17 ± 0.62
  v_ev3_after = 1.52 ± 0.65
  t_ev4 = 1.8358951242479817
  v_ev4_before = -1.52 ± 0.65
  v_ev4_after = 1.06 ± 0.61
  

  which agrees with

  julia> e =  0.7 ± 0.1
  0.7 ± 0.1
  
  julia> v_ev1_before = -4.42944691807223 ± 0.0
  -4.42944691807223 ± 0.0
  
  julia> v_ev1_after = -e * v_ev1_before
  3.1 ± 0.44
  
  julia> v_ev2_before = -v_ev1_after
  -3.1 ± 0.44
  
  julia> v_ev2_after = -e * v_ev2_before
  2.17 ± 0.62
  
  julia> v_ev3_before = -v_ev2_after
  -2.17 ± 0.62
  
  julia> v_ev3_after = -e * v_ev3_before
  1.52 ± 0.65
  
  julia> v_ev4_before = -v_ev3_after
  -1.52 ± 0.65
  
  julia> v_ev4_after = -e * v_ev4_before
  1.06 ± 0.61

  Does this make sense?

• if you look at the output of my example above you can note that when the ball touches the ground, v_before has opposite sign compared to the v_after at the previous iteration, but the derivative has the same sign, instead of having the sign flipped as well. I have no idea what's going on inside the solver. As a workaround I can suggest you the following hack: set v_before as

  v_before = Measurements.result(integrator.u[2].val, -1, integrator.u[2])

  in order to flip the sign of the derivative. In this way I get the following result:

  Event at time = 0.45152364098595665 ± 0.0, h = -0.0 ± 0.2
  v_before = -4.42944691807223 ± 0.0
  v_after = 3.1 ± 0.44
  
  Event at time = 1.0836567383662938 ± 0.0, h = -0.0 ± 0.34
  v_before = -3.1 ± 0.44
  v_after = 2.17 ± 0.62
  
  Event at time = 1.5261499065325295 ± 0.0, h = -0.0 ± 0.59
  v_before = -2.17 ± 0.62
  v_after = 1.52 ± 0.65
  
  Event at time = 1.8358951242488952 ± 0.0, h = -0.0 ± 0.78
  v_before = -1.52 ± 0.65
  v_after = 1.06 ± 0.61
  

I'd tentatively say that there is nothing wrong in Measurements.jl but the solver is missing to correctly propagate the derivate for some weird reasons, which is not unreasonable since that package doesn't know anything about Measurements.jl and it's already crazy that they can mostly work together.

Maybe we can convince @ChrisRackauckas to have a look at this? 🙂

[giordano]

The more I look at this, the more I'm convinced that the solver is working fine, but maybe something weird is happening when the callback is invoked?

[ChrisRackauckas]

What does Measurements.jl do during the rootfind? My guess is that rootfinding is not something that plays well with Measurements? We're calling Roots.jl here.

[MartinOtter]

Here is the code to compute t_ev1 with Roots.jl:

using Roots
f(t)  = -g/2*t^2 + h0
t_ev1 = find_zero(f,(0.1±0.0, 0.5±0.0))
@show t_ev1
end

The result is

t_ev1 = 0.452 ± 0.045

so exactly the result of the analytic solution. This indicates that Measurements.jl is doing it right.

A potential reason of the issue could be that the bounds t1,t2 given to find_zero(.., (t1, t2)) in DifferentialEquations are Float64 instead of Measurement{Float64}. But maybe, it is also something else.

[giordano]

What does Measurements.jl do during the rootfind? My guess is that rootfinding is not something that plays well with Measurements? We're calling Roots.jl here.

It looks like we have the same guess 🙂

If in the definition of zero_func at https://github.com/JuliaDiffEq/DiffEqBase.jl/blob/2386f211a2228bb302a0145c039ffa9ca2d2b450/src/callbacks.jl#L539-L542 I add a @show Θ I see the following output:

Θ = 0.0 ± 0.0
Θ = 1.0 ± 0.0
Θ = 0.5 ± 0.0
Θ = 0.513 ± 0.085
Θ = 0.509 ± 0.088
Θ = 0.509 ± 0.088
Θ = 0.509 ± 0.088
Θ = 0.509 ± 0.088
Θ = 0.509 ± 0.088
Θ = 0.509 ± 0.088
Θ = 0.5089510673864017
Θ = 0.508951067386402

I verified these lines come from the find_root at https://github.com/JuliaDiffEq/DiffEqBase.jl/blob/2386f211a2228bb302a0145c039ffa9ca2d2b450/src/callbacks.jl#L561. The fact that Θ is changing type is suspicious 🤔

I need to investigate more, but it looks like find_root is the culprit here.

[MartinOtter]

When changing the line in https://github.com/JuliaDiffEq/DiffEqBase.jl/blob/2386f211a2228bb302a0145c039ffa9ca2d2b450/src/callbacks.jl#L561 in the following way:

# Θ = prevfloat(find_zero(zero_func, (bottom_θ,top_Θ), Roots.AlefeldPotraShi(), atol = callback.abstol/100))
 Θ = find_zero(zero_func, (bottom_θ,top_Θ), Roots.AlefeldPotraShi(), atol = callback.abstol/100)

so removing the prevfloat call, then the numeric solution with DifferentialEquations.jl gives the result of the analytic solution:

Event at time = 0.452 ± 0.045, h = -9.99e-13 ± 3.5e-15
   v_before = -4.43 ± 0.44, e = 0.7 ± 0.1
   v_after  = 3.1 ± 0.54
Event at time = 1.08 ± 0.14, h = -9.982e-13 ± 2.1e-15
   v_before = -3.1 ± 0.54, e = 0.7 ± 0.1
   v_after  = 2.17 ± 0.66
Event at time = 1.53 ± 0.27, h = -1.0002e-12 ± 3.7e-15
   v_before = -2.17 ± 0.66, e = 0.7 ± 0.1
   v_after  = 1.52 ± 0.67
Event at time = 1.84 ± 0.39, h = -1.0e-12 ± 1.1e-14
   v_before = -1.52 ± 0.67, e = 0.7 ± 0.1
   v_after  = 1.06 ± 0.62

So, this is clearly an issue of DifferentialEquations.jl and I will close this issue here.