nutopy.ivp.djexp

djexp(d2f, tf, tfd, t0, t0d, x0, x0d, dx0, dx0d, pars=None, parsd=None, *, options=None, time_steps=None, args=(), kwargs={})[source]

The djexp function computes a numerical approximation of the derivative of the solution of the linearized initial value problem

\[\begin{split}\dot{x} &= f(t, x, pars), \\ \dot{y} &= \frac{\partial f}{\partial x}(t, x, pars) \cdot y, \\ x(t_0) &= x_0, \\ y(t_0) &= y_0,\end{split}\]

where \(pars \in \mathbf{R}^k\) is a vector of parameters, \(t\) is the time, \(x\) is the state variable and \(y\) its linearization. The time \(t_0\) is refered as the initial time, \(x_0\) as the initial condition, \(y_0\) as the ìnitial condition` of the linear part and \(f\) as the dynamics.

The solution at the final time \(t_f\) is then denoted by \(x(t_f, t_0, x_0, pars)\), \(y(t_f, t_0, x_0, y_0, pars)\).

Hence, djexp approximates

\[\begin{split}x'(t_f, t_0, x_0, pars) \cdot (\delta t_f, \delta t_0, \delta x_0, \delta pars) =& \phantom{+} \frac{\partial x}{\partial t_f}(t_f, t_0, x_0, pars) \cdot \delta t_f \\ & + \frac{\partial x}{\partial t_0}(t_f, t_0, x_0, pars) \cdot \delta t_0 \\ & + \frac{\partial x}{\partial x_0}(t_f, t_0, x_0, pars) \cdot \delta x_0 \\ & + \frac{\partial x}{\partial pars}(t_f, t_0, x_0, pars) \cdot \delta pars\end{split}\]

and

\[\begin{split}y'(t_f, t_0, x_0, y_0, pars) \cdot (\delta t_f, \delta t_0, \delta x_0, \delta y_0, \delta pars) =& \phantom{+} \frac{\partial y}{\partial t_f}(t_f, t_0, x_0, pars) \cdot \delta t_f \\ & + \frac{\partial y}{\partial t_0}(t_f, t_0, x_0, pars) \cdot \delta t_0 \\ & + \frac{\partial y}{\partial x_0}(t_f, t_0, x_0, pars) \cdot \delta x_0 \\ & + \frac{\partial y}{\partial y_0}(t_f, t_0, x_0, pars) \cdot \delta y_0 \\ & + \frac{\partial y}{\partial pars}(t_f, t_0, x_0, pars) \cdot \delta pars.\end{split}\]

Let us introduce \(\delta x(t) = x'(t, t_0, x_0, pars) \cdot (0, \delta t_0, \delta x_0, \delta pars)\) and \(\delta y(t) = y'(t, t_0, x_0, y_0, pars) \cdot (0, \delta t_0, \delta x_0, \delta y_0, \delta pars)\). Then, \(\delta x(t_f)\) and \(\delta y(t_f)\) are computed solving the following variational equations (we omit \((t, x, pars)\)):

\[\begin{split}\dot{x} &= f, \\ \dot{\delta x} &= \frac{\partial f}{\partial x} \cdot \delta x + \frac{\partial f}{\partial pars} \cdot \delta pars, \\ \dot{y} &= \frac{\partial f}{\partial x} \cdot y, \\ \dot{\delta y} &= \frac{\partial^2 f}{\partial x^2} \cdot (y, \delta x) + \frac{\partial^2 f}{\partial pars \partial x} \cdot (y, \delta pars) + \frac{\partial f}{\partial x} \cdot \delta y, \\ x(t_0) &= x_0, \\ \delta x(t_0) &= \delta x_0 - \delta t_0\, f(t_0, x_0, pars), \\ y(t_0) &= y_0, \\ \delta y(t_0) &= \delta y_0 - \delta t_0\, \left( \frac{\partial f}{\partial x}(t_0, x_0, pars) \cdot y_0 \right).\end{split}\]

The missing parts are given by:

\[\frac{\partial x}{\partial t_f}(t_f, t_0, x_0, pars) \cdot \delta t_f = \delta t_f\, f(t_f, x(t_f), pars)\]

and

\[\frac{\partial y}{\partial t_f}(t_f, t_0, x_0, y_0, pars) \cdot \delta t_f = \delta t_f\, \left( \frac{\partial f}{\partial x}(t_f, x(t_f), pars) \cdot y(t_f) \right)\]

where \(x(t_f)\) stands for \(x(t_f, t_0, x_0, pars)\) and \(y(t_f)\) stands for \(y(t_f, t_0, x_0, y_0, pars)\).

Parameters

  • d2f (callable) –

    A vector function

    z, dz, d2z = d2f(t, x, dx1, dx2) if (pars, parsd) is None, where

    z = f(t, x), dz = df/dx(t, x) . dx1 and

    d2z = d2f/dx2(t, x) . (dx1, dx2)

    or

    z, dz, d2z = d2f(t, x, dx1, dx2, pars, parsd) if (pars, parsd) is not None, where

    z = f(t, x, pars),

    dz = df/dx(t, x, pars) . dx1 + df/dpars(t, x, pars) . parsd and

    d2z = d2f/dx2(t, x, pars) . (dx1, dx2) + d2f/dparsdx(t, x, pars) . (parsd, dx2)

  • tf (float) – Final time

  • tfd (float) – Increment \(\delta t_f\)

  • t0 (float) – Initial time

  • t0d (float) – Increment \(\delta t_0\)

  • x0 (float vector) – Initial condition

  • x0d (float vector) – Increment \(\delta x_0\)

  • dx0 (float vector or matrix) – Initial condition \(y_0\) of the linear part of the equation

  • dx0d (float vector or matrix) – Initial increment \(\delta y_0\) of the linear part of the equation

  • pars (float vector, optional, positional) – Parameters

  • parsd (float vector, optional, positional) – Increment \(\delta pars\)

  • options (nutopy.ivp.Options, optional, key only) – A dictionnary of solver options. The integrator scheme is chosen setting SolverMethod.

  • time_steps (float vector, optional, key only) –

    If time_steps is not None, then it is a vector of time steps, that is sol.tout = time_steps. If time_steps is None, then sol.tout = [t0, tf] if a fixed-step integrator scheme is used while sol.tout is given by the integrator if a variable step-size scheme is used.

    Remark. Note that in the case of a variable step-size integrator scheme, time_steps is used only for output, that is the steps are computed by the integrator but dense output is used to provide by interpolation the solution x at times in time_steps.

  • args (tuple, key only) – optional arguments for d2f

  • kwargs (dictionnary, key only) – keywords optional arguments for d2f

Returns

sol – A dictionary/struct of outputs:

  • sol.xf or sol.get('xf') : final point x(tf, t0, x0, pars);

  • sol.xfd : derivative x'(tf, t0, x0, pars) . (tfd, t0d, x0d, parsd);

  • sol.dxf : final point y(tf, t0, x0, dx0, pars);

  • sol.dxfd : derivative y'(tf, t0, x0, dx0, pars) . (tfd, t0d, x0d, dx0d, parsd);

  • sol.tout : the integration time-steps: tout=time_steps if provided;

  • sol.xout : the integration points: xout[i] = x(tout[i], t0, x0, pars);

  • sol.xdout : derivative at integration time-steps: xdout[i] = x'(tout[i], t0, x0, pars) . (tfd, t0d, x0d, parsd);

  • sol.dxout : the integration points: dxout[i] = y(tout[i], t0, x0, dx0, pars);

  • sol.dxdout : derivative at integration time-steps: dxdout[i] = y'(tout[i], t0, x0, dx0, pars) . (tfd, t0d, x0d, dx0d, parsd);

  • sol.success : whether or not the integrator exited successfully;

  • sol.status : termination status of the integrator;

  • sol.message : description of the cause of the termination;

  • sol.nfev : number of evaluations of the dynamics;

  • sol.nsteps : number of integration time-steps.

See also

nutopy.ivp.exp, nutopy.ivp.dexp, nutopy.ivp.jexp, nutopy.ivp.Options, nutopy.errors

Notes

Examples

Example: Second derivative w.r.t to the initial condition d^2x/dx0^2

>>> import numpy as np
>>> from nutopy import ivp

>>> def d2f(t, x, dx1, dx2):
...     y       = t / np.sin(x)
...    dfdx    = - t * np.cos(x) / (np.sin(x) ** 2)
...    dy      = dfdx * dx1
...    d2fdx2  = t  * (np.sin(x) ** 2 + 2.0 * np.cos(x) ** 2) / (np.sin(x) ** 3)
...    d2y     = d2fdx2 * dx1 * dx2
...    return y, dy, d2y

>>> t0   = 0.1
>>> tf   = np.sqrt(1.0 + 2.0 * np.cos(np.pi / 2.0) + 0.0 ** 2)
>>> x0   = np.pi / 2.0
>>> dx0  = 1.0
>>> tfd  = 0.0
>>> t0d  = 0.0
>>> x0d  = 1.0
>>> dx0d = 0.0
>>> sol = ivp.djexp(d2f,tf, tfd, t0, t0d, x0, x0d, dx0, dx0d)
>>> print('d^2x/dx0^2(tf, t0, x0) = ', sol.dxfd)
d^2x/dx0^2(tf, t0, x0) =  0.7545817755798825
Raises

previous

nutopy.ivp.jexp

next

nutopy.nle