siconos.kernel.EulerMoreauOSI (Python class)

class siconos.kernel.EulerMoreauOSI(*args)[source]

Bases: siconos.kernel.OneStepIntegrator

One Step time Integrator for First Order Dynamical Systems.

This integrator is the work horse of the event–capturing time stepping schemes for first order systems. It is mainly based on some extensions of the Backward Euler and \(\theta-\gamma\) schemes proposed in the pionnering work of J.J. Moreau for the sweeping process

J.J. Moreau. Evolution problem associated with a moving convex set in a Hilbert space. Journal of Differential Equations, 26, pp 347–374, 1977.

Variants are now used to integrate LCS, Relay systems, Higher order sweeping process see for instance

Consistency of a time-stepping method for a class of piecewise linear networks

M.K. Camlibel, W.P.M.H. Heemels, and J.M. Schumacher IEEE Transactions on Circuits and Systems I, 2002, 49(3):349–357

Numerical methods for nonsmooth dynamical systems: applications in mechanics and electronics

V Acary, B Brogliato Springer Verlag 2008

Convergence of time-stepping schemes for passive and extended linear complementarity systems L. Han, A. Tiwari, M.K. Camlibel, and J.-S. Pang SIAM Journal on Numerical Analysis 2009, 47(5):3768-3796

On preserving dissipativity properties of linear complementarity dynamical systems with the &theta-method

Greenhalgh Scott, Acary Vincent, Brogliato Bernard Numer. Math., , 2013.

Main time–integration schemes are based on the following \(\theta-\gamma\) scheme

\begin{cases} \label{eq:toto1} M x_{k+1} = M x_{k} +h\theta f(x_{k+1},t_{k+1})+h(1-\theta) f(x_k,t_k) + h \gamma r(t_{k+1}) + h(1-\gamma)r(t_k) \\[2mm] y_{k+1} = h(t_{k+1},x_{k+1},\lambda _{k+1}) \\[2mm] r_{k+1} = g(x_{k+1},\lambda_{k+1},t_{k+1})\\[2mm] \mbox{nslaw} ( y_{k+1} , \lambda_{k+1}) \end{cases}

where \(\theta = [0,1]\) and \(\gamma \in [0,1]\). As in Acary & Brogliato 2008, we call the previous problem the ``one–step nonsmooth problem’‘.

Another variant can also be used (FullThetaGamma scheme)

\begin{cases} M x_{k+1} = M x_{k} +h f(x_{k+\theta},t_{k+1}) + h r(t_{k+\gamma}) \\[2mm] y_{k+\gamma} = h(t_{k+\gamma},x_{k+\gamma},\lambda _{k+\gamma}) \\[2mm] r_{k+\gamma} = g(x_{k+\gamma},\lambda_{k+\gamma},t_{k+\gamma})\\[2mm] \mbox{nslaw} ( y_{k+\gamma} , \lambda_{k+\gamma}) \end{cases}

EulerMoreauOSI class is used to define some time-integrators methods for a list of first order dynamical systems. A EulerMoreauOSI instance is defined by the value of theta and possibly gamma and the list of concerned dynamical systems.

Each DynamicalSystem is associated to a SiconosMatrix, named “W”, which is the “iteration” matrix. W matrices are initialized and computed in initializeIterationMatrixW and computeW. Depending on the DS type, they may depend on time t and DS state x.

For first order systems, the implementation uses _r for storing the the input due to the nonsmooth law. This EulerMoreauOSI scheme assumes that the relative degree is zero or one and one level for _r is sufficient

Main functions:

  • computeFreeState(): computes xfree (or vfree), dynamical systems state without taking non-smooth part into account
  • updateState(): computes x (q,v), the complete dynamical systems states.

See User’s guide, for details.

Generated class (swig), based on C++ header Program listing for file kernel/src/simulationTools/EulerMoreauOSI.hpp.

Constructors

EulerMoreauOSI(double theta)

constructor from theta value only

Parameters:theta – value for all DS.
EulerMoreauOSI(double theta, double gamma)

constructor from theta value only

Parameters:
  • theta – value for all linked DS.
  • gamma – value for all linked DS.
W(DynamicalSystem ds) -> array_like (np.float64, 2D)[source]

get W corresponding to DynamicalSystem ds

Parameters:ds – a pointer to DynamicalSystem
Returns:pointer to a SiconosMatrix
WBoundaryConditions(DynamicalSystem ds) -> array_like (np.float64, 2D)[source]

get WBoundaryConditions corresponding to DynamicalSystem ds

Parameters:ds – a pointer to DynamicalSystem, optional, default = NULL. get WBoundaryConditions[0] in that case
Returns:pointer to a SiconosMatrix
computeFreeOutput(InteractionsGraph::VDescriptor vertex_inter, OneStepNSProblem *osnsp) → None[source]

integrates the Interaction linked to this integrator, without taking non-smooth effects into account

Parameters:
  • vertex_inter – of the interaction graph
  • osnsp – pointer to OneStepNSProblem
computeFreeState() → None[source]

Perform the integration of the dynamical systems linked to this integrator without taking into account the nonsmooth input r.

computeKhat(Interaction inter, array_like (np.float64, 2D) m, VectorOfSMatrices workM, double h) → None[source]
computeResidu() → double[source]

Computes the residuFree and residu of all the DynamicalSystems.

Returns:the maximum of the 2-norm over all the residu
computeResiduInput(double time, InteractionsGraph indexSet) → double[source]

compute the residu of the input of the relation (R or p) This computation depends on the type of OSI

Parameters:
  • time – time of computation
  • indexSet – the index set of the interaction that are concerned
computeResiduOutput(double time, InteractionsGraph indexSet) → double[source]

compute the residu of the output of the relation (y) This computation depends on the type of OSI

Parameters:
  • time – time of computation
  • indexSet – the index set of the interaction that are concerned
computeW(double time, DynamicalSystem ds, DynamicalSystemsGraph::VDescriptor dsv, array_like (np.float64, 2D) W) → None[source]

compute W EulerMoreauOSI matrix at time t

Parameters:
  • time – the current time
  • ds – the DynamicalSystem
  • dsv – a descriptor of the ds on the graph (redundant to avoid invocation)
  • W – the matrix to compute
computeWBoundaryConditions(DynamicalSystem ds) → None[source]

compute WBoundaryConditionsMap[ds] EulerMoreauOSI matrix at time t

Parameters:ds – a pointer to DynamicalSystem
display() → None[source]

Displays the data of the EulerMoreauOSI’s integrator.

gamma() → double[source]

get gamma

Returns:a double
getW(DynamicalSystem ds=DynamicalSystem()) -> array_like (np.float64, 2D)[source]

get the value of W corresponding to DynamicalSystem ds

Parameters:ds – a pointer to DynamicalSystem, optional, default = NULL. get W[0] in that case
Returns:SimpleMatrix
getWBoundaryConditions(DynamicalSystem ds=DynamicalSystem()) -> array_like (np.float64, 2D)[source]

get the value of WBoundaryConditions corresponding to DynamicalSystem ds

Parameters:ds – a pointer to DynamicalSystem, optional, default = NULL. get WBoundaryConditions[0] in that case
Returns:SimpleMatrix
initializeIterationMatrixW(double time, DynamicalSystem ds) → None[source]

initialize iteration matrix W EulerMoreauOSI matrix at time t

Parameters:
  • time – the time (double)
  • ds – a pointer to DynamicalSystem
initializeIterationMatrixWBoundaryConditions(DynamicalSystem ds) → None[source]

initialize iteration matrix WBoundaryConditionsMap[ds] EulerMoreauOSI

Parameters:ds – a pointer to DynamicalSystem
initializeWorkVectorsForDS(double t, DynamicalSystem ds) → None[source]

initialization of the EulerMoreauOSI integrator; for linear time invariant systems, we compute time invariant operator (example : W)

initialization of the work vectors and matrices (properties) related to one dynamical system on the graph and needed by the osi

Parameters:
  • t – time of initialization
  • ds – the dynamical system
initializeWorkVectorsForInteraction(Interaction inter, InteractionProperties interProp, DynamicalSystemsGraph DSG) → None[source]

initialization of the work vectors and matrices (properties) related to one interaction on the graph and needed by the osi

Parameters:
  • inter – the interaction
  • interProp – the properties on the graph
  • DSG – the dynamical systems graph
integrate(double tinit, double tend, double tout, int useless) → None[source]

integrate the system, between tinit and tend (->iout=true), with possible stop at tout (->iout=false)

Parameters:
  • tinit – initial time
  • tend – end time
  • tout – real end time
  • useless – flag (for EulerMoreauOSI, used in LsodarOSI)
numberOfIndexSets() → int[source]

get the number of index sets required for the simulation

Returns:unsigned int
prepareNewtonIteration(double time) → None[source]

computes all the W matrices

Parameters:time – current time
setGamma(double newGamma) → None[source]

set the value of gamma

Parameters:newGamma – a double
setTheta(double newTheta) → None[source]

set the value of theta

Parameters:newTheta – a double
setUseGamma(bool b) → None[source]

set the boolean to indicate that we use gamma

Parameters:b – true if gamma has to be used, false otherwise
setUseGammaForRelation(bool newUseGammaForRelation) → None[source]

set the boolean to indicate that we use gamma for the relation

Parameters:newUseGammaForRelation – a bool
theta() → double[source]

get theta

Returns:a double
updateInput(*args)[source]

Warning - Overloaded function : multiple signatures available, check prototypes below.

updateInput(double time) → None[source]

update the input of the Interaction attached to this Integrator

updateInput(double time, int level) → None[source]

update the input of the Interaction attached to this Integrator

Parameters:
  • time – current time
  • level – level of interest for the dynamics
updateOutput(*args)[source]

Warning - Overloaded function : multiple signatures available, check prototypes below.

updateOutput(double time) → None[source]

update the output of the Interaction attached to this Integrator

updateOutput(double time, int level) → None[source]

update the output of the Interaction attached to this Integrator

Parameters:
  • time – current time
  • level – level of interest for the dynamics
updateState(int level) → None[source]

updates the state of the Dynamical Systems

Parameters:level – the level of interest for the dynamics: not used at the time
useGamma() → bool[source]

get bool useGamma

Returns:a bool
useGammaForRelation() → bool[source]

get bool gammaForRelation for the relation

Returns:a