Trying to use Measurement.jl with DifferentialEquations.jl

Specific Domains Modelling & Simulations
, ,
1

In the example below I am trying to use Measurement.jl with DiffernetialEquations.jl to solve an ODE with error propagation.

clearconsole()
@time using DifferentialEquations, Measurements, MATLAB


function Buckle(du,u,p,t)
R0,  f,  ps, kai, ks, Rbreak, sr0,alphapa=p
rhoL=998;
vL=1e-3;
c=1481;
pinf0=1.01e5;
k=7/5;

SR=0
if u[1]>=(R0/sqrt(sr0/kai +1))
    SR=@. kai*(((u[1]/(R0/sqrt(sr0/kai +1)))^2)-1)
elseif u[1]>= Rbreak*R0
    SR=S
end


du[1]=u[2];
du[2]=(1/(rhoL*u[1]))*(-1.5*rhoL*u[2]^2+(pinf0+2*sr0/R0)*((R0/u[1])^(3*k))*(1-3*k*u[2]/c)-4*vL*u[2]/u[1]-4*(ks/(1+alphapa*abs(u[2])/u[1]))*u[2]/(u[1]^2)-2*SR*exp(-kais*(u[1]/R0-1))/u[1]-pinf0-ps*sin(2*pi*f*t));

end
#################
### Parameters ##
#################
R0=1e-6
ps=10e3
f=1e6±0.001e-5

sr0=0.05
kai=1.8
kais=0
ks=1.8e-8
Rbreak=1.01
alphapa=0
###############
###############
###############
p=[R0,  f,  ps, kai, ks, Rbreak, sr0,alphapa,kais]
u0=[R0,0]
tspan=(0,1/f)
prob = ODEProblem(Buckle,u0,tspan,p)
@time sol = solve(prob,SSPRK104(),dt=1e-10,saveat=0.01/f,reltol=1e-10,abstol=1e-10)

I can only use ± on R0 in my parameter list. I encounter an error stating

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

when I try to include ± for ps, f, sr0, kai, kais, ks, Rbreak and alphapa. Am I doing something wrong or is my model incompatible with this option?

2

Make your initial condition be a Measurement so that the uncertainties can propagate along the solution.

3

Like u0=[R0±0, 0±0] ?

4

Yes, and probably also the time span:

tspan = (0 ± 0, 1 / f)
5

You shouldn’t need to do that unless you’re using continuous event handling.

6

f is already a Measurement though

7

Oh, then it should just promote correctly anyways.

1 Like
8

Thank you. This solved the issue.

9

this works well without a callback now. what should I change to make it compatible with DiffEqCallbacks?

using DiffEqCallbacks
u0=[R0±0,0±0]
tspan=(0,1/f)
saved_values = SavedValues(Float64, Array{Float64})
cb = SavingCallback((u,t,integrator)->integrator(t,Val{1}), saved_values, saveat = 0:0.01/f:1/f)
prob = ODEProblem(Buckle,u0,tspan,p)
@time sol = solve(prob,SSPRK104(),dt=1e-11,saveat=0.01/f,reltol=1e-10,abstol=1e-10,callback=cb)
10

You just said right here that the solution should be a Float64, but you want it to be a Measurement.

1 Like
11

So how do I change it to Measurement? It might be obvious but I’m a rookie at Julia now.

12

I think I found it

saved_values = SavedValues(Float64, Array{Measurement{Float64}})

fixes it