Problem Class

Note

Here we describe the standard problem classes. From these users might derive problem classes for their particular needs. This is mainly the case when the right hand side function (rhs / res) has to be provided with additional information passed as instance attributes. For an example see cvode_with_disc.py.

class assimulo.problem.cProblem
finalize()

Method for specifying the finalization options when the simulation have finished.

handle_event()

Defines how to handle a discontinuity. This functions gets called when a discontinuity has been found in the supplied event functions. The solver is the solver attribute while the event_info is a list of length 2 where the first element is a list containing information about state events and the second element is a Boolean for indicating if there has been a time event. If there has not been a state event the first element is an empty list. The state event list contains a set of integers of values (-1,0,1), the values indicates which state event has triggered (determined from state_event(…) ) and the value indicates to where the state event is ‘headed’.

initialize()

Method for specializing initiation.

reset()

Resets a problem to its default values.

Explicit Problem

class assimulo.problem.Explicit_Problem

Problem for our explicit integrators (ODEs). A problem consists of the right-hand-side and some initial conditions.

Parameters:

rhs 
    Function that calculates the right-hand-side. Depending on
    the problem and the support of the solver, this function has
    the following input parameters:
    
        rhs(t,y)      - Normal ODE
        rhs(t,y,sw)   - An ODE with different modes, sw is a list of
                        switches (Boolean list) which should be held
                        constant during the integration and only be
                        changed when an event have occured. Used together
                        with event functions.
        rhs(t,y,p)  - An ODE with parameters for which sensitivities
                        should be calculated.
        rhs(t,y,sw,p) - An ODE with both parameters and switches.
        
        Returns:
            A numpy array of size len(y).

y0
    Defines the starting values 
t0
    Defines the starting time
p0 (Depending on if the solver supports sensitivity calculations)
    Parameters for which sensitivites are to be calculated
sw0 (Depending on if the solver supports state events)
    Defines the starting values of the switches. 
    Should be a list of Booleans.

Signature of default or user provided methods. Their use is solver dependent.

def state_events(self ,t ,y, sw)
    Defines the event (root) functions.
    
    Returns:
        A numpy array.
        
def time_events(self, t, y, sw)
    Defines the time events. This function should return
    the next time-point for a time event. At a time-event
    the usual method handle_event is called for the specific
    handling. If there are no more time events. This function
    should return None.
    
    Returns:
        Float
            The time-point for the next time-event.
        None
            No time-event.
    
def jac(self, t, y, sw=None)
    Defines the jacobian. J=df/dx.
    
    Returns:
        A numpy matrix of size len(y)*len(y).
        
def jacv(self, t, y, fy, v)
    Defines a Jacobian Vector product. df/dx*v.
    
    Returns:
        A numpy vector of size len(y).

def handle_result(self, solver, t, y)
    Method for specifying how the result is handled. 
    By default the data is stored in two vectors, solver.(t_sol/y_sol). If
    the problem to be solved also involve sensitivities these results are
    stored in p_sol
    
def handle_event(self, object solver, event_info):
    Defines how to handle a discontinuity. This functions is called when
    a discontinuity has been found in the supplied event functions. The solver
    is the solver attribute. The event_info is a list of length 2 where
    the first element is a list containing information about state events and
    the second element is a Boolean  indicating if there has been a time
    event. If there has not been a state event the first element is an empty
    list. The state event list contains a set of integers of values (-1,0,1),
    the values indicate which state event has triggered (determined from 
    state_event(...) ) and the value indicates to where the state event is 'headed'.

Implicit Problem

class assimulo.problem.Implicit_Problem

Problem for our implicit integrators (DAEs). A problem consists of the residual function and some initial conditions.

Parameters

res   
    Function that calculates the residual. Depending on
    the problem and the support of the solver, this function can
    have the following input parameters.
    
        res(t,y,yd)   - Normal DAE
        res(t,y,yd,sw)   - An DAE with different modes, sw is a list of
                           switches (boolean list) which should be held
                           constant during the integration and only be
                           changed when an event have occured. Used together
                           with event functions.
        res(t,y,yd,p)   - An DAE with parameters for which sensitivities
                           should be calculated.
        res(t,y,yd,sw,p) - An DAE with both parameters and switches.
        
        Returns:
            A numpy array of size len(y).
y0
    Defines the starting values of y0.
yd0
    Defines the starting values of yd0.
t0
    Defines the starting time.
p0 (Depending on if the solver supports sensitivity calculations)
    Parameters for which sensitivites are to be calculated
sw0 (Depending on if the solver supports state events)
    Defines the starting values of the switches. 
    Should be a list of Booleans.

Parameters (optionally contained in class)

algvar
    Defines the differential and algebraic components of the problem.
    Should be a list of integers. For more information, see the
    property algvar in IDA.

Signature of default or user provided methods. Their use is solver dependent.

def state_events(self ,t ,y ,yd, sw)
    Defines the event (root) functions.
    
    Returns:
        A numpy array.
    
def time_events(self, t, y, yd, sw)
    Defines the time events. This function should return
    the next time-point for a time event. At a time-event
    the usual method handle_event is called for the specific
    handling. If there are no more time events. This function
    should return None.
    
    Returns:
        Float
            The time-point for the next time-event.
        None
            No time-event.
    
def jac(self, c, t, y, yd, sw)
    Defines the Jacobian, which should be of the form
    J = dF/dx + c*dF/dx'.
    
    Returns:
        A numpy array of size len(y)*len(y).
        
def handle_result(self, solver, t, y, yd)
    Method for specifying how the result is  handled. 
    By default the data is stored in three vectors, solver.(t_sol/y_sol/yd_sol). 
    If the problem to be solved also involve sensitivities these results are
    stored in p_sol
    
def handle_event(self, object solver, event_info):
    Defines how to handle a discontinuity. This functions gets called when
    a discontinuity has been found in the supplied event functions. The solver
    is the solver attribute. The event_info is a list of length 2 where
    the first element is a list containing information about state events and
    the second element is a Boolean  indicating if there has been a time
    event. If there has not been a state event the first element is an empty
    list. The state event list contains a set of integers of values (-1,0,1),
    the values indicate which state event has triggered (determined from 
    state_event(...) ) and the value indicates to where the state event is 'headed'.

Mechanical Problem

class assimulo.special_systems.cMechanical_System

Special problem class for (constrained) mechanical systems:

\begin{eqnarray*} \dot{p} & = & v \\ M(p) \dot{v} & = & f(t,p,v)-G(p)^\mathrm{T} \lambda \\ 0 & = & g(p) \end{eqnarray*}

Parameters

n_p   number of position variables
forces(t, p, v) applied forces with len(y)=2*np
n_la  number of constraints
GT(p)  n_p x n_la constraint matrix, p=y[:np]
      (required if n_la > 0)
pos0, vel0, lam0
    Defines the initial values for positions, velocities 
    and constraint forces
posd0, veld0
    Defines the initial derivatives for positions and velocities 
t0
    Defines the initial time
mass_matrix(p)  n_p x n_p nonsingular mass matrix
      (if not defined it is assumed to be the identity matrix)
constr_3(t,y)  n_la index-3 constraints
constr_2(t,y)  n_la index-2 constraints
      (optional)
constr_1(t,y) n_la index-1 constraints
index  sets the type of equations to be solved
     'ind3' Index-3 DAE (position constraints)
     'ind2' Index-2 DAE (velocity constraints)
     'ind1' Index-1 DAE (lambda constraints)
     'ostab2' overdetermined stabilized index 2
     'ostab1' overdetermined stabilized index 1