Hi all. Is there a way to deal with this problem in Julia?
function explicit(y,p,t)
sqrt(1-y[1]^2)
end
tspan = (0.0,π)
y0 = 0.0
problem = ODEProblem(explicit,y0,tspan)
sol = solve(problem, Rosenbrock23())
In MATLAB, Cleve Moler solve this problem by add a bound like this:
f = @(t,y) sqrt(1-min(y,1)^2).
This could also be used in Python with Scipy.odeint.
def f(t,y): return np.sqrt(1-np.min(y,1)**2)
But Julia don't allow this.
Hi everyone, first off, thanks for developing an awesome DE library/framework. It's really impressive just how much functionality you all have crammed into this library. I hope I'm in the right place to ask this.
I'm attempting to install the diffeqpy package to utilize the DifferentialEquations.jl library as a sort of "drop-in" solver to our already existing systems biology models written in Python since we need a speedup and easier parallelization options.
I'm on an Ubuntu 18.04 system and I've
- Installed Julia through the official binaries
- Installed the
DifferentialEquations.jlpackage throughPkg - Installed
PyCall.jland built it with the Python binary I'm using with mypipenvvrirtual environment. - Installed the
diffeqpypackage to the virtual environment viapipenv install(backended with pip, just in the virtual environment) - Written benchmarking scripts for a system of Lorenz Attractors as described in https://github.com/SciML/diffeqpy written using:
a. SciPy via Python 3.6
b.DifferentialEquations.jlvia Julia
c.diffeqpyviapython-jl
d.diffeqpyviapython-jlwith only benchmarking the solving of the problem, not the specification viade.ODEProblem()
e.diffeqpyviapython-jlwith only benchmarking the solving of the problem as in (d.) and also pre-compiling the problem vianumba - Executed the scripts and plotted the results.
These are the results I've obtained:
- Diffeqpy took 13.09746799999997 seconds to complete
- Diffeqpy Precompiled took 12.395007999999985 seconds to complete
- Diffeqpy Numba Compiled took 12.467692999999949 seconds to complete
- SciPy ODEINT took 0.004549000000000001 seconds to complete
- DifferentialEquations.jl took 0.050286285 seconds to complete
(tried to share graphically but couldn't figure out how to do that)
You can see the full code at https://github.com/dcolli23/benchmarking_diffeqpy
I'm almost certain I've done something wrong here since this performance is vastly different and I know DifferentialEquations.jl to be some of the fastest implementations of ODE solvers out there. Do you all have any idea what I could have done wrong in my setup of the tech stack or if I'm just interfacing with the solver incorrectly?
I seriously appreciate your all's help in all of this! Thank you so much! And I'm happy to provide any extra information as needed!
Hello, I am trying to fit an ODE with 8 external inputs from measurement data which are interpolated and I find BFGS very slow compared to BlackBoxOptim. I also tried ADAM from DiffEqFlux and it is also very slow. Does anyone here have any thoughts on why ?
I am using LinearInterpolation() objects to get the inputs at each time t
I have more details here
https://discourse.julialang.org/t/bfgs-very-slow-compared-to-blackboxoptim-how-to-improve-performance/49253/14
_itp are LinearInterpolation() objects from Interpolations.jl, the rest are constants enclosed in a wrapper function.function thermal_model!(du, u, p, t)
# scaling parameters
p_friction, p_forcedconv, p_tread2road, p_deflection,
p_natconv, p_carcass2air, p_tread2carcass, p_air2ambient = p
fxtire = fxtire_itp(t)
fytire = fytire_itp(t)
fztire = fztire_itp(t)
vx = vx_itp(t)
alpha = alpha_itp(t)
kappa = kappa_itp(t)
r_loaded = r_loaded_itp(t)
h_splitter = h_splitter_itp(t)
# arc length of tread area
theta_1 = acos(min(r_loaded - h_splitter, r_unloaded) / r_unloaded)
theta_2 = acos(min(r_loaded, r_unloaded) / r_unloaded)
area_tread_forced_air = r_unloaded * (theta_1 - theta_2) * tire_width
area_tread_contact = tire_width * 2 * sqrt(max(r_unloaded^2 - r_loaded^2, 0))
q_friction = p_friction * vx * (abs(fytire * tan(alpha)) + abs(fxtire * kappa))
q_tread2ambient_forcedconv = p_forcedconv * h_forcedconv * area_tread_forced_air * (t_tread - t_ambient) * vx^0.805
q_tread2ambient_natconv = p_natconv * h_natconv * (area_tread - area_tread_contact) * (t_tread - t_ambient)
q_tread2carcass = p_tread2carcass * h_tread2carcass * area_tread * (t_tread - t_carcass)
q_carcass2air = p_carcass2air * h_carcass2air * area_tread * (t_carcass - t_air)
q_carcass2ambient_natconv = p_natconv * h_natconv * area_sidewall * (t_carcass - t_ambient)
q_tread2road = p_tread2road * h_tread2road * area_tread_contact * (t_tread - t_track)
q_deflection = p_deflection * h_deflection * vx * abs(fztire)
q_air2ambient = p_air2ambient * h_natconv * area_rim * (t_air - t_ambient)
du[1] = der_t_tread = (q_friction - q_tread2carcass - q_tread2road - q_tread2ambient_forcedconv - q_tread2ambient_natconv)/(m_tread * cp_tread)
du[2] = der_t_carcass = (q_tread2carcass + q_deflection - q_carcass2air - q_carcass2ambient_natconv)/(m_carcass * cp_carcass)
du[3] = der_t_air = (q_carcass2air - q_air2ambient)/(m_air * cp_air)
end
function explicit(y,p,t)
sqrt(1-y^2)
end
[slack] <torkel.loman> Is there something I can do to help solve https://discourse.julialang.org/t/strange-maxiters-numeric-instability-occured-when-solving-certain-sde/49392/9 ?
Happy to do some digging, but unsure what I should be looking for.
((((((1.0 * α₁ * (0.0 + (10000.0 * (1.0 - (1.0 / (1.0 + exp(-20.0 * (t - 24.0)))))))) / (α₂ + (0.0 + (10000.0 * (1.0 - (1.0 / (1.0 + exp(-20.0 * (t - 24.0))))))))) - (((1.0 * α₄) * x₁(t)) / (α₅ + x₁(t)))) - ((x₁(t) * α₆₂) / 0.1048)) + (((α₆₅ * x₃(t)) / ((α₆₆ + x₃(t)) * ((0.1048 * ((x₁(t) + x₇(t)) + ((((x₉(t) + x₁₀(t)) + x₁₁(t)) + x₁₂(t)) * 2))) / α₇₁))) / 0.1048)) + ((((x₉(t) + x₁₀(t)) + x₁₁(t)) + x₁₂(t)) * α₂₃)) - ((x₁(t) * x₅(t)) * α₂₄)ERROR: Failed to apply rule ~~(z::_isone) * ~~x => ~x on expression (1.0 * (0.0 + (10000.0 * (1.0 - (1.0 / (1.0 + exp(-20.0 * (t - 24.0)))))))) * α₁This is for sensitivity analysis of ODE's. and this is the only equation in my ODE system that explicitly depends on t, gradients of other equations work fine.
[slack] <Peter J> My code does a lot of parameter-parallel ODE solving, and I'd like to do them on the GPU, but since a lot of julia is not supported on the GPU by DiffEqGPU (broadcast, matrix multiply...) would this idea work?
Given an ode described by f, an IC u_0 and a list of parameters [p_1, p_2,... p_n]. Create a new ode function big_f(du,u,p,t) , and w_0. Where w_0 is a cuarray consisting of u_0 concatenated n times, and big_f applies f seperately to each copy of u_0, each with a different p_i.
ifelse(stillon,1,...)
stillon == false causes the callback to not fire no matter what.
[slack] <Andreas Schlaginhaufen> Is it possible that adjoints methods stop working when the ODE function includes a Zygote.gradient() call? For me this always leads to strange errors like "MethodError: no method matching mul!(::Array{Float64,1}, ::Array{Float64,2}, ::LinearAlgebra.Adjoint{Float64,Array{Float64,1}}, ::Bool, ::Bool)".
Here is a simple example (gradients of ODE function work, but adjoint method is failing):: https://files.slack.com/files-pri/T68168MUP-F01EEJNVDR9/download/mwe
[slack] <Andreas Schlaginhaufen> Is it possible that adjoints methods stop working when the ODE function includes a Zygote.gradient() call? For me this always leads to problems like "MethodError: no method matching mul!(::Array{Float64,1}, ::Array{Float64,2}, ::LinearAlgebra.Adjoint{Float64,Array{Float64,1}}, ::Bool, ::Bool)".
Here is a simple example (gradients of ODE function work, but adjoint method is failing):