"""
ivp module description
"""
import numpy as np
from .mod_ivpsolve import mod_ivpsolve
from .errors import *
from . import options
def directionalFiniteDiff(fun, t, x, dx, pars):
veps = np.finfo(float).eps
if(isinstance(x,float)):
h = np.sqrt(veps)*np.sqrt(np.max(1.0, np.abs(x)))
j = (fun(t, x+h*dx, pars)-fun(t, x, pars))/h
else:
n = x.size
j = np.zeros((n,1))
f = fun(t, x, pars)
h = np.sqrt(veps)*np.sqrt(np.maximum(1.0, np.abs(np.dot(x, dx))))
j[:,0] = (fun(t, x+h*dx, pars)-f)/h # on laisse dans ce sens car c pour fortran
return j
# -----------------------------------------------------------------------------------
# -----------------------------------------------------------------------------------
[docs]def exp(f, tf, t0, x0, pars=None, *, options=None, time_steps=None, dfdx=None, \
args=(), kwargs={}):
r"""
The ``exp`` function computes a numerical approximation of the solution of the initial
value problem
.. math::
\frac{d}{d t} x = \dot{x} = f(t, x, pars), \quad x(t_0) = x_0,
where :math:`pars \in \mathbf{R}^k` is a vector of ``parameters``, :math:`t` is the
``time`` and :math:`x` is the ``state`` variable. The time :math:`t_0` is refered as
the ``initial time``, :math:`x_0` as the ``initial condition`` and :math:`f` as
the ``dynamics``. The solution of the initial value problem at the ``final time``
:math:`t_f` is denoted :math:`x(t_f, t_0, x_0, pars)`.
Parameters
----------
f : callable
A vector function representing the dynamics: ``f(t, x)`` or ``f(t, x, pars)``
tf : :class:`float`
Final time
t0 : :class:`float`
Initial time
x0 : float vector
Initial condition
pars : float vector, optional, positional
If ``pars`` is ``None`` then provide ``f(t,x)`` else ``f(t,x,pars)``.
options : :meth:`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``.
dfdx : callable, optional, key only
a vector function ``dfdx(t, x, dx)``
or ``dfdx(t, x, dx, pars)`` if ``pars`` is not ``None``, computing
.. math::
\frac{\partial f}{\partial x}(t, x, pars) \cdot dx,
where ``dx`` is a float vector of the same dimension as ``x``.
**Remark.** ``dfdx`` is used only if ``SolverMethod`` is an implicit Runge-Kutta
scheme.
args : tuple, key only
optional arguments for ``f`` (and ``dfdx`` if provided)
kwargs : dictionnary, key only
keywords optional arguments for ``f`` (and ``dfdx`` if provided)
Returns
-------
sol
A dictionary/struct of outputs:
- ``sol.xf`` or ``sol.get('xf')`` : final point ``x(tf, t0, x0, pars)``;
- ``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.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.dexp
nutopy.ivp.jexp
nutopy.ivp.djexp
nutopy.ivp.Options
nutopy.errors
Notes
-----
Use :class:`numpy.ndarray` to define float vectors.
Examples
--------
.. plot::
:format: doctest
:include-source: True
>>> import matplotlib.pyplot as plt
>>> import numpy as np
>>> from nutopy import ivp
>>> def f(t, x):
... dx = np.zeros(x.size)
... dx[0] = -x[0] + x[1]
... dx[1] = x[1]
... return dx
>>> x0 = np.array([-1.0, 1.0])
>>> t0 = 0.0
>>> tf = 1.0
>>> sol = ivp.exp(f, tf, t0, x0)
>>> print(sol.xf)
[0.80732175 2.71828183]
>>> plt.plot(sol.tout, sol.xout[:,0], 'b', sol.tout, sol.xout[:,1], 'r')
>>> plt.show()
Raises
------
ArgumentTypeError
Variable tf must be a float or cast to float.
ArgumentTypeError
Variable t0 must be a float or cast to float.
ArgumentTypeError
Variable options must be a Options object.
ArgumentDimensionError
Variable x0 must be a one dimensional array.
ArgumentDimensionError
Variable pars must be a one dimensional array.
InputArgumentError
You must have time_steps[0]=t0 and time_steps[-1]=tf.
InputArgumentError
The argument time_steps cannot be of size 1.
"""
# -----------------------------------------------------------------------------------
#
# None
#
# -----------------------------------------------------------------------------------
pars_is_none=False
if pars is None:
pars_is_none=True # in this case, we should have f(t,x) and not f(t,x,pars)
pars = np.array([])
if time_steps is None:
time_steps = np.array([])
# -----------------------------------------------------------------------------------
#
# Cast of input variables
#
# -----------------------------------------------------------------------------------
if(not isinstance(t0, float)):
try:
t0 = float(t0)
except Exception as exc:
raise ArgumentTypeError('Variable t0 must be a float or cast to float.')
if(not isinstance(tf, float)):
try:
tf = float(tf)
except Exception as exc:
raise ArgumentTypeError('Variable tf must be a float or cast to float.')
if(not isinstance(x0, float)):
x0 = np.asarray(x0)
if(not isinstance(time_steps, float)):
time_steps = np.asarray(time_steps)
else:
msgErreur = 'The argument time_steps cannot be of size 1.';
raise InputArgumentError(msgErreur)
if(not isinstance(pars, float)):
pars = np.asarray(pars)
# -----------------------------------------------------------------------------------
#
# Tests on input variables
#
# -----------------------------------------------------------------------------------
# -----------------------------------------------------------------------------------
if options==None:
options = Options()
# -----------------------------------------------------------------------------------
#Variable x0 must be a vector or a float
msgErreur = 'Variable x0 must be a one dimensional numpy array.';
if(isinstance(x0,np.ndarray)):
#Variable x0 must be a vector
if(x0.ndim==1):
n = x0.size
elif(x0.ndim==2):
(lig,col) = x0.shape
if(lig==1 or col==1):
n = x0.size
else :
raise ArgumentDimensionError(msgErreur)
else :
raise ArgumentDimensionError(msgErreur)
elif(isinstance(x0,float)):
n = 1
#else :
#raise ArgumentTypeError(msgErreur)
# -----------------------------------------------------------------------------------
#Variable pars must be a vector or a float
msgErreur = 'Variable pars must be a one dimensional numpy array.';
if(isinstance(pars,np.ndarray)):
#Variable pars must be a vector
if(pars.ndim==1):
npars = pars.size
elif(pars.ndim==2):
(lig,col) = pars.shape
if(lig==1 or col==1):
npars = pars.size
else :
raise ArgumentDimensionError(msgErreur)
else :
raise ArgumentDimensionError(msgErreur)
elif(isinstance(pars,float)):
npars = 1
#else :
# raise ArgumentTypeError(msgErreur)
# -----------------------------------------------------------------------------------
#Variable options must be a Options object
msgErreur = 'Variable options must be a Options object.';
if(not isinstance(options,Options)):
raise ArgumentTypeError(msgErreur)
# -----------------------------------------------------------------------------------
#
# prepare call to Fortran solver
#
# -----------------------------------------------------------------------------------
# -----------------------------------------------------------------------------------
[dw, iw, lsw, sw] = options.getAllOptions()
ndw = len(dw)
niw = len(iw)
nsw = len(lsw)
nswint = len(sw)
swint = [ord(c) for c in sw] #Code ascii de la chaine de caractere
# -----------------------------------------------------------------------------------
if(not isinstance(x0,float)):
x0.astype(dtype='d',order='F')
if(not isinstance(pars,float)):
pars.astype(dtype='d',order='F')
# -----------------------------------------------------------------------------------
if time_steps.size==0:
tspan = np.array([t0, tf])
elif time_steps.size==1:
msgErreur = 'The argument time_steps cannot be of size 1.';
raise InputArgumentError(msgErreur)
else:
if(t0==time_steps[0]) and (tf==time_steps[-1]):
tspan = time_steps
else:
msgErreur = 'You must have time_steps[0]=t0 and time_steps[-1]=tf.';
raise InputArgumentError(msgErreur)
nt = tspan.size
tspan.astype(dtype='d',order='F')
# -----------------------------------------------------------------------------------
# args
if not isinstance(args,tuple):
args = (args,)
# f
if pars_is_none:
myf = lambda t, x, pars: f(t, x, *args, **kwargs)
else:
myf = lambda t, x, pars: f(t, x, pars, *args, **kwargs)
# -----------------------------------------------------------------------------------
# df
if dfdx is None:
mydfdx = lambda t, x, dx, pars: directionalFiniteDiff(myf, t, x, dx, pars)
else:
if pars_is_none:
mydfdx = lambda t, x, dx, pars: dfdx(t, x, dx, *args, **kwargs)
else:
mydfdx = lambda t, x, dx, pars: dfdx(t, x, dx, pars, *args, **kwargs)
# -----------------------------------------------------------------------------------
[flag, nfev, nsteps, messageFortran] = mod_ivpsolve.expfunpy(myf,mydfdx,x0,pars, \
tspan,dw,iw,swint,lsw,n,npars,nt,ndw,niw,nsw,nswint)
message = messageFortran.strip().decode()
# -----------------------------------------------------------------------------------
#
# post computations, prepare outputs
#
# -----------------------------------------------------------------------------------
# -----------------------------------------------------------------------------------
tout = mod_ivpsolve.tout
xout_for = mod_ivpsolve.xout
# tensorize
if(not isinstance(x0,float)):
x0shape = x0.shape
if time_steps.size==2:
tout = np.array([tout[0], tout[-1]])
if(not isinstance(x0,float)):
x0_out = np.reshape(xout_for[:, 0], x0shape)
xf_out = np.reshape(xout_for[:, -1], x0shape)
else:
x0_out = xout_for[0, 0]
xf_out = xout_for[0, -1]
xout = np.concatenate(( x0_out[None], xf_out[None] ), axis=0)
else:
if(not isinstance(x0,float)):
xout = np.reshape(xout_for[:, 0], x0shape)[None]
else:
xout = xout_for[0, 0][None]
for i in range(nsteps-1):
if(not isinstance(x0,float)):
x = np.reshape(xout_for[:, i+1], x0shape)
else:
x = xout_for[0, i+1]
xout = np.concatenate((xout, x[None]), axis=0)
#
xf = xout[-1]
status = flag
if(status==1 or status==0):
success=True
else:
success=False
# -----------------------------------------------------------------------------------
sol = dict(xf=xf, tout=tout, xout=xout, status=status, success=success, \
message=message, nfev=nfev, nsteps=nsteps)
expsol = ExpSol(sol)
return expsol
# -----------------------------------------------------------------------------------
# -----------------------------------------------------------------------------------
class ExpSol:
def __init__(self, *args, **kwargs):
"""
Doc de __init__
"""
#args -- tuple of anonymous arguments
#kwargs -- dictionary of named arguments
self.sol = dict()
self.update(*args, **kwargs)
self.xf = self.sol['xf']
self.tout = self.sol['tout']
self.xout = self.sol['xout']
self.status = self.sol['status']
self.success = self.sol['success']
self.message = self.sol['message']
self.nfev = self.sol['nfev']
self.nsteps = self.sol['nsteps']
def update(self, *args, **kwargs):
"""
Doc de update
"""
for dico in args:
if(isinstance(dico,dict)):
self.__private_update(dico)
else:
raise ArgumentTypeError('ExpSol requires arguments of type dict or ' + \
'a list of param=value, see help(ExpSol)!')
self.__private_update(kwargs)
def __private_update(self, dictionary):
"""
Doc de __private_update
"""
for key, val in dictionary.items():
if key in self.sol:
if(isinstance(self.sol[key],val.__class__)):
self.sol[key] = val
else:
raise ArgumentTypeError('ExpSol ' + key + ' must be of type: ' + \
str(self.sol[key].__class__) + ' instead of type: ' + \
str(val.__class__) + '!')
else:
self.sol[key] = val
def get(self, key):
"""
Doc de get
"""
try:
value = self.sol[key]
except KeyError as err:
raise OptionNameError(key + ' is not a ExpSol, see help(ExpSol)!\n')
return value
def __str__(self):
"""
Doc de __str__
"""
keys = sorted(self.sol.keys())
string = self.__class__.__name__ + ':\n\n'
for key in keys:
string = string + "\t" + key.ljust(20) + ' = ' + str(self.sol[key]) + '\n'
return string
# -----------------------------------------------------------------------------------
# -----------------------------------------------------------------------------------
[docs]def dexp(df, tf, tfd, t0, t0d, x0, x0d, pars=None, parsd=None, *, \
options=None, time_steps=None, args=(), kwargs={}):
r"""
The ``dexp`` function computes a numerical approximation of the derivative of the
solution of the initial value problem
.. math::
\frac{d}{d t} x = \dot{x} = f(t, x, pars), \quad x(t_0) = x_0,
where :math:`pars \in \mathbf{R}^k` is a vector of ``parameters``, :math:`t` is the
``time`` and :math:`x` is the ``state`` variable. The time :math:`t_0` is refered as
the ``initial time``, :math:`x_0` as the ``initial condition`` and :math:`f` as the
``dynamics``. The solution of the initial value problem at the ``final time``
:math:`t_f` is denoted :math:`x(t_f, t_0, x_0, pars)`.
Hence, ``dexp`` approximates
.. math::
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.
Let us introduce :math:`\delta x(t) = x'(t, t_0, x_0, pars) \cdot
(0, \delta t_0, \delta x_0, \delta pars)`.
Then, :math:`\delta x(t_f)` is computed solving the following `variational equations`:
.. math::
\dot{x} &= f(t, x, pars), \\
\dot{\delta x} &= \frac{\partial f}{\partial x}(t, x, pars) \cdot \delta x +
\frac{\partial f}{\partial pars}(t, x, pars) \cdot \delta pars, \\
x(t_0) &= x_0, \\
\delta x(t_0) &= \delta x_0 - f(t_0, x_0, pars)\, \delta t_0.
The missing part is given by:
.. math::
\frac{\partial x}{\partial t_f}(t_f, t_0, x_0, pars) \cdot \delta t_f =
f(t_f, x(t_f, t_0, x_0, pars), pars)\, \delta t_f.
Parameters
----------
df : callable
A vector function
``y, yd = df(t, x, xd)`` if ``(pars, parsd)`` is ``None``, where
``y = f(t, x)`` and ``yd = df/dx(t, x) . xd``
or
``y, yd = df(t, x, xd, pars, parsd)`` if ``(pars, parsd)`` is not
``None``, where ``y = f(t, x, pars)`` and
``yd = df/dx(t, x, pars) . xd + df/dpars(t, x, pars) . parsd``.
tf : :class:`float`
Final time
tfd : :class:`float`
Increment :math:`\delta t_f`
t0 : :class:`float`
Initial time
t0d : :class:`float`
Increment :math:`\delta t_0`
x0 : float vector
Initial condition
x0d : float vector
Increment :math:`\delta x_0`
pars : float vector, optional, positional
Parameters
parsd : float vector, optional, positional
Increment :math:`\delta pars`
options : :meth:`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 ``df``
kwargs : dictionnary, key only
keywords optional arguments for ``df``
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.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.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.jexp
nutopy.ivp.djexp
nutopy.ivp.Options
nutopy.errors
Notes
-----
- Use :class:`numpy.ndarray` to define float vectors.
- :meth:`nutopy.ivp.dexp` is the derivative of :meth:`nutopy.ivp.exp`
Examples
--------
Example: dx/dx2, x = (x1, x2)
>>> import numpy as np
>>> from nutopy import ivp
>>> def df(t, x, xd):
... y = np.zeros(x.size)
... y[0] = -x[0] + x[1]
... y[1] = x[1]
... yd = np.zeros(x.size)
... yd[0] = -xd[0] + xd[1]
... yd[1] = xd[1]
... return (y, yd)
>>> tf = 1.0
>>> t0 = 0.0
>>> x0 = np.array([-1.5, 0.5])
>>> tfd = 0.0
>>> t0d = 0.0
>>> x0d = np.array([0.0, 1.0])
>>> sol = ivp.dexp(df, tf, tfd, t0, t0d, x0, x0d)
>>> print('dx/dx2(tf, t0, x0) = ', sol.xfd)
dx/dx2(tf, t0, x0) = [1.17520119 2.71828183]
Raises
------
ArgumentTypeError
Variable tf must be a float or cast to float.
ArgumentTypeError
Variable tfd must be a float or cast to float.
ArgumentTypeError
Variable t0 must be a float or cast to float.
ArgumentTypeError
Variable t0d must be a float or cast to float.
ArgumentTypeError
Variable options must be a Options object.
ArgumentDimensionError
Variable x0 must be a one dimensional array.
ArgumentDimensionError
Variable x0d must be a one dimensional array.
ArgumentDimensionError
Variable x0d must have the same length than x0.
ArgumentDimensionError
Variable pars must be a one dimensional array.
ArgumentDimensionError
Variable parsd must be a one dimensional array.
ArgumentDimensionError
Variable parsd must have the same length than pars.
InputArgumentError
Either give pars and parsd or any.
InputArgumentError
You must have time_steps[0]=t0 and time_steps[-1]=tf.
InputArgumentError
The argument time_steps cannot be of size 1.
"""
# -----------------------------------------------------------------------------------
#
# None
#
# -----------------------------------------------------------------------------------
pars_is_none=False
if pars is None:
pars_is_none=True # in this case, we should have f(t,x) and not f(t,x,pars)
pars = np.array([])
parsd_is_none=False
if parsd is None:
parsd_is_none=True
parsd = np.array([])
if pars_is_none != parsd_is_none:
raise InputArgumentError('Either give pars and parsd or any.')
if time_steps is None:
time_steps = np.array([])
# -----------------------------------------------------------------------------------
#
# Cast of input variables
#
# -----------------------------------------------------------------------------------
if(not isinstance(t0, float)):
try:
t0 = float(t0)
except Exception as exc:
raise ArgumentTypeError('Variable t0 must be a float or cast to float.')
if(not isinstance(tf, float)):
try:
tf = float(tf)
except Exception as exc:
raise ArgumentTypeError('Variable tf must be a float or cast to float.')
if(not isinstance(t0d, float)):
try:
t0d = float(t0d)
except Exception as exc:
raise ArgumentTypeError('Variable t0d must be a float or cast to float.')
if(not isinstance(tfd, float)):
try:
tfd = float(tfd)
except Exception as exc:
raise ArgumentTypeError('Variable tfd must be a float or cast to float.')
if(not isinstance(x0, float)):
x0 = np.asarray(x0)
if(not isinstance(x0d, float)):
x0d = np.asarray(x0d)
if(not isinstance(pars, float)):
pars = np.asarray(pars)
if(not isinstance(parsd, float)):
parsd = np.asarray(parsd)
if(not isinstance(time_steps, float)):
time_steps = np.asarray(time_steps)
else:
msgErreur = 'The argument time_steps cannot be of size 1.';
raise InputArgumentError(msgErreur)
# -----------------------------------------------------------------------------------
#
# Tests on input variables
#
# -----------------------------------------------------------------------------------
# -----------------------------------------------------------------------------------
if options==None:
options = Options()
# -----------------------------------------------------------------------------------
#Variable x0 must be a vector
msgErreur = 'Variable x0 must be a one dimensional numpy array!';
if(isinstance(x0,np.ndarray)):
#Variable x0 must be a vector
if(x0.ndim==1):
n = x0.size
elif(x0.ndim==2):
(lig,col) = x0.shape
if(lig==1 or col==1):
n = x0.size
else :
raise ArgumentDimensionError(msgErreur)
else :
raise ArgumentDimensionError(msgErreur)
elif(isinstance(x0,float)):
n = 1
# -----------------------------------------------------------------------------------
# Variable x0d must be a vector or a matrix, with the same number of rows than x0
# and with k columns, with k to determine
msgErreur = 'Variable x0d must be a one dimensional numpy array!';
nx = n
if(isinstance(x0d,np.ndarray)):
#Variable x0d must be a vector
if(x0d.ndim==1):
ndx = x0d.size
elif(x0d.ndim==2):
(lig,col) = x0d.shape
if(lig==1 or col==1):
ndx = x0d.size
else :
raise ArgumentDimensionError(msgErreur)
else :
raise ArgumentDimensionError(msgErreur)
elif(isinstance(x0d,float)):
ndx = 1
#Test if ndx = nx
if(nx!=ndx):
raise ArgumentDimensionError('Variable x0d must have the same length than x0!')
# -----------------------------------------------------------------------------------
#Variable pars must be a vector or a float
msgErreur = 'Variable pars must be a one dimensional numpy array!';
if(isinstance(pars,np.ndarray)):
#Variable pars must be a vector
if(pars.ndim==1):
npars = pars.size
elif(pars.ndim==2):
(lig,col) = pars.shape
if(lig==1 or col==1):
npars = pars.size
else :
raise ArgumentDimensionError(msgErreur)
else :
raise ArgumentDimensionError(msgErreur)
elif(isinstance(pars,float)):
npars = 1
# -----------------------------------------------------------------------------------
#Variable parsd must be a vector or a float
msgErreur = 'Variable parsd must be a one dimensional numpy array!';
if(isinstance(parsd,np.ndarray)):
#Variable parsd must be a vector
if(parsd.ndim==1):
nparsd = parsd.size
elif(parsd.ndim==2):
(lig,col) = parsd.shape
if(lig==1 or col==1):
nparsd = parsd.size
else :
raise ArgumentDimensionError(msgErreur)
else :
raise ArgumentDimensionError(msgErreur)
elif(isinstance(parsd,float)):
nparsd = 1
#Test if npars = nparsd
if(npars!=nparsd):
raise ArgumentDimensionError('Variable parsd must have the same length ' + \
'than pars!')
# -----------------------------------------------------------------------------------
#Variable options must be a Options object
msgErreur = 'Variable options must be a Options object!';
if(not isinstance(options,Options)):
raise ArgumentTypeError(msgErreur)
# -----------------------------------------------------------------------------------
#
# prepare call to Fortran solver
#
# -----------------------------------------------------------------------------------
# -----------------------------------------------------------------------------------
[dw, iw, lsw, sw] = options.getAllOptions()
ndw = len(dw)
niw = len(iw)
nsw = len(lsw)
nswint = len(sw)
swint = [ord(c) for c in sw] #Code ascii de la chaine de caractere
# -----------------------------------------------------------------------------------
if(not isinstance(x0,float)):
x0.astype(dtype='d',order='F')
if(not isinstance(x0d,float)):
x0d.astype(dtype='d',order='F')
if(not isinstance(pars,float)):
pars.astype(dtype='d',order='F')
if(not isinstance(parsd,float)):
parsd.astype(dtype='d',order='F')
# -----------------------------------------------------------------------------------
if time_steps.size==0:
tspan = np.array([t0, tf])
elif time_steps.size==1:
msgErreur = 'The argument time_steps cannot be of size 1.';
raise InputArgumentError(msgErreur)
else:
if(t0==time_steps[0]) and (tf==time_steps[-1]):
tspan = time_steps
else:
msgErreur = 'You must have time_steps[0]=t0 and time_steps[-1]=tf!';
raise InputArgumentError(msgErreur)
nt = tspan.size
tspan.astype(dtype='d',order='F')
# -----------------------------------------------------------------------------------
# args
if not isinstance(args,tuple):
args = (args,)
# df
if pars_is_none:
mydf = lambda t, x, xd, pars, parsd: df(t, x, xd, *args, **kwargs)
else:
mydf = lambda t, x, xd, pars, parsd: df(t, x, xd, pars, parsd, *args, **kwargs)
# utile pour les derivees
def myf(t, x, pars):
if(not isinstance(x,float)):
xd = np.zeros(x.size)
else:
xd = 0.0
if(not isinstance(pars,float)):
parsd = np.zeros(pars.size)
else:
parsd = 0.0
dx, dxd = mydf(t, x, xd, pars, parsd)
return dx
# -----------------------------------------------------------------------------------
# -----------------------------------------------------------------------------------
# derivative wrt to t0
if t0d!=0:
x0d = x0d - myf(t0, x0, pars) * t0d
# -----------------------------------------------------------------------------------
[flag, nfev, nsteps, mesFor] = mod_ivpsolve.dexpfunpy(mydf,x0,x0d,pars,parsd,tspan, \
dw,iw,swint,lsw,n,npars,nt,ndw,niw,nsw,nswint)
# -----------------------------------------------------------------------------------
#
# post computations, prepare outputs
#
# -----------------------------------------------------------------------------------
# -----------------------------------------------------------------------------------
message = mesFor.strip().decode()
tout = mod_ivpsolve.tout
xout_for = mod_ivpsolve.xout
xdout_for = mod_ivpsolve.xdout
# tensorize
if(not isinstance(x0,float)):
x0shape = x0.shape
x0dshape = x0d.shape
if time_steps.size==2:
# tout
tout = np.array([tout[0], tout[-1]])
# xout
if(not isinstance(x0,float)):
x0_out = np.reshape(xout_for[:, 0], x0shape)
xf_out = np.reshape(xout_for[:, -1], x0shape)
else:
x0_out = xout_for[0, 0]
xf_out = xout_for[0, -1]
xout = np.concatenate(( x0_out[None], xf_out[None] ), axis=0)
# xdout
if(not isinstance(x0d,float)):
x0d_out = np.reshape(xdout_for[:, 0], x0dshape)
xfd_out = np.reshape(xdout_for[:, -1], x0dshape)
else:
x0d_out = xdout_for[0, 0]
xfd_out = xdout_for[0, -1]
# derivative wrt to tf
if tfd!=0:
xfd_out = xfd_out + myf(tout[-1], xf_out, pars) * tfd
# stack
xdout = np.concatenate(( x0d_out[None], xfd_out[None] ), axis=0)
else:
# xout
if(not isinstance(x0,float)):
xout = np.reshape(xout_for[:, 0], x0shape)[None]
else:
xout = xout_for[0, 0][None]
# xdout
if(not isinstance(x0d,float)):
xdout = np.reshape(xdout_for[:, 0], x0dshape)[None]
else:
xdout = xdout_for[0, 0][None]
#
for i in range(nsteps-1):
# xout
if(not isinstance(x0,float)):
x = np.reshape(xout_for[:, i+1], x0shape)
else:
x = xout_for[0, i+1]
# xdout
if(not isinstance(x0d,float)):
xd = np.reshape(xdout_for[:, i+1], x0dshape)
else:
xd = xdout_for[0, i+1]
# derivative wrt to tf
if tfd!=0:
xd = xd + myf(tout[i+1], x, pars) * tfd
xout = np.concatenate((xout, x[None]), axis=0)
xdout = np.concatenate((xdout, xd[None]), axis=0)
#
xf = xout[-1]
xfd = xdout[-1]
status = flag
if(status==1 or status==0):
success=True
else:
success=False
# -----------------------------------------------------------------------------------
sol = dict(xf=xf, xfd=xfd, tout=tout, xout=xout, xdout=xdout, status=status, \
success=success, message=message, nfev=nfev, nsteps=nsteps)
dexpsol = dExpSol(sol)
return dexpsol
# -----------------------------------------------------------------------------------
# -----------------------------------------------------------------------------------
class dExpSol(ExpSol):
def __init__(self, *args, **kwargs):
"""
Doc de __init__
"""
#args -- tuple of anonymous arguments
#kwargs -- dictionary of named arguments
ExpSol.__init__(self, *args, **kwargs)
self.xfd = self.sol['xfd']
self.xdout = self.sol['xdout']
# -----------------------------------------------------------------------------------
# -----------------------------------------------------------------------------------
[docs]def jexp(df, tf, t0, x0, dx0, pars=None, *, options=None, time_steps=None, \
args=(), kwargs={}):
r"""
The ``jexp`` function computes a numerical approximation of the solution of the
linearization of the initial value problem
.. math::
\frac{d}{d t} x = \dot{x} = f(t, x, pars), \quad x(t_0) = x_0,
where :math:`pars \in \mathbf{R}^k` is a vector of ``parameters``, :math:`t` is the
``time`` and :math:`x` is the ``state`` variable. The time :math:`t_0` is refered as
the ``initial time``, :math:`x_0` as the ``initial condition`` and :math:`f` as the
``dynamics``.
Hence, ``dexp`` approximates the solution of
.. math::
\dot{x} &= f(t, x, pars), \\
\dot{\delta x} &= \frac{\partial f}{\partial x}(t, x, pars) \cdot \delta x, \\
x(t_0) &= x_0, \\
\delta x(t_0) &= \delta x_0,
where :math:`\delta x_0` may be a vector of :math:`\mathbf{R}^n` (:math:`n` being
the dimension of the state :math:`x`) or a real matrix of :math:`n` rows and
:math:`k` columns.
The solution at the ``final time`` :math:`t_f` is denoted
:math:`(x(t_f, t_0, x_0, pars), \delta x(t_f, t_0, x_0, \delta x_0, pars))`.
Parameters
----------
df : callable
A vector function
``y, dy = df(t, x, dx)`` if ``pars`` is ``None``, where
``y = f(t, x)`` and ``dy = df/dx(t, x) . dx``
or
``y, dy = df(t, x, dx, pars)`` if ``pars`` is not
``None``, where ``y = f(t, x, pars)`` and
``dy = df/dx(t, x, pars) . dx``.
**Remark.** In the definition of ``df``, the argument ``dx`` is a vector
of the same shape as ``x``.
tf : :class:`float`
Final time
t0 : :class:`float`
Initial time`
x0 : float vector
Initial condition
dx0 : float vector or matrix
Initial condition :math:`\delta x_0` of the linear part of the equation
pars : float vector, optional, positional
Parameters
options : :meth:`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 ``df``
kwargs : dictionnary, key only
keywords optional arguments for ``df``
Returns
-------
sol
A dictionary/struct of outputs:
- ``sol.xf`` or ``sol.get('xf')`` : final point ``x(tf, t0, x0, pars)``;
- ``sol.dxf`` : principal linear part ``dx(tf, t0, dx0, x0, pars)``;
- ``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.dxout`` : principal linear part at integration time-steps:
``dxout[i] = dx(tout[i], t0, x0, dx0, pars)``;
- ``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.djexp
nutopy.ivp.Options
nutopy.errors
Notes
-----
- Use :class:`numpy.ndarray` to define float vectors.
- If ``dx0`` is not a matrix but a vector, then
``jexp(df, tf, t0, x0, dx0, pars) = dexp(df, tf, 0, t0, 0, x0, dx0, pars, 0)``.
Examples
--------
Example: dx/dx0
>>> import numpy as np
>>> from nutopy import ivp
>>> def df(t, x, dx):
... y = np.zeros(x.size)
... y[0] = -x[0] + x[1]
... y[1] = x[1]
... dy = np.zeros(x.size)
... dy[0] = -dx[0] + dx[1]
... dy[1] = dx[1]
... return (y, dy)
>>> tf = 1.0
>>> t0 = 0.0
>>> x0 = np.array([-1.5, 0.5])
>>> dx0 = np.eye(x0.size)
>>> sol = ivp.jexp(df, tf, t0, x0, dx0)
>>> print('dx/dx0(tf, t0, x0) = ', sol.dxf)
dx/dx0(tf, t0, x0) = [[0.36787944 1.17520119]
[0. 2.71828183]]
Raises
------
ArgumentTypeError
Variable tf must be a float or cast to float.
ArgumentTypeError
Variable t0 must be a float or cast to float.
ArgumentTypeError
Variable options must be a Options object.
ArgumentDimensionError
Variable x0 must be a one dimensional array.
ArgumentDimensionError
Variable dx0 must be either a one or two dimensional numpy array.
ArgumentDimensionError
Variable dx0 must have the same number of columns than x0.
ArgumentDimensionError
Variable pars must be a one dimensional array.
InputArgumentError
You must have time_steps[0]=t0 and time_steps[-1]=tf.
InputArgumentError
The argument time_steps cannot be of size 1.
"""
# -----------------------------------------------------------------------------------
#
# None
#
# -----------------------------------------------------------------------------------
pars_is_none=False
if pars is None:
pars_is_none=True # in this case, we should have f(t,x) and not f(t,x,pars)
pars = np.array([])
if time_steps is None:
time_steps = np.array([])
# -----------------------------------------------------------------------------------
#
# Cast of input variables
#
# -----------------------------------------------------------------------------------
if(not isinstance(t0, float)):
try:
t0 = float(t0)
except Exception as exc:
raise ArgumentTypeError('Variable t0 must be a float or cast to float.')
if(not isinstance(tf, float)):
try:
tf = float(tf)
except Exception as exc:
raise ArgumentTypeError('Variable tf must be a float or cast to float.')
if(not isinstance(x0, float)):
x0 = np.asarray(x0)
if(not isinstance(dx0, float)):
dx0 = np.asarray(dx0)
if(not isinstance(pars, float)):
pars = np.asarray(pars)
if(not isinstance(time_steps, float)):
time_steps = np.asarray(time_steps)
else:
msgErreur = 'The argument time_steps cannot be of size 1.';
raise InputArgumentError(msgErreur)
# -----------------------------------------------------------------------------------
#
# Tests on input variables
#
# -----------------------------------------------------------------------------------
# -----------------------------------------------------------------------------------
if options==None:
options = Options()
# -----------------------------------------------------------------------------------
#Variable x0 must be a vector
msgErreur = 'Variable x0 must be a one dimensional numpy array.';
if(isinstance(x0,np.ndarray)):
#Variable x0 must be a vector
if(x0.ndim==1):
n = x0.size
elif(x0.ndim==2):
(lig,col) = x0.shape
if(lig==1 or col==1):
n = x0.size
else :
raise ArgumentDimensionError(msgErreur)
else :
raise ArgumentDimensionError(msgErreur)
elif(isinstance(x0,float)):
n = 1
# else :
# raise ArgumentTypeError(msgErreur)
# -----------------------------------------------------------------------------------
# Variable dx0 must be a vector or a matrix, with the same number of rows than x0
# and with k columns, with k to determine
msgErreur = 'Variable dx0 must be either a one or two dimensional numpy array.';
if(isinstance(dx0,np.ndarray)):
#Test if dx0 is a vector or a matrix
if(dx0.ndim==1):
ndx = dx0.size
k = 1
elif(dx0.ndim==2):
(ndx,k) = dx0.shape
else :
raise ArgumentDimensionError(msgErreur)
elif(isinstance(dx0,float)):
ndx = 1
k = 1
#Test if ndx = nx
nx = n
if(nx!=ndx):
raise ArgumentDimensionError('Variable dx0 must have the same number of ' + \
'columns than x0.')
# -----------------------------------------------------------------------------------
#Variable pars must be a vector or a float
msgErreur = 'Variable pars must be a one dimensional numpy array.';
if(isinstance(pars,np.ndarray)):
#Variable pars must be a vector
if(pars.ndim==1):
npars = pars.size
elif(pars.ndim==2):
(lig,col) = pars.shape
if(lig==1 or col==1):
npars = pars.size
else :
raise ArgumentDimensionError(msgErreur)
else :
raise ArgumentDimensionError(msgErreur)
elif(isinstance(pars,float)):
npars = 1
# -----------------------------------------------------------------------------------
#Variable options must be a Options object
msgErreur = 'Variable options must be a Options object.';
if(not isinstance(options,Options)):
raise ArgumentTypeError(msgErreur)
# -----------------------------------------------------------------------------------
#
# prepare call to Fortran solver
#
# -----------------------------------------------------------------------------------
# -----------------------------------------------------------------------------------
[dw, iw, lsw, sw] = options.getAllOptions()
ndw = len(dw)
niw = len(iw)
nsw = len(lsw)
nswint = len(sw)
swint = [ord(c) for c in sw] #Code ascii de la chaine de caractere
# -----------------------------------------------------------------------------------
if(not isinstance(x0,float)):
x0.astype(dtype='d',order='F')
if(not isinstance(dx0,float)):
dx0.astype(dtype='d',order='F')
if(not isinstance(pars,float)):
pars.astype(dtype='d',order='F')
# -----------------------------------------------------------------------------------
if time_steps.size==0:
tspan = np.array([t0, tf])
elif time_steps.size==1:
msgErreur = 'The argument time_steps cannot be of size 1.';
raise InputArgumentError(msgErreur)
else:
if(t0==time_steps[0]) and (tf==time_steps[-1]):
tspan = time_steps
else:
msgErreur = 'You must have time_steps[0]=t0 and time_steps[-1]=tf!';
raise InputArgumentError(msgErreur)
nt = tspan.size
tspan.astype(dtype='d',order='F')
# -----------------------------------------------------------------------------------
# args
if not isinstance(args,tuple):
args = (args,)
#
# parsd is dummy here but we choose to call dexpfunpy
# we could have chosen to call expfunpy but in this case we have to adapt the rhs
#
if pars_is_none:
mydf = lambda t, x, dx, pars: df(t, x, dx, *args, **kwargs)
else:
if isinstance(args,tuple):
mydf = lambda t, x, dx, pars: df(t, x, dx, pars, *args, **kwargs)
else:
mydf = lambda t, x, dx, pars: df(t, x, dx, pars, args, **kwargs)
# -----------------------------------------------------------------------------------
[flag, nfev, nsteps, messageFortran] = mod_ivpsolve.jexpfunpy(mydf,x0,dx0,pars, \
tspan,dw,iw,swint,lsw,n,k,npars,nt,ndw,niw,nsw,nswint)
message = messageFortran.strip().decode()
# -----------------------------------------------------------------------------------
#
# post computations, prepare outputs
#
# -----------------------------------------------------------------------------------
# -----------------------------------------------------------------------------------
tout = mod_ivpsolve.tout
xout_for = mod_ivpsolve.xout
dxout_for = mod_ivpsolve.dxout
# tensorize
if(not isinstance(x0,float)):
x0shape = x0.shape
dx0shape = dx0.shape
if time_steps.size==2:
# tout
tout = np.array([tout[0], tout[-1]])
# xout
if(not isinstance(x0,float)):
x0_out = np.reshape(xout_for[:, 0], x0shape)
xf_out = np.reshape(xout_for[:, -1], x0shape)
else:
x0_out = xout_for[0, 0]
xf_out = xout_for[0, -1]
xout = np.concatenate(( x0_out[None], xf_out[None] ), axis=0)
# dxout
if(not isinstance(dx0,float)):
dx0_out = np.reshape(dxout_for[:, 0:k], dx0shape)
dxf_out = np.reshape(dxout_for[:, -k: ], dx0shape)
else:
dx0_out = dxout_for[0, 0]
dxf_out = dxout_for[0, -1]
dxout = np.concatenate(( dx0_out[None], dxf_out[None] ), axis=0)
else:
# xout
if(not isinstance(x0,float)):
xout = np.reshape(xout_for[:, 0], x0shape)[None]
else:
xout = xout_for[0, 0][None]
# dxout
icur = 0
if(not isinstance(dx0,float)):
dxout = np.reshape(dxout_for[:, icur:icur+k], dx0shape)[None]
else:
dxout = dxout_for[0, icur][None]
icur = icur + k
for i in range(nsteps-1):
# xout
if(not isinstance(x0,float)):
x = np.reshape(xout_for[:, i+1], x0shape)
else:
x = xout_for[0, i+1]
# dxout
if(not isinstance(dx0,float)):
dx = np.reshape(dxout_for[:, icur:icur+k], dx0shape)
else:
dx = dxout_for[0, icur]
icur = icur + k
#
xout = np.concatenate((xout, x[None]), axis=0)
dxout = np.concatenate((dxout, dx[None]), axis=0)
#
xf = xout[-1]
dxf = dxout[-1]
status = flag
if(status==1 or status==0):
success=True
else:
success=False
# -----------------------------------------------------------------------------------
sol = dict(xf=xf, dxf=dxf, tout=tout, xout=xout, dxout=dxout, status=status, \
success=success, message=message, nfev=nfev, nsteps=nsteps)
jexpsol = jexpSol(sol)
return jexpsol
# -----------------------------------------------------------------------------------
# -----------------------------------------------------------------------------------
class jexpSol(ExpSol):
def __init__(self, *args, **kwargs):
"""
Doc de __init__
"""
#args -- tuple of anonymous arguments
#kwargs -- dictionary of named arguments
ExpSol.__init__(self, *args, **kwargs)
self.dxf = self.sol['dxf']
self.dxout = self.sol['dxout']
# -----------------------------------------------------------------------------------
# -----------------------------------------------------------------------------------
[docs]def djexp(d2f, tf, tfd, t0, t0d, x0, x0d, dx0, dx0d, pars=None, parsd=None, *, \
options=None, time_steps=None, args=(), kwargs={}):
r"""
The ``djexp`` function computes a numerical approximation of the derivative of the
solution of the linearized initial value problem
.. math::
\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,
where :math:`pars \in \mathbf{R}^k` is a vector of ``parameters``, :math:`t` is the
``time``, :math:`x` is the ``state`` variable and :math:`y` its linearization.
The time :math:`t_0` is refered as the ``initial time``, :math:`x_0` as the
``initial condition``, :math:`y_0` as the `ìnitial condition`` of the linear part
and :math:`f` as the ``dynamics``.
The solution at the ``final time`` :math:`t_f` is then denoted by
:math:`x(t_f, t_0, x_0, pars)`, :math:`y(t_f, t_0, x_0, y_0, pars)`.
Hence, ``djexp`` approximates
.. math::
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
and
.. math::
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.
Let us introduce :math:`\delta x(t) = x'(t, t_0, x_0, pars) \cdot
(0, \delta t_0, \delta x_0, \delta pars)` and
:math:`\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, :math:`\delta x(t_f)` and :math:`\delta y(t_f)` are computed solving the
following `variational equations` (we omit :math:`(t, x, pars)`):
.. math::
\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).
The missing parts are given by:
.. math::
\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
.. math::
\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 :math:`x(t_f)` stands for :math:`x(t_f, t_0, x_0, pars)` and
:math:`y(t_f)` stands for :math:`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 : :class:`float`
Final time
tfd : :class:`float`
Increment :math:`\delta t_f`
t0 : :class:`float`
Initial time
t0d : :class:`float`
Increment :math:`\delta t_0`
x0 : float vector
Initial condition
x0d : float vector
Increment :math:`\delta x_0`
dx0 : float vector or matrix
Initial condition :math:`y_0` of the linear part of the equation
dx0d : float vector or matrix
Initial increment :math:`\delta y_0` of the linear part of the equation
pars : float vector, optional, positional
Parameters
parsd : float vector, optional, positional
Increment :math:`\delta pars`
options : :meth:`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
-----
- Use :class:`numpy.ndarray` to define float vectors.
- :meth:`nutopy.ivp.djexp` is the derivative of :meth:`nutopy.ivp.jexp`
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
------
ArgumentTypeError
Variable tf must be a float or cast to float.
ArgumentTypeError
Variable tfd must be a float or cast to float.
ArgumentTypeError
Variable t0 must be a float or cast to float.
ArgumentTypeError
Variable t0d must be a float or cast to float.
ArgumentTypeError
Variable options must be a Options object.
ArgumentDimensionError
Variable x0 must be a one dimensional array.
ArgumentDimensionError
Variable x0d must be a one dimensional array.
ArgumentDimensionError
Variable x0d must have the same length than x0.
ArgumentDimensionError
Variable dx0 must be either a one or two dimensional numpy array.
ArgumentDimensionError
Variable dx0 must have the same number of columns than x0.
ArgumentDimensionError
Variable dx0d must be either a one or two dimensional numpy array.
ArgumentDimensionError
Variable dx0d must have the same number of columns than x0.
ArgumentDimensionError
Variable dx0d must have the same number of rows than dx0.
ArgumentDimensionError
Variable pars must be a one dimensional array.
ArgumentDimensionError
Variable parsd must be a one dimensional array.
ArgumentDimensionError
Variable parsd must have the same length than pars.
InputArgumentError
Either give pars and parsd or any.
InputArgumentError
You must have time_steps[0]=t0 and time_steps[-1]=tf.
InputArgumentError
The argument time_steps cannot be of size 1.
"""
# -----------------------------------------------------------------------------------
#
# None
#
# -----------------------------------------------------------------------------------
pars_is_none=False
if pars is None:
pars_is_none=True # in this case, we should have f(t,x) and not f(t,x,pars)
pars = np.array([])
parsd_is_none=False
if parsd is None:
parsd_is_none=True
parsd = np.array([])
if pars_is_none != parsd_is_none:
raise InputArgumentError('Either give pars and parsd or any.')
if time_steps is None:
time_steps = np.array([])
# -----------------------------------------------------------------------------------
#
# Cast of input variables
#
# -----------------------------------------------------------------------------------
if(not isinstance(t0, float)):
try:
t0 = float(t0)
except Exception as exc:
raise ArgumentTypeError('Variable t0 must be a float or cast to float.')
if(not isinstance(tf, float)):
try:
tf = float(tf)
except Exception as exc:
raise ArgumentTypeError('Variable tf must be a float or cast to float.')
if(not isinstance(t0d, float)):
try:
t0d = float(t0d)
except Exception as exc:
raise ArgumentTypeError('Variable t0d must be a float or cast to float.')
if(not isinstance(tfd, float)):
try:
tfd = float(tfd)
except Exception as exc:
raise ArgumentTypeError('Variable tfd must be a float or cast to float.')
if(not isinstance(x0, float)):
x0 = np.asarray(x0)
if(not isinstance(x0d, float)):
x0d = np.asarray(x0d)
if(not isinstance(dx0, float)):
dx0 = np.asarray(dx0)
if(not isinstance(dx0d, float)):
dx0d = np.asarray(dx0d)
if(not isinstance(pars, float)):
pars = np.asarray(pars)
if(not isinstance(parsd, float)):
parsd = np.asarray(parsd)
if(not isinstance(time_steps, float)):
time_steps = np.asarray(time_steps)
else:
msgErreur = 'The argument time_steps cannot be of size 1.';
raise InputArgumentError(msgErreur)
# -----------------------------------------------------------------------------------
#
# Tests on input variables
#
# -----------------------------------------------------------------------------------
# -----------------------------------------------------------------------------------
if options==None:
options = Options()
# -----------------------------------------------------------------------------------
#Variable x0 must be a vector
msgErreur = 'Variable x0 must be a one dimensional numpy array.';
if(isinstance(x0,np.ndarray)):
#Variable x0 must be a vector
if(x0.ndim==1):
n = x0.size
elif(x0.ndim==2):
(lig,col) = x0.shape
if(lig==1 or col==1):
n = x0.size
else :
raise ArgumentDimensionError(msgErreur)
else :
raise ArgumentDimensionError(msgErreur)
elif(isinstance(x0,float)):
n = 1
# else :
# raise ArgumentTypeError(msgErreur)
# -----------------------------------------------------------------------------------
#Variable x0d must be a vector
msgErreur = 'Variable x0d must be a one dimensional numpy array.';
if(isinstance(x0d,np.ndarray)):
#Variable x0d must be a vector
if(x0d.ndim==1):
n = x0d.size
elif(x0d.ndim==2):
(lig,col) = x0d.shape
if(lig==1 or col==1):
n = x0d.size
else :
raise ArgumentDimensionError(msgErreur)
else :
raise ArgumentDimensionError(msgErreur)
elif(isinstance(x0d,float)):
n = 1
# else :
# raise ArgumentTypeError(msgErreur)
# -----------------------------------------------------------------------------------
# Variable dx0 must be a vector or a matrix, with the same number of rows than x0
# and with k columns, with k to determine
msgErreur = 'Variable dx0 must be either a one or two dimensional numpy array.';
if(isinstance(dx0,np.ndarray)):
#Test if dx0 is a vector or a matrix
if(dx0.ndim==1):
ndx = dx0.size
k = 1
elif(dx0.ndim==2):
(ndx,k) = dx0.shape
else :
raise ArgumentDimensionError(msgErreur)
elif(isinstance(dx0,float)):
ndx = 1
k = 1
# else :
# raise ArgumentTypeError(msgErreur)
# -----------------------------------------------------------------------------------
#Test if ndx = nx
nx = n
if(nx!=ndx):
raise ArgumentDimensionError('Variable dx0 must have the same number of ' + \
'columns than x0.')
# -----------------------------------------------------------------------------------
# Variable dx0d must be a vector or a matrix, with the same number of rows than x0
# and with k columns, with k to determine
msgErreur = 'Variable dx0d must be either a one or two dimensional numpy array.';
if(isinstance(dx0d,np.ndarray)):
#Test if dx0d is a vector or a matrix
if(dx0d.ndim==1):
ndx = dx0d.size
kbis= 1
elif(dx0d.ndim==2):
(ndx,kbis) = dx0d.shape
else :
raise ArgumentDimensionError(msgErreur)
elif(isinstance(dx0d,float)):
ndx = 1
kbis= 1
# else :
# raise ArgumentTypeError(msgErreur)
# -----------------------------------------------------------------------------------
#Test if ndx = nx
if(nx!=ndx):
raise ArgumentDimensionError('Variable dx0d must have the same number of ' + \
'columns than x0.')
if(k!=kbis):
raise ArgumentDimensionError('Variable dx0d must have the same number of ' + \
'rows than dx0.')
# -----------------------------------------------------------------------------------
#Variable pars must be a vector or a float
msgErreur = 'Variable pars must be a one dimensional numpy array.';
if(isinstance(pars,np.ndarray)):
#Variable pars must be a vector
if(pars.ndim==1):
npars = pars.size
elif(pars.ndim==2):
(lig,col) = pars.shape
if(lig==1 or col==1):
npars = pars.size
else :
raise ArgumentDimensionError(msgErreur)
else :
raise ArgumentDimensionError(msgErreur)
elif(isinstance(pars,float)):
npars = 1
# else :
# raise ArgumentTypeError(msgErreur)
# -----------------------------------------------------------------------------------
#Variable parsd must be a vector or a float
msgErreur = 'Variable parsd must be a one dimensional numpy array.';
if(isinstance(parsd,np.ndarray)):
#Variable parsd must be a vector
if(parsd.ndim==1):
nparsd = parsd.size
elif(parsd.ndim==2):
(lig,col) = parsd.shape
if(lig==1 or col==1):
nparsd = parsd.size
else :
raise ArgumentDimensionError(msgErreur)
else :
raise ArgumentDimensionError(msgErreur)
elif(isinstance(parsd,float)):
nparsd = 1
# else :
# raise ArgumentTypeError(msgErreur)
#Test if npars = nparsd
if(npars!=nparsd):
raise ArgumentDimensionError('Variable parsd must have the same length than ' + \
'pars!')
# -----------------------------------------------------------------------------------
#Variable options must be a Options object
msgErreur = 'Variable options must be a Options object.';
if(not isinstance(options,Options)):
raise ArgumentTypeError(msgErreur)
# -----------------------------------------------------------------------------------
#
# prepare call to Fortran solver
#
# -----------------------------------------------------------------------------------
# args
if not isinstance(args,tuple):
args = (args,)
# -----------------------------------------------------------------------------------
if pars_is_none:
myd2f = lambda t, x, dx1, dx2, pars, parsd: d2f(t, x, dx1, dx2, \
*args, **kwargs)
else:
myd2f = lambda t, x, dx1, dx2, pars, parsd: d2f(t, x, dx1, dx2, \
pars, parsd, *args, **kwargs)
# utile pour les derivees
def myf(t, x, pars):
if(not isinstance(x,float)):
dx1 = np.zeros(x.size)
dx2 = np.zeros(x.size)
else:
dx1 = 0.0
dx2 = 0.0
if(not isinstance(pars,float)):
parsd = np.zeros(pars.size)
else:
parsd = 0.0
y, dy, d2y = myd2f(t, x, dx1, dx2, pars, parsd)
return y
def mydf_pars0(t, x, dx, pars):
if(not isinstance(x,float)):
dx2 = np.zeros(x.size)
else:
dx2 = 0.0
if(not isinstance(pars,float)):
parsd = np.zeros(pars.size)
else:
parsd = 0.0
y, dy, d2y = myd2f(t, x, dx, dx2, pars, parsd)
return dy
def d2f_djexp(t, x, xd, dx, dxd, pars, parsd):
z, dz, d2z = myd2f(t, x, xd, dx, pars, parsd)
y = z
yd = dz
dyd = d2z
if(not isinstance(x,float)):
dx2 = np.zeros(x.size)
else:
dx2 = 0.0
if(not isinstance(pars,float)):
ld = np.zeros(pars.size)
else:
ld = 0.0
z, dz, d2z = myd2f(t, x, dx, dx2, pars, ld)
dy = dz
z, dz, d2z = myd2f(t, x, dxd, dx2, pars, ld)
dyd = dyd + dz
return y, yd, dy, dyd
# -----------------------------------------------------------------------------------
#
# x0
# x0d - t0d x f(t0, x0, pars))
# dx0
# dx0d - t0d x ( df/dx(t0, x0, pars) . dx0)
#
# Remark: dx0d and dx0 may be matrices
#
if t0d!=0:
x0d = x0d - t0d * myf(t0, x0, pars)
if((isinstance(dx0d,np.ndarray)) and (dx0d.ndim==2)):
for i in range(0,k):
dx0d[:, i] = dx0d[:, i] - t0d * mydf_pars0(t0, x0, dx0[:,i], pars)
else:
dx0d = dx0d - t0d * mydf_pars0(t0, x0, dx0, pars)
# -----------------------------------------------------------------------------------
[dw, iw, lsw, sw] = options.getAllOptions()
ndw = len(dw)
niw = len(iw)
nsw = len(lsw)
nswint = len(sw)
swint = [ord(c) for c in sw] #Code ascii de la chaine de caractere
# -----------------------------------------------------------------------------------
if(not isinstance(x0,float)):
x0.astype(dtype='d',order='F')
if(not isinstance(x0d,float)):
x0d.astype(dtype='d',order='F')
if(not isinstance(dx0,float)):
dx0.astype(dtype='d',order='F')
if(not isinstance(dx0d,float)):
dx0d.astype(dtype='d',order='F')
if(not isinstance(pars,float)):
pars.astype(dtype='d',order='F')
if(not isinstance(parsd,float)):
parsd.astype(dtype='d',order='F')
# -----------------------------------------------------------------------------------
if time_steps.size==0:
tspan = np.array([t0, tf])
elif time_steps.size==1:
msgErreur = 'The argument time_steps cannot be of size 1.';
raise InputArgumentError(msgErreur)
else:
if(t0==time_steps[0]) and (tf==time_steps[-1]):
tspan = time_steps
else:
msgErreur = 'You must have time_steps[0]=t0 and time_steps[-1]=tf!';
raise InputArgumentError(msgErreur)
nt = tspan.size
tspan.astype(dtype='d',order='F')
# -----------------------------------------------------------------------------------
[flag, nfev, nsteps, mesFor] = mod_ivpsolve.djexpfunpy(d2f_djexp,x0,x0d,dx0,dx0d, \
pars,parsd,tspan,dw,iw,swint,lsw,n,k,npars,nt,ndw,niw,nsw,nswint)
message = mesFor.strip().decode()
# -----------------------------------------------------------------------------------
tout = mod_ivpsolve.tout
xout_for = mod_ivpsolve.xout
xdout_for = mod_ivpsolve.xdout
dxout_for = mod_ivpsolve.dxout
dxdout_for = mod_ivpsolve.dxdout
# tensorize
if(not isinstance(x0,float)):
x0shape = x0.shape
x0dshape = x0d.shape
dx0shape = dx0.shape
dx0dshape = dx0d.shape
if time_steps.size==2:
# tout
tout = np.array([tout[0], tout[-1]])
# xout
if(not isinstance(x0,float)):
x0_out = np.reshape(xout_for[:, 0], x0shape)
xf_out = np.reshape(xout_for[:, -1], x0shape)
else:
x0_out = xout_for[0, 0]
xf_out = xout_for[0, -1]
xout = np.concatenate(( x0_out[None], xf_out[None] ), axis=0)
# xdout
if(not isinstance(x0d,float)):
x0d_out = np.reshape(xdout_for[:, 0], x0dshape)
xfd_out = np.reshape(xdout_for[:, -1], x0dshape)
else:
x0d_out = xdout_for[0, 0]
xfd_out = xdout_for[0, -1]
# derivative wrt to tf
if tfd!=0:
xfd_out = xfd_out + myf(tout[-1], xf_out, pars) * tfd
xdout = np.concatenate(( x0d_out[None], xfd_out[None] ), axis=0)
# dxout
if(not isinstance(dx0,float)):
dx0_out = np.reshape(dxout_for[:, 0:k], dx0shape)
dxf_out = np.reshape(dxout_for[:, -k: ], dx0shape)
else:
dx0_out = dxout_for[0, 0]
dxf_out = dxout_for[0, -1]
dxout = np.concatenate(( dx0_out[None], dxf_out[None] ), axis=0)
# dxdout
if(not isinstance(dx0d,float)):
dx0d_out = np.reshape(dxdout_for[:, 0:k], dx0dshape)
dxfd_out = np.reshape(dxdout_for[:, -k: ], dx0dshape)
else:
dx0d_out = dxdout_for[0, 0]
dxfd_out = dxdout_for[0, -1]
# derivative wrt to tf
if tfd!=0:
if((isinstance(dx0d,np.ndarray)) and (dx0d.ndim==2)):
for i in range(0,k):
dxfd_out[:,i] = dxfd_out[:,i] + tfd * mydf_pars0(tout[-1], xf_out, \
dxf_out[:,i], pars)
else:
dxfd_out = dxfd_out + tfd * mydf_pars0(tout[-1], xf_out, dxf_out, pars)
dxdout = np.concatenate(( dx0d_out[None], dxfd_out[None] ), axis=0)
else:
# xout
if(not isinstance(x0,float)):
xout = np.reshape(xout_for[:, 0], x0shape)[None]
else:
xout = xout_for[0, 0][None]
# xdout
if(not isinstance(x0d,float)):
xdout = np.reshape(xdout_for[:, 0], x0dshape)[None]
else:
xdout = xdout_for[0, 0][None]
icur = 0
# dxout
if(not isinstance(dx0,float)):
dxout = np.reshape(dxout_for[:, icur:icur+k], dx0shape)[None]
else:
dxout = dxout_for[0, icur][None]
# dxdout
if(not isinstance(dx0d,float)):
dxdout = np.reshape(dxdout_for[:, icur:icur+k], dx0dshape)[None]
else:
dxdout = dxdout_for[0, icur][None]
icur = icur + k
for i in range(nsteps-1):
# xout
if(not isinstance(x0,float)):
x = np.reshape(xout_for[:, i+1], x0shape)
else:
x = xout_for[0, i+1]
# xdout
if(not isinstance(x0d,float)):
xd = np.reshape(xdout_for[:, i+1], x0dshape)
else:
xd = xdout_for[0, i+1]
# derivative wrt to tf: at time t (for t0, we already have it,
# that's why we start at i=1)
if tfd!=0:
xd = xd + myf(tout[i+1], x, pars) * tfd
# dxout
if(not isinstance(dx0,float)):
dx = np.reshape(dxout_for[:, icur:icur+k], dx0shape)
else:
dx = dxout_for[0, icur]
# dxdout
if(not isinstance(dx0d,float)):
dxd = np.reshape(dxdout_for[:, icur:icur+k], dx0dshape)
else:
dxd = dxdout_for[0, icur]
# derivative wrt to tf
if tfd!=0:
if((isinstance(dx0d,np.ndarray)) and (dx0d.ndim==2)):
for i in range(0,k):
dxd[:,i] = dxd[:,i] + tfd * mydf_pars0(tout[i+1], x, dx[:,i], pars)
else:
dxd = dxd + tfd * mydf_pars0(tout[i+1], x, dx, pars)
#
icur = icur + k
#
xout = np.concatenate((xout, x[None]), axis=0)
xdout = np.concatenate((xdout, xd[None]), axis=0)
dxout = np.concatenate((dxout, dx[None]), axis=0)
dxdout = np.concatenate((dxdout, dxd[None]), axis=0)
#
xf = xout[-1]
xfd = xdout[-1]
dxf = dxout[-1]
dxfd = dxdout[-1]
# -----------------------------------------------------------------------------------
status = flag
if(status==1 or status==0):
success=True
else:
success=False
# -----------------------------------------------------------------------------------
sol = dict( xf=xf, xfd=xfd, dxf=dxf, dxfd=dxfd, tout=tout, \
dxout=dxout, dxdout=dxdout, xout=xout, xdout=xdout, \
status=status, success=success, message=message, \
nfev=nfev, nsteps=nsteps)
djexpsol = djexpSol(sol)
return djexpsol
# -----------------------------------------------------------------------------------
# -----------------------------------------------------------------------------------
class djexpSol(ExpSol):
def __init__(self, *args, **kwargs):
"""
Doc de __init__
"""
#args -- tuple of anonymous arguments
#kwargs -- dictionary of named arguments
ExpSol.__init__(self, *args, **kwargs)
self.xfd = self.sol['xfd']
self.xdout = self.sol['xdout']
self.dxf = self.sol['dxf']
self.dxout = self.sol['dxout']
self.dxfd = self.sol['dxfd']
self.dxdout = self.sol['dxdout']
# -----------------------------------------------------------------------------------
# -----------------------------------------------------------------------------------
[docs]class Options(options.Options):
"""
Options for :meth:`nutopy.ivp.exp`, :meth:`nutopy.ivp.dexp`,
:meth:`nutopy.ivp.jexp` and :meth:`nutopy.ivp.djexp`.
Attributes
----------
MaxSteps : :class:`int`, default is ``100000``
Maximum number of integration steps.
MaxStepSize : :class:`float`, default is ``0.0``
Maximum step size during integration.
If it is equal to ``0.0`` then the integrator ``SolverMethod``
has no constraint on the maximal step size.
SolverMethod : :class:`str`, default is ``'dopri5'``
Integrator name.
See Refs [1]_ and [2]_ for a description of the methods.
The integrators with step size control are interfaces of Fortran
softwares available
`here <http://www.unige.ch/~hairer/software.html>`_.
- With step size control:
- Explicit:
dopri5, dop853.
- Implicit:
radau5, radau9, radau13, radau (adaptative order 5, 9 and 13).
- No step size control:
- Explicit:
euler_explicite, runge, heun, rk4, rk5.
- Implicit:
euler_implicite or euler, midpoint,
gauss4, gauss6, gauss8,
radauia1, radauia3, radauia5,
radauiia1, radauiia3, radauiia5,
radaus (symplectic radau of order 5),
lobatto4, lobatto6,
lobattoiiia2, lobattoiiia4, lobattoiiia6,
lobattoiiib2, lobattoiiib4, lobattoiiib6,
lobattoiiic2, lobattoiiic4, lobattoiiic6,
sdirk3, sdirk4l, sdirk4a, dirk5.
TolAbs : :class:`float`, default is ``1e-10``
Absolute error tolerance.
TolRel : :class:`float`, default is ``1e-8``
Relative error tolerance.
Examples
--------
>>> from nutopy import ivp
Constructor usage
>>> options = ivp.Options()
>>> options = ivp.Options(SolverMethod='dopri5', TolAbs=1e-8)
>>> options = ivp.Options({'SolverMethod' : 'dopri5', 'TolAbs' : 1e-8})
Update
>>> options.update(SolverMethod='dopri5')
Get
>>> solver = options.get('SolverMethod')
References
----------
.. [1] E. Hairer, S. P. Norsett & G. Wanner, Solving Ordinary Differential Equations I, Nonstiff Problems,
vol 8 of Springer Serie in Computational Mathematics, Springer-Verlag, second edn (1993).
.. [2] E. Hairer & G. Wanner, Solving Ordinary Differential Equations II, Stiff and Differential-Algebraic Problems,
vol 14 of Springer Serie in Computational Mathematics, Springer-Verlag, second edn (1996).
"""
#On peut initialiser les options de 4 manieres :
# soit avec dictionnaire
# soit avec une liste de dictionnaire
# soit avec une liste de param=value
# soit avec une liste de dictionnaire et de param=value
#
# Exemple
#
# def f(*args, **kwargs):
# print 'args: ', args, ' kwargs: ', kwargs
#
# >>> f('a')
# args: ('a',) kwargs: {}
# >>> f(ar='a')
# args: () kwargs: {'ar': 'a'}
# >>> f(1,2,param=3)
# args: (1, 2) kwargs: {'param': 3}
#
def __init__(self, *args, **kwargs):
"""
Doc de __init__
Parameters
----------
args : tuple
Anonymous arguments
kwargs : dictionary
Named arguments
"""
self.name = 'ivp'
#On donne les valeurs par defaut de options
self.options = dict(MaxSteps = 100000,
MaxStepSize = 0.0,
SolverMethod = 'dopri5',
TolAbs = 1e-10,
TolRel = 1e-8)
#args -- tuple of anonymous arguments
#kwargs -- dictionary of named arguments
self.update(*args, **kwargs)
def getAllOptions(self):
"""
Doc de getAllOptions
"""
# Cette methode permet de regrouper les options a donner en parametre aux
# fonctions fortran faisant l interface
dw = [ self.options['MaxStepSize'] ,
self.options['TolAbs'] ,
self.options['TolRel'] ]
iw = [ self.options['MaxSteps'] ]
lsw = [ len(self.options['SolverMethod']) ]
sw = self.options['SolverMethod']
return [dw, iw, lsw, sw]