tangent_integrator not compatible with stiff solvers

[jamblejoe]

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)

# IIP Tangent Space dynamics
function create_tangent(@nospecialize(f::F), @nospecialize(jacobian::JAC), J::JM,
    ::Val{true}, ::Val{k}) where {F, JAC, JM, k}
    J = deepcopy(J)
    tangentf = (du, u, p, t) -> begin
        uv = @view u[:, 1]
        f(view(du, :, 1), uv, p, t)
        jacobian(J, uv, p, t)
        mul!((@view du[:, 2:(k+1)]), J, (@view u[:, 2:(k+1)]))
        nothing
    end
    return tangentf
end

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.

---

Want to back this issue? Post a bounty on it! We accept bounties via Bountysource.

[Datseris]

https://discourse.julialang.org/t/chaostools-lyapunovs-with-stiff-solvers/35961

[rusandris]

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())