title: Time components of numdifftools.Jacobian should not be equal to zero
body: <p>I am solving differential equations numerically using <code>solve_ivp</code> in order to integrate positions and velocities from an initial condition. This is straighfoward in Python and easy to do. But I am also calculating the partials of final positions and velocities with respect to the initial positions/velocities and time using <code>numdifftools.Jacobian</code>. This I have more trouble mastering. Let me explain.</p>
  <p>I have a set of differential equations (see <code>rhs_CR3BP</code> function) that I numerically integrate over a timespan <code>[0,T]</code> using the Python solver <code>solve_ivp</code> given an initial state vector <code>X_0</code> at t=0 (position and velocity). At the end of the integration, I have all the state vectors from t=0 to t=T stored in <code>solve_ivp.y</code>, with the final state vector <code>X_f</code> at t=T stored in <code>solve_ivp.y[:, -1]</code>.</p>
  <p>I am trying to compute numerically partial derivatives of the final state vector <code>X_f</code> with respect to the initial state vector <code>X_0</code> and period <code>T</code>, which are in fact the components of a Jacobian matrix.</p>
  <p>Consider the following vector:</p>
  <pre><code>X = [x_0, y_0, z_0, vx_0, vy_0, vz_0, T] # this is X_0 and T
  </code></pre>
  <p>The vectorial partial derivative of <code>X_f</code> relative to <code>X</code> is given by the following Jacobian matrix:</p>
  <p><a href="https://i.stack.imgur.com/SXZVu.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/SXZVu.png" alt="enter image description here" /></a></p>
  <p>To compute this matrix, I made a function (see <code>propagation</code> function) that takes as a variable a vector composed of <code>X_0</code> and <code>T</code> called <code>X</code>, and returns <code>X_f</code>. Then, this function is combined with <code>numdifftools.Jacobian</code> to compute the Jacobian matrix. When I run my code (see below), the Jacobian matrix is properly calculated for the partial derivatives of <code>X_f</code> relative to <code>X_0</code> (first 6 columns), but the very last column with partial derivatives of <code>X_f</code> relative to <code>T</code> are equal to zero. I know this cannot be because those derivatives are actually the final velocities and accelerations at t=T.</p>
  <p><a href="https://i.stack.imgur.com/Uvpcz.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/Uvpcz.png" alt="enter image description here" /></a></p>
  <p>How come the initial condition <code>X_0</code> is taken into account in the evolution of <code>X_f</code> and the time <code>T</code> is not? That does not make any sense. In <code>solve_ivp</code>, it looks as if the timespan (<code>tspan=(0,T)</code>) is completely ignored but the initial condition (<code>y0=X_0</code>) is not. Why is that? Am I not using <code>solve_ivp</code> and <code>numdifftools.Jacobian</code> properly?</p>
  <p>If my question needs more details, please let me know, I'll happily provide more explanations by the end of day.</p>
  <pre><code>import numpy as np
  from scipy.integrate import solve_ivp
  import numdifftools as nd
  
  def rhs_CR3BP(t, X0, mu):
      &quot;&quot;&quot;Integrates the CR3BP equations of motion&quot;&quot;&quot;
      
      x, y, z, v_x, v_y, v_z = X0
      
      r = np.sqrt((x-1+mu)**2+y**2+z**2)
      d = np.sqrt((x+mu)**2+y**2+z**2)
      
      a_x = 2*v_y + x - (1-mu)*(x+mu)/d**3 - mu/r**3*(x-1+mu)
      a_y = -2*v_x + y - (1-mu)*y/d**3 - mu*y/r**3
      a_z = -(1-mu)*z/d**3 - mu*z/r**3
      
      return np.array([v_x, v_y, v_z, a_x, a_y, a_z])
  
  def propagation(X, mu):
      
      X_0 = X[0:6]
      T = X[6].real
      
      sol = solve_ivp(fun=rhs_CR3BP, t_span=(0, T), y0=X_0, method='RK45', rtol=1e-10, atol=1e-10, args=(mu,), dense_output=True)
      X_f = sol.y[:, -1]
      
      return X_f
  
  T = 0.732 # period
  x, y, z, vx, vy, vz = 1.085, 0, 0, 0, -0.464, 0 # position and velocity components
  X0 = [x, y, z, vx, vy, vz] # initial state vector
  X = np.array([x, y, z, vx, vy, vz, T]) # X_0 and T
  
  mu = 0.012 # arbitrary parameter
  f = lambda X: propagation(X, mu)
  f = nd.Jacobian(f, method='complex')
  DF = np.real(f(X))
  </code></pre>