# Sensitivities with respect to initial conditions in DifferentialEquations.jl

**URL:** <https://discourse.julialang.org/t/sensitivities-with-respect-to-initial-conditions-in-differentialequations-jl/25555>\
**Category:** Numerics\
**Tags:** diffeq\
**Created:** [June 22, 2019, 9:41pm UTC](https://discourse.julialang.org/t/sensitivities-with-respect-to-initial-conditions-in-differentialequations-jl/25555 "2019-06-22T21:41:54Z")\
**Posts on this page:** 18\
**Page:** 1

<div class="post-metadata">

**Author:** ![rjpower4](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rjpower4/32/9016_2.png) [@rjpower4](https://discourse.julialang.org/u/rjpower4)\
**Post date:** [June 22, 2019, 9:41pm UTC](https://discourse.julialang.org/t/sensitivities-with-respect-to-initial-conditions-in-differentialequations-jl/25555/1 "2019-06-22T21:41:54Z")

</div>

For my research, I perform many differential corrections processes. Often times, for these processes I want the transition matrix, \Phi(t\_1, t\_0), which is obtained by propagating the linearization of the system alongside the nonlinear differential equations such that the following approximation can be made:

\mathbf{\delta x}(t\_1) = \Phi(t\_1, t\_0) \mathbf{\delta x}(t\_0)

where \mathbf{\delta x}(t) is the variation with respect to the reference solution at time t. The linear differential equations are obviously

\dot{\Phi}(t, t\_0) = A(t)\Phi(t, t\_0)

where A(t) is the jacobian of the nonlinear sytem at time t.

Now to the actual question, I see that there is infrastructure in DifferentialEquations.jl for finding sensitivity of the state to the parameters, but is there any ergonomic way to find this transition matrix given the equations of motion and Jacobian without gluing them together in some function? My goal is being able to do this in general for many different systems which only need to provide their EOM and Jacobian definitions. A contrived small example would be

```julia

abstract type AbstractSystem end

struct SystemA <: AbstractSystem
    alpha::Float64
end

struct SystemB <: AbstractSystem
    beta::Float64
end

function eom!(du, u, p::SystemA, t)
    du[1] = u[1] + exp(u[2])
    du[2] = -p.alpha* u[2]
end

function jac!(J, u, p::SystemA, t)
    J[1, 1] = 1
    J[1, 2] = exp(u[2])
    J[2, 1] = 0
    J[2, 2] = -p.alpha
end

function eom!(du, u, p::SystemB, t)
    du[1] = u[2] - u[2]^3 - p.beta* u[1]
    du[2] = u[1] - u[2] - u[1] * u[2]
end

function jac!(J, u, p::SystemB, t)
    J[1,1] = -p.beta
    J[1,2] = 1 - 3 * u[2]^2
    J[2,1] = 1 - u[2]
    J[2,2] = -1 - u[1]
end

function stm(u0, system::T, tspan) where {T <: AbstractSystem}
    # Do something here to get the `sol` as well as the 
    # transition matrices over timespan.
end

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 22, 2019, 10:04pm UTC](https://discourse.julialang.org/t/sensitivities-with-respect-to-initial-conditions-in-differentialequations-jl/25555/2 "2019-06-22T22:04:43Z")

</div>

The best way is to probably follow this [http://docs.juliadiffeq.org/latest/analysis/sensitivity.html#Examples-using-ForwardDiff.jl-1](http://docs.juliadiffeq.org/latest/analysis/sensitivity.html#Examples-using-ForwardDiff.jl-1) except put partials on `u0` instead of `p`.

---

<div class="post-metadata">

**Author:** ![rjpower4](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rjpower4/32/9016_2.png) [@rjpower4](https://discourse.julialang.org/u/rjpower4)\
**Post date:** [June 24, 2019, 9:03pm UTC](https://discourse.julialang.org/t/sensitivities-with-respect-to-initial-conditions-in-differentialequations-jl/25555/3 "2019-06-24T21:03:34Z")

</div>

Ok, thanks. Do you forsee this ever being incorporated into DifferentialEquations or is the use case not common enough?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 24, 2019, 11:35pm UTC](https://discourse.julialang.org/t/sensitivities-with-respect-to-initial-conditions-in-differentialequations-jl/25555/4 "2019-06-24T23:35:30Z")

</div>

it’s in there, but it’s just documented poorly. `adjoint_sensitivities_u0` exists too which gives the sensitivites w.r.t. p and u0, but I need to find time to add docs.

Would you mind adding a part to the ForwardDiff example showing how to calculate the forward sensitivities using dual numbers? And add `adjoint_sensitivities_u0` which just returns the tuple `(du0,dp)`? The signature is here and the same as the other adjoints: [https://github.com/JuliaDiffEq/DiffEqSensitivity.jl/blob/master/src/adjoint\_sensitivity.jl#L325-L336](https://github.com/JuliaDiffEq/DiffEqSensitivity.jl/blob/master/src/adjoint_sensitivity.jl#L325-L336) . Going to catch a flight though but there’s both forward and backwards mode, so hopefully these notes are parsable.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 24, 2019, 11:55pm UTC](https://discourse.julialang.org/t/sensitivities-with-respect-to-initial-conditions-in-differentialequations-jl/25555/5 "2019-06-24T23:55:47Z")

</div>

Or if you can’t add the docs on this, ping me in like a day so I get a reminder

---

<div class="post-metadata">

**Author:** ![rjpower4](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rjpower4/32/9016_2.png) [@rjpower4](https://discourse.julialang.org/u/rjpower4)\
**Post date:** [June 25, 2019, 12:51am UTC](https://discourse.julialang.org/t/sensitivities-with-respect-to-initial-conditions-in-differentialequations-jl/25555/6 "2019-06-25T00:51:29Z")

</div>

Yeah, I will see if I can get an example up and running.

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [June 25, 2019, 3:06am UTC](https://discourse.julialang.org/t/sensitivities-with-respect-to-initial-conditions-in-differentialequations-jl/25555/7 "2019-06-25T03:06:07Z")

</div>

You can do this with jet transport, eg using TaylorIntegration.jl. Cc @PerezHz, @lbenet

---

<div class="post-metadata">

**Author:** ![mleprovost](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mleprovost/32/7166_2.png) [@mleprovost](https://discourse.julialang.org/u/mleprovost)\
**Post date:** [January 17, 2020, 1:36am UTC](https://discourse.julialang.org/t/sensitivities-with-respect-to-initial-conditions-in-differentialequations-jl/25555/8 "2020-01-17T01:36:05Z")

</div>

Hello,

I am interested as well in computing the Jacobian of the solution with respect to initial conditions. i.e the matrix dx(t=t)/dx(t =0).

I have found this example from the documentation [https://docs.juliadiffeq.org/stable/analysis/sensitivity/](https://docs.juliadiffeq.org/stable/analysis/sensitivity/)  
However the sensitivity analysis is only done with respect to the parameter `p`, how can I modify this code to construct the jacobian with respect to u0?

Another question, I will be using interpolated function for the RHS of the ODE. Since Interpolations.jl and ForwardDiff.jl are not compatible, how can I go around this problem?

```julia
function f(du,u,p,t)
  du[1] = dx = p[1]*u[1] - p[2]*u[1]*u[2]
  du[2] = dy = -p[3]*u[2] + u[1]*u[2]
end

u0 = [1.0;1.0]
p = [1.5,1.0,3.0]
prob = ODEForwardSensitivityProblem(f,u0,(0.0,10.0),p)

sol = solve(prob,DP8())
x,dp = extract_local_sensitivities(sol)

```

Best,

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [January 17, 2020, 5:01am UTC](https://discourse.julialang.org/t/sensitivities-with-respect-to-initial-conditions-in-differentialequations-jl/25555/9 "2020-01-17T05:01:33Z")

</div>

It’s shown here: [https://docs.juliadiffeq.org/dev/analysis/sensitivity/#concrete\_solve-Examples-1](https://docs.juliadiffeq.org/dev/analysis/sensitivity/#concrete_solve-Examples-1)

---

<div class="post-metadata">

**Author:** ![rjpower4](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rjpower4/32/9016_2.png) [@rjpower4](https://discourse.julialang.org/u/rjpower4)\
**Post date:** [January 30, 2020, 11:41pm UTC](https://discourse.julialang.org/t/sensitivities-with-respect-to-initial-conditions-in-differentialequations-jl/25555/10 "2020-01-30T23:41:57Z")

</div>

That example provides the sensitivity of a scalar valued function of the states, \frac{\partial f(t\_f, \bar{x}\_f)}{\partial \bar{x}\_0} and \frac{\partial f(t\_f, \bar{x}\_f)}{\partial \bar{p}}, however I can not get it to provide the sensitivities of all of the final states with respect to the initial states, i.e., \frac{\partial \bar{x}\_f}{\partial \bar{x}\_0} – more generally, \frac{\partial \bar{f}(t\_f, \bar{x}\_f)}{\partial \bar{x}\_0}. After reading through all of the Zygote.jl documentation I am unsure as to whether this is possible. However, clearly, to determine the sensitivities of the function f(t\_f, \bar{x}\_f), these must be known by the library, correct?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [January 31, 2020, 12:26am UTC](https://discourse.julialang.org/t/sensitivities-with-respect-to-initial-conditions-in-differentialequations-jl/25555/11 "2020-01-31T00:26:50Z")

</div>

I don’t understand the question. You noticed that using `concrete_solve` where you calculate only the end will give `df(xf)/dx0`?

---

<div class="post-metadata">

**Author:** ![rjpower4](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rjpower4/32/9016_2.png) [@rjpower4](https://discourse.julialang.org/u/rjpower4)\
**Post date:** [January 31, 2020, 12:45am UTC](https://discourse.julialang.org/t/sensitivities-with-respect-to-initial-conditions-in-differentialequations-jl/25555/12 "2020-01-31T00:45:58Z")

</div>

The `Zygote.gradient` function does will not provide the sensitivities for a vector valued f. For example, using the same problem on the linked page:

```julia
using DiffEqSensitivity, OrdinaryDiffEq, Zygote

function fiip(du,u,p,t)
  du[1] = dx = p[1]*u[1] - p[2]*u[1]*u[2]
  du[2] = dy = -p[3]*u[2] + p[4]*u[1]*u[2]
end
p = [1.5,1.0,3.0,1.0]; u0 = [1.0;1.0]
prob = ODEProblem(fiip,u0,(0.0,10.0),p)
sol = concrete_solve(prob,Tsit5())

```

Attempting to get the partials of the final state with respect to the initial state from the `Zygote.gradient` function yields the following error:

```julia
du01,dp1 = Zygote.gradient((u0,p)->last(concrete_solve(prob,Tsit5(),u0,p,saveat=0.1,sensealg=QuadratureAdjoint())),u0,p)
# ERROR: LoadError: Output should be scalar; gradients are not defined for output [1.0337581393337192, 0.9063703701433584]

```

This makes sense given the function is named `gradient` and what I desire is a Jacobian.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [January 31, 2020, 7:34am UTC](https://discourse.julialang.org/t/sensitivities-with-respect-to-initial-conditions-in-differentialequations-jl/25555/13 "2020-01-31T07:34:08Z")

</div>

Instead of gradient use `pullback` to build the whole Jacobian:

```julia
out = Zygote.pullback((u0,p)->Array(concrete_solve(prob,Tsit5(),u0,p,saveat=0.1,sensealg=QuadratureAdjoint()))[:,end],u0,p)
[out[2]([1.0,0.0])[2] out[2]([0.0,1.0])[2]]

4×2 Array{Float64,2}:
 2.19418 -6.20335 
 0.199812 -0.703454
 0.575124 -1.70302 
 0.946907 -2.71072 

```

---

<div class="post-metadata">

**Author:** ![bkalpert](https://avatars.discourse-cdn.com/v4/letter/b/e36b37/32.png) [@bkalpert](https://discourse.julialang.org/u/bkalpert)\
**Post date:** [May 26, 2021, 12:29am UTC](https://discourse.julialang.org/t/sensitivities-with-respect-to-initial-conditions-in-differentialequations-jl/25555/14 "2021-05-26T00:29:01Z")

</div>

Was this problem addressed satisfactorily? I.e., Jacobian w.r.t. initial condition u0, at each time step, in addition to that for parameter p, as provided by ODEForwardSensitivityProblem. I’m looking at the suggested documentation, but missing it. I do see how to obtain it solely for the integration endpoint.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [May 26, 2021, 1:20am UTC](https://discourse.julialang.org/t/sensitivities-with-respect-to-initial-conditions-in-differentialequations-jl/25555/15 "2021-05-26T01:20:09Z")

</div>

What’s the question?

---

<div class="post-metadata">

**Author:** ![bkalpert](https://avatars.discourse-cdn.com/v4/letter/b/e36b37/32.png) [@bkalpert](https://discourse.julialang.org/u/bkalpert)\
**Post date:** [May 26, 2021, 1:29am UTC](https://discourse.julialang.org/t/sensitivities-with-respect-to-initial-conditions-in-differentialequations-jl/25555/16 "2021-05-26T01:29:06Z")

</div>

Can one obtain the derivative of the ODE solution, at each time step, w.r.t. both the initial condition u0 and the parameter p?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [May 26, 2021, 10:26am UTC](https://discourse.julialang.org/t/sensitivities-with-respect-to-initial-conditions-in-differentialequations-jl/25555/17 "2021-05-26T10:26:24Z")

</div>

I’d just do it the same old way as always:

```julia
function f(theta)
  _prob = remake(prob,p=theta[1:n],u0=[n+1:end])
  solve(_prob,alg)
end
ForwardDiff.jacobian(f,x)

```

---

<div class="post-metadata">

**Author:** ![bkalpert](https://avatars.discourse-cdn.com/v4/letter/b/e36b37/32.png) [@bkalpert](https://discourse.julialang.org/u/bkalpert)\
**Post date:** [May 26, 2021, 10:36am UTC](https://discourse.julialang.org/t/sensitivities-with-respect-to-initial-conditions-in-differentialequations-jl/25555/18 "2021-05-26T10:36:11Z")

</div>

Thank you for the reply.