Dear Julia users,
I have a complex ODE systems in which I use Gaussian overlaps, and I obtain wrong results when I estimate the Jaccobian of my ODE system using ForwardDiff.jacobian. I get a diagonal matrix, with zero for all values out of the diagonal, which is obviously wrong.
Here is a reproducible example:
using ForwardDiff
using Distributions
using QuadGK
using OrdinaryDiffEq
ns=10
epsilon= 0.1
p = epsilon,ns
function LV_model2(u,p,t)
epsilon,ns = p
du=similar(u)
#loop to define derivatives
for i in 1:ns
#partial derivatives that will be use to calculate the derivative
function pop_derivative_p(x)
overlaps = Vector{}(undef, ns)
for v in 1:ns
mu1=u[v]
function inte(theta)
t1= pdf.(Normal.(x,0.1),theta)
t2= pdf.(Normal.(mu1,0.1),theta)
return(min(t1,t2))
end
overlaps[v] = quadgk(inte, -10, 10, rtol=1e-6)[1] #this is the part that is the problem, if I use inte(x) instead it seems it works
end
d_pop= 1 .- (sum(overlaps))
return d_pop
end
∂f_∂x_p(x) = ForwardDiff.derivative(x->pop_derivative_p(x),x)
#### Final derivatives
du[i]= epsilon * ∂f_∂x_p(u[i])
end
return du
end
uinit=rand(ns)
jacob=ForwardDiff.jacobian(u -> LV_model2(u,p,nothing),uinit)
Here is the jacobian I get:
It seems that the problem comes from the calculation of overlaps, as if I simply replace overlaps[v] = quadgk(inte, -10, 10, rtol=1e-6)[1] by overlaps[v] = inte(x) I obtain a Jacobian that makes sense.
I tried to calculate the overlap manually, using sum(inte.([-10:0.01:10;]))*0.01 but I get the same problem as when using quadgk
I don’t understand what causes this wrong result. Did you already encounter that problem? Would you have a solution?