Using the Numerical Integration Routine

Using the Numerical Integration Routine

To use the numerical integration routine, you must first create a Butcher Table object and an ODE object. You will then pass pointers to these objects to the CNumericalIntegrator constructor when you instantiate it. You then set the initial conditions and set up the integrator to handle events or any other special requirements. The actual integration is performed with a few function calls, during which time you can monitor the progress or check for errors. Finally, you can retrieve the output. Each of these steps is detailed below.

Instantiating the Integrator

The integrator can be instantiated in either of two manners as with any C++ object. The first is to create the object directly, passing the required arguments (see Section ). The second method is to create a pointer to an object of type CNumericalIntegrator and instantiate it using the new function. Both of these methods are demonstrated in Figure 1.

In either case, pointers to the Butcher Table and ODE objects must be included in the instantiation call. Additional parameters may be passed as desired, including (see Section ):

Desired output times

Initial conditions

Initial step size

Absolute tolerance

Relative tolerance

Any of the parameters not set in the call to the constructor can be set using a variety of accessor functions described below.

Once you have created the integrator object, you are ready to set it up to perform the integration.

Setting up for Integration

A large number of accessor functions have been provided to prepare the integrator for use. These include such things as the initial conditions, desired output times, step size and norm control, event handling, and values to limit the integrator should a problem arise. See the sections below for descriptions of these functions. There are a few important issues that will be addressed here, however.

The default is to perform integration using variable step sizes. Methods which have an embedded formula (e.g., DOPRI5) will use this formula for step size control. Other methods (e.g., Euler) which do not have an embedded formula require the use of Richardson Extrapolation (see Section for details). Richardson Extrapolation is more accurate but is a factor of 3 more computationally expensive than the embedded formula. You may force an embedded method to use Richardson Extrapolation.

If you require constant step size integration then you must use UseFixedStepSize( to set this. In this case, each step will have the same length, with no guarantee of the error tolerance of the solution. The step size can be set at any time, but should generally not be reset after integration has begun.

When using variable step size control, the size of each step is controlled by the size of the error norm. The integrator can use either the -norm, which monitors the largest error of the solution, or the two-norm, which monitors the root-sum-squared error of the solution. The default is to use the -norm.

You must supply a list of times at which to output the solution. Generally, there are at least two times in this list (the start and finish times). Any number of additional times may be supplied in between, as long as the list is monotonic. The integrator is capable of integrating forward or backward in time; the direction is sensed by the direction of the values requested. You are also allowed to specify a starting time separate from the list. This time must be consistent with the request list and outside its interval. This allows you to start the integration at some arbitrary time and integrate for some time before logging output. This is useful, for instance, where there is a transient solution you wish to have settled before recording output. You can also use this to start integration at one time, integrate to the time specified in the output list, and report only the solution at that time.

Generally, the integrator will keep track of the value of the solution at each of the output times. After integration is finished these values are read using GetYOut( . However, this may not be desirable because of memory limitations. In this case, you can turn off output logging using LogOutput(false), and use the OutputStep() function defined in the ODE object (see Section ), which will then be responsible for output.

If event handling is required for the ODE, then it must be turned on before integration starts. See Sections , , and for details on events.

Performing the Integration

Once the integrator is set up, performing the integration is straightforward. You may proceed by individual integration steps (TakeOneStep), to the next specified output time (StepToNextTOut), to the final specified time (Integrate or StepToTFinal), or by any combination of these. This provides you with control of how integration proceeds, and allows you to monitor the output in order to change or reset the ODE.

The CDifferentialEquation class also defines a number of functions which are called by the integrator to reset or update the ODE after each step, at events, or when output occurs. See Section for details.

Monitoring the Progress and Detecting Errors

You can monitor the progress of the integration when using TakeOneStep( or StepToNextTOut() using GetPercentComplete(), GetFractionComplete(), and IntegrationComplete(). The first two functions tell you how much to the integration has been completed, while the last reports if a terminal event, the final time, or an integration error has been reached.

Additionally, there are a number of functions that provide an indication of how much effort has gone into the integration. There are functions to tell you how many steps have been taken, how many times the derivative function has been called, and let you access the current state of the integrator.

Use IntegratorHasBadInput( or IntegratorHasError() to see if an error has occurred because of illegal input or a problem with integration, respectively. Should either of these indicate that a problem exists, you can inquire to see if a particular error has been set or get a string describing the nature of the problem. Note that if the integrator receives illegal input (such as an inconsistent time vector, for instance), it will raise a flag but will not stop integration. An integrator error, on the other hand, will throw an exception and cease execution until the problem is corrected or the integrator reset.

Getting the Output

Once integration is complete (or at any point along the way) you can retrieve the results using the output accessor functions. These functions can give you the state of the system, the time vector, when events occurred, what the state is at those times, and which event occurred. Note that the integrator calls the ODE's OutputStep( function for every requested output time so it can output or save any values it requires.