Open jamblejoe opened 5 years ago
My problem is that I have to use a stiff solver. This leaves me two possibilities: using autodifferentiation, which at the moment, doesn't work or don't use autodifferentiation, which works for example with the Lorentz system, but with larger systems with complicated Jacobians (like the one I described in the discourse post) it is not an option for me because it runs practically forever (I waited 8+ hours before I stopped the process). The MWE:
function lorenz(du,u,p,t)
sigma = p[1]
rho = p[2]
beta = p[3]
# eom
du[1] = sigma*(u[2]-u[1])
du[2] = u[1]*(rho-u[3])-u[2]
du[3] = u[1]*u[2]-beta*u[3]
end
function lorenz_jac(J,u,p,t)
sigma = p[1]
rho = p[2]
beta = p[3]
J[1,1] = -sigma
J[1,2] = sigma
J[1,3] = 0
J[2,1] = rho-u[3]
J[2,2] = -1
J[2,3] = -u[1]
J[3,1] = u[2]
J[3,2] = u[1]
J[3,3] = -beta
end
state = rand(3)
param = [16.,45.92,4.]
ds = ContinuousDynamicalSystem(lorenz,state,param,lorenz_jac)
lyapunovs(ds,2000.;Ttr = 2000.,alg = Rosenbrock23())
Hi, as shown here https://discourse.julialang.org/t/chaostools-tangent-integrator-stiff-solvers/26464, the tangent_integrator errors, if one provides an algorithm, which uses ForwardDiff.jl to calculate the Jacobian, etc. This is the case for all stiff solvers. One can avoid using auto-differentiation by passing
autodiff=false
to the algorithm, e.g.Rodas5(autodiff=false)
but we often want this feature!The error is due to caching the Jacobi-matrix in the
create_tangent
function (https://github.com/JuliaDynamics/DynamicalSystemsBase.jl/blob/master/src/dynamicalsystem.jl#L343-L355)Especially line 346:
J = deepcopy(J)
.A workaround is described here http://docs.juliadiffeq.org/latest/basics/faq.html#I-get-Dual-number-errors-when-I-solve-my-ODE-with-Rosenbrock-or-SDIRK-methods.