Ordinary Differential Equation (ODE) Object

In order to use the numerical integrator, you must supply it with an ODE object derived from the abstract class CDifferentialEquation. This class defines a number of routines that allow the ODE to communicate with the integrator to perform such functions as returning the state derivative, supplying output, handling events, providing default initial conditions, or resetting for a new integration run.

This section describes the various variables and member functions which you must override in order to provide the appropriate behavior for the ODE object.

Creating a New ODE Object

You must create a new object to represent the differential equation. This object must be derived from the abstract class CDifferentialEquation, as a public base class. You must define a number of member functions and initialize the member variables based on the capabilities you require your ODE to handle. This section describes some of the considerations in creating your ODE object. The following sections describe the member functions and variables in more detail.

At a minimum you must provide a constructor and the derivative function. The constructor is responsible for setting the number of equations and sizing the derivative vector. The derivative function is where you place the ODE you are integrating. It must be in first-order form, and return the derivative given the current time and state vector.

There are a number of other functions that you can override as needed. These functions allow you to handle events, provide default values to the integrator, respond to reset and update requests from the integrator, and output values at the current output time.

Function Overview

The base ODE class has a large number of virtual member functions defined for various purposes. There are functions to handle events, provide default values to the integrator, respond to reset and update requests from the integrator, and output values at the current output time.

Table 20 gives a brief description of the general nature of each group of functions. The actual functions are defined in the sections to follow. For each, the usage is defined, along with any required or optional inputs, and a description of the purpose of the routine and any special considerations in using it. The final sections list the member functions that you cannot override, and also the member variables you may need to set based on your needs.

Table 20. Differential Equation Function Classifications.

Constructor

The constructor is called automatically when the numerical integration object is instantiated. You must define the number of equations and size the derivative vector. You also need to set the values for any other internal variables, such as for event handling and default values.

Definition

CMyODE();

Usage

CMyODE MyODE;

CMyODE* pMyODE;

pMyODE = new CMyODE;

Arguments

None.

Return Value

None.

Notes

The constructor is automatically called when the ODE object is instantiated. You must define the number of equations and events, and must size the derivative and event vectors (see below). By default, each of the internal variables is set to zero for real and integer variables, and false for all Boolean variables. The constructor is where you must reset any of these values that will be required by your ODE.

ODE Function to Solve

Definition

virtual colvec OdeFcn(const double t, const colvec& y);

Usage

This function is called by the integrator to provide the derivative for the specified time and state vector. This function will be called numerous times in the computation of a single step. It should not set any internal variables which depend on time; wait for the UpdateStep() call to set any of these (for instance, updating discrete filter states).

Arguments

Return Value

Notes

The ODE function returns the derivative of the state vector for the given time and state vector. This function should place the derivative in the mYdot internal variable (this is to save creation and destruction of temporary variables, which in turn saves time) and then return it. You must override this pure virtual function.

Update and Output Functions

Table 21 lists the update and output functions that are called by the integrator. Override these routines if your ODE needs to handle these events.

Table 21. ODE Update and Output Functions.

Definition

virtual void Setup() ;

Usage

pODE->Setup();

Arguments

None.

Return Value

None.

Notes

This empty function is provided to you as a placeholder for code which needs to be called at the start of the simulation run; it is intended primarily for use with Monte Carlo simulations. Override this routine to perform any initialization which needs to be done at the start (reset discrete states, etc.). You must call this routine manually.

Definition

virtual void Cleanup() ;

Usage

pODE-> Cleanup ();

Arguments

None.

Return Value

None.

Notes

This empty function is provided to you as a placeholder for code which needs to be called at the end of the simulation run; it is intended primarily for use with Monte Carlo simulations. Override this routine to perform any cleanup which needs to be done at the end (delete temporary storage, etc.). You must call this routine manually.

Definition

virtual void Reset();

Usage

This function is called by the integrator when a reset occurs (i.e., the integrator Reset() function is called). This allows you to reset any internal variables to prepare the ODE for a new integration.

Arguments

None.

Return Value

None.

Notes

None.

Definition

virtual void UpdateStep(double t, colvec& y, bool fixedStep);

Usage

This function is called by the integrator at the start of every integration step to allow the ODE to update any internal variables it requires (a discrete filter, for example).

Arguments

Return Value

None.

Notes

This function allows you to update any internal states at the start of the new step. The fixedStep variable is provided for such instances where you need to change the ODE behavior based on step size type (for instance, choose between discrete or continuous filter). The flag is true for fixed step mode, false for variable step mode.

Definition

virtual void OutputStep(double t, colvec& y, bool fixedStep)

Usage

This function is called by the integrator at every requested output time. This allows you to save any internal states or perform any output you require.

Arguments

Return Value

None.

Notes

The function allows you to handle output for the ODE internal states. All output will be handled by this function if you call the integrator function LogOutput(false). This routine can be handy when you require a large number of output points such that memory would be prohibitive.

Definition

virtual void OutputEvent(double t, const colvec& y, int state)

Usage

This function is called by the integrator at every event. This allows you to save any internal states or perform any output you require.

Arguments

Return Value

None.

Notes

The function allows you to handle output for the ODE internal states. All event output will be handled by this function if you call the integrator function LogOutput(false). Note that this routine will be called even if you have set OutputEventsInYOut(false). This routine can be handy when you wish to output additional information about the event other than the states.

Event Handling Routines

These routines supply your ODE with the capability to have the integrator look for different types of events during the course of integration. An event is an equation who's value is zero at the instant of the event. You can tell the integrator to look for events as they cross zero going from negative to positive, positive to negative, or in either direction. You also tell the integrator how to handle the event: record the event for informational purposes only; reset the state values and continue integration; or, terminate integration. You can have any number of events for the integrator to look for, each with its own direction and type.

Two enumerated and matrix types are defined for event handling. The type dirvec is a column vector whose elements are of the enumerated type EEventDirection. Table 22 lists the members of the EEventDirection enumerated type and their values. The CDifferentialEquation base class defines the variable mEventDirection of type dirvec for you to define the direction of each event. The type typevec is a column vector whose elements are of the enumerated type EEventType. Table 23 lists the members of the EEventType enumerated type and their values. The CDifferentialEquation base class defines the variable mEventType of type typevec for you to define the direction of each event. You should size and initialize both of these variables in the constructor of your derived class.

Table 22. Values of the EEventDirection enumerated type.

Table 23. Values of the EEventType enumerated type.

At each step, the integrator will query the ODE to get the event states. It looks for zero crossings for each of the event states. When a zero crossing occurs in the correct direction, the integrator then performs a zero search to find the time at which each event occurs. These are then sorted in order and processed one at a time. A terminal or reset event will cause the integrator to halt at that time and reset or halt the integrator as required.

Your event equations can be of any form that returns a value of zero at the event, and has positive and negative values on either side. For instance, you can look for extrema of an oscillator for finding when the velocity crosses zero. You can halt integration at a particular time, t₀, by searching for when (t-t₀) crosses zero. In fact, any equation involving time and the states can be used.

Table 24 lists the event handling routines which you can override to allow your ODE to process events. These variables are listed in the sections to follow.

Table 24. ODE Event Handling Routines.

Definition

virtual colvec& GetEventState(double t, const colvec& y);

virtual double GetEventState(double t, const colvec& y, int state);

Usage

This function is called by the integrator to get the event states for the specified time and state vector. A particular state can be queried by supplying a value for state.

Arguments

Return Value

Notes

You should override the first version of this function to provide the integrator with your event equations (the second version calls the first). In it you must compute the event values and place them in the variable mEventStates, which you then set as the return value. All equations must be of type double.

Definition

virtual void ResetStatesAfterEvent(double t, colvec& y, int state);

Usage

This function is called by the integrator to get the new integrator states after a reset event.

Arguments

Return Value

None.

Notes

Override this function to compute the new states after a reset event. The integrator provides the time of the event and a reference to the current state vector. You must return the new state vector in y. The variable state is provided to let you know which state triggered the event in order that you may reset the states appropriately.

Values of Default Integration Arguments

Table 25 lists the functions which allow you to provide default arguments to the numerical integrator.

Table 25. ODE Functions That Set the Values of Default Integration Arguments.

Definition

virtual double DefaultH0() const;

Usage

This function provides the integrator with a default initial step size.

Arguments

None.

Return Value

Notes

None.

Definition

virtual double DefaultHMax() const;

Usage

This function provides the integrator with a default maximum step size.

Arguments

None.

Return Value

Notes

None.

Definition

virtual double DefaultHMin() const;

Usage

This function provides the integrator with a default minimum step size.

Arguments

None.

Return Value

Notes

None.

Definition

virtual double DefaultAbsTol() const

Usage

This function provides the integrator with a default absolute tolerance.

Arguments

None.

Return Value

Notes

This function allows you to provide a default absolute tolerance. Note that only a single value is returned; this value will be used for all elements of the error vector. If you wish to use different tolerances for each element of the state vector you will need to set them manually with the SetAbsTol() integrator function.

Definition

virtual double DefaultRelTol() const;

Usage

This function provides the integrator with a default relative tolerance.

Arguments

None.

Return Value

Notes

This function allows you to provide a default relative tolerance. Note that only a single value is returned; this value will be used for all elements of the error vector. If you wish to use different tolerances for each element of the state vector you will need to set them manually with the SetRelTol() integrator function.

Definition

virtual colvec DefaultTOut() const;

Usage

This function provides the integrator with a default list of requested output times.

Arguments

None.

Return Value

Notes

None.

Definition

virtual colvec DefaultICs() const;

Usage

This function provides the integrator with a default vector of initial conditions.

Arguments

None.

Return Value

Notes

None.

Exception Handling

When an error occurs in the evaluation of your ODE, you can use the C++ exception mechanism to pass this information on. There is a special purpose exception type defined in CDifferentialEquation for this purpose. This type is defined to be:

class ODEerror : public std::runtime_error

public:

explicit ODEerror(const char* what_arg) : std::runtime_error(what_arg)

explicit ODEerror(const string& what_arg) : std::runtime_error(what_arg)

When an error occurs, you must throw an exception of type ODEerror. For example,

if (y == 0.0)

throw(ODEerror('Divide by zero in MyODERoutine.cpn'));

else

p = x / y;

If the error occurred in a call from the integrator, this exception will be caught and re-thrown as an exception of type CNumericalIntegrator:: IntegratorError (see Section ). This new exception will pass on the text you supplied to ODEerror. If the error occurs in a call you have made directly, then you are responsible for catching and responding to the exception.

try

catch (CDifferentialEquation::ODEerror theErr)

Note that the ODEerror exception returns one value, a string describing the error. You can access this string using the function ODEerror::what() (inherited from std::runtime_error).

You can turn off the exception handling feature by defining the macro "_MSL_NO_EXCEPTIONS" in a global prefix file.

Functions That Cannot Be Overridden

This section describes the functions defined in the base class which will be inherited by your ODE object. These functions are called by the integrator at various times to provide information about the ODE object. While you cannot override these functions, they will reference internal variables that you will need to set in the other routines listed above. The functions are listed below grouped according to purpose.

Internal Variables

This section lists the internal variables defined in the base class. These are the variables which you must modify to prepare the ODE to communicate with the integrator. These variables are listed below, grouped by purpose.

Example Derived ODE Class

const double kCOR = 0.8; // define the coefficient of restitution

const double kG = 9.8; // define gravity (m/s^2)

class CBouncingBallODE : public CDifferentialEquation

public:

CBouncingBallODE (); // constructor

~CBouncingBallODE () // destructor

virtual colvec OdeFcn (const double t, const colvec& y);

virtual colvec& GetEventState (double t, colvec& y);

virtual void ResetStatesAfterEvent (double t, colvec& y, int state);

virtual double DefaultAbsTol () const

virtual double DefaultRelTol () const

CBouncingBallODE::CBouncingBallODE()

colvec CBouncingBallODE::OdeFcn(const double t, const colvec& y)

colvec& CBouncingBallODE::GetEventState(double t, colvec& y)

void CBouncingBallODE::ResetStatesAfterEvent(double t, colvec& y, int state)