How to efficiently calculate the Jacobian of the solution of a differential equation with respect to the initial conditions?

General Usage Performance
,
1

I need to numerically perform an integral of the form

\int_{0}^{t_{max}} F \left( \mathbf{u}(t), \frac{\partial \mathbf{u}(t)}{ \partial \mathbf{u}(0)} \right) dt

where F is a certain function, \mathbf{u}(t) is the solution of a simple ODE (it is a 1D Hamiltonian system) and \frac{\partial \mathbf{u}(t)}{ \partial \mathbf{u}(0)} is the Jacobian of \mathbf{u}(t) with respect to the initial condition \mathbf{u}(0). For that reason, I need a fast way of evaluating \mathbf{u}(t) and \frac{\partial \mathbf{u}(t)}{ \partial \mathbf{u}(0)} for many values of t.

The following code allows me to calculate these values, but it isn’t efficient:

using DifferentialEquations, ForwardDiff

function f!(du,u,p,t)
    du[1] = exp(-u[2])*(1-exp(-u[2]))
    du[2] = 2*u[1]
end

function get_trajectory(u0,t_max)
    problem = ODEProblem(f!,u0,(0.0,t_max))
    solve(problem)
end

function get_jacobian(u0,t_max)
    problem(u) = ODEProblem(f!,u,(0.0,t_max))
    t->ForwardDiff.jacobian(u->solve(problem(u))(t),u0)
end

u = get_trajectory([1.0,0.0], 1.0)
u(0.5)

jac = get_jacobian([1.0,0.0], 1.0)
jac(0.5)

It works fine for calculating \mathbf{u}, as, when I call get_trajectory, the differential equation is solved once and then each call of u(t) is fast. On the other hand, each time I call jac(t) the differential equation is solved once again, which makes the call slow. I was wondering if there is a way to solve the differential equation only once and get all the values I need.

On the DifferentialEquations.jl documentation there is a page on Sensitivity Analysis that appears to treat a similar question, but the Jacobian is calculated with respect to parameters. I also found this topic which may have a solution to my problem, but the links containing the references to the solution appear broke to me.

2

Hm, the problem of URL stability. Here is what I found.