13.1 Overview

The ode package provides classes to solve Ordinary Differential Equations problems.

This package solves Initial Value Problems of the form y'=f(t,y) with t₀ and y(t₀)=y₀ known. The provided integrators compute an estimate of y(t) from t=t₀ to t=t₁.

All integrators provide dense output. This means that besides computing the state vector at discrete times, they also provide a cheap mean to get both the state and its derivative between the time steps. They do so through classes extending the StepInterpolator abstract class, which are made available to the user at the end of each step.

All integrators handle multiple discrete events detection based on switching functions. This means that the integrator can be driven by user specified discrete events (occurring when the sign of user-supplied switching function changes). The steps are shortened as needed to ensure the events occur at step boundaries (even if the integrator is a fixed-step integrator). When the events are triggered, integration can be stopped (this is called a G-stop facility), the state vector can be changed, or integration can simply go on. The latter case is useful to handle discontinuities in the differential equations gracefully and get accurate dense output even close to the discontinuity.

All integrators support setting a maximal number of evaluations of differential equations function. If this number is exceeded, an exception will be thrown during integration. This can be used to prevent infinite loops if for example error control or discrete events create a really large number of extremely small steps. By default, the maximal number of evaluation is set to Integer.MAX_VALUE (i.e. 2³¹-1 or 2147483647). It is recommended to set this maximal number to a value suited to the ODE problem, integration range, and step size or error control settings.

The user should describe his problem in his own classes which should implement the FirstOrderDifferentialEquations interface. Then he should pass it to the integrator he prefers among all the classes that implement the FirstOrderIntegrator interface. The following example shows how to implement the simple two-dimensional problem:

• y'₀(t) = ω × (c₁ - y₁(t))
• y'₁(t) = ω × (y₀(t) - c₀)

private static class CircleODE implements FirstOrderDifferentialEquations {

    private double[] c;
    private double omega;

    public CircleODE(double[] c, double omega) {
        this.c     = c;
        this.omega = omega;
    }

    public int getDimension() {
        return 2;
    }

    public void computeDerivatives(double t, double[] y, double[] yDot) {
        yDot[0] = omega * (c[1] - y[1]);
        yDot[1] = omega * (y[0] - c[0]);
    }

}
        

Computing the state y(16.0) starting from y(0.0) = (0.0, 1.0) and integrating the ODE is done as follows (using Dormand-Prince 8(5,3) integrator as an example):

FirstOrderIntegrator dp853 = new DormandPrince853Integrator(1.0e-8, 100.0, 1.0e-10, 1.0e-10);
FirstOrderDifferentialEquations ode = new CircleODE(new double[] { 1.0, 1.0 }, 0.1);
double[] y = new double[] { 0.0, 1.0 }; // initial state
dp853.integrate(ode, 0.0, y, 16.0, y); // now y contains final state at time t=16.0
        

13.2 Continuous Output

The solution of the integration problem is provided by two means. The first one is aimed towards simple use: the state vector at the end of the integration process is copied in the y array of the FirstOrderIntegrator.integrate method, as shown by previous example. The second one should be used when more in-depth information is needed throughout the integration process. The user can register an object implementing the StepHandler interface or a StepNormalizer object wrapping a user-specified object implementing the FixedStepHandler interface into the integrator before calling the FirstOrderIntegrator.integrate method. The user object will be called appropriately during the integration process, allowing the user to process intermediate results. The default step handler does nothing. Considering again the previous example, we want to print the trajectory of the point to check it really is a circle arc. We simply add the following before the call to integrator.integrate:

StepHandler stepHandler = new StepHandler() {
    public void init(double t0, double[] y0, double t) {
    }
            
    public void handleStep(StepInterpolator interpolator, boolean isLast) {
        double   t = interpolator.getCurrentTime();
        double[] y = interpolator.getInterpolatedState();
        System.out.println(t + " " + y[0] + " " + y[1]);
    }
};
integrator.addStepHandler(stepHandler);
        

ContinuousOutputModel is a special-purpose step handler that is able to store all steps and to provide transparent access to any intermediate result once the integration is over. An important feature of this class is that it implements the Serializable interface. This means that a complete continuous model of the integrated function throughout the integration range can be serialized and reused later (if stored into a persistent medium like a file system or a database) or elsewhere (if sent to another application). Only the result of the integration is stored, there is no reference to the integrated problem by itself.

Other default implementations of the StepHandler interface are available for general needs (DummyStepHandler, StepNormalizer) and custom implementations can be developed for specific needs. As an example, if an application is to be completely driven by the integration process, then most of the application code will be run inside a step handler specific to this application.

Some integrators (the simple ones) use fixed steps that are set at creation time. The more efficient integrators use variable steps that are handled internally in order to control the integration error with respect to a specified accuracy (these integrators extend the AdaptiveStepsizeIntegrator abstract class). In this case, the step handler which is called after each successful step shows up the variable stepsize. The StepNormalizer class can be used to convert the variable stepsize into a fixed stepsize that can be handled by classes implementing the FixedStepHandler interface. Adaptive stepsize integrators can automatically compute the initial stepsize by themselves, however the user can specify it if he prefers to retain full control over the integration or if the automatic guess is wrong.

13.3 Discrete Events Handling

ODE problems are continuous ones. However, sometimes discrete events must be taken into account. The most frequent case is the stop condition of the integrator is not defined by the time t but by a target condition on state y (say y[0] = 1.0 for example).

Discrete events detection is based on switching functions. The user provides a simple g(t, y) function depending on the current time and state. The integrator will monitor the value of the function throughout integration range and will trigger the event when its sign changes. The magnitude of the value is almost irrelevant. For the sake of root finding, it should however be continuous (but not necessarily smooth) at least in the roots vicinity. The steps are shortened as needed to ensure the events occur at step boundaries (even if the integrator is a fixed-step integrator).

When an event is triggered, the event time, current state and an indicator whether the switching function was increasing or decreasing at event time are provided to the user. Several different options are available to him:

• integration can be stopped (this is called a G-stop facility),
• the state vector or the derivatives can be changed,
• or integration can simply go on.

The first case, G-stop, is the most common one. A typical use case is when an ODE must be solved up to some target state is reached, with a known value of the state but an unknown occurrence time. As an example, if we want to monitor a chemical reaction up to some predefined concentration for the first substance, we can use the following switching function setting:

public double g(double t, double[] y) {
  return y[0] - targetConcentration;
}

public int eventOccurred(double t, double[] y, boolean increasing) {
  return STOP;
}
       

The second case, change state vector or derivatives is encountered when dealing with discontinuous dynamical models. A typical case would be the motion of a spacecraft when thrusters are fired for orbital maneuvers. The acceleration is smooth as long as no maneuvers are performed, depending only on gravity, drag, third body attraction, radiation pressure. Firing a thruster introduces a discontinuity that must be handled appropriately by the integrator. In such a case, we would use a switching function setting similar to this:

public double g(double t, double[] y) {
  return (t - tManeuverStart) * (t - tManeuverStop);
}

public int eventOccurred(double t, double[] y, boolean increasing) {
  return RESET_DERIVATIVES;
}
        

The third case is useful mainly for monitoring purposes, a simple example is:

public double g(double t, double[] y) {
  return y[0] - y[1];
}

public int eventOccurred(double t, double[] y, boolean increasing) {
  logger.log("y0(t) and y1(t) curves cross at t = " + t);
  return CONTINUE;
}
        

13.4 Available Integrators

The tables below show the various integrators available for non-stiff problems. Note that the implementations of Adams-Bashforth and Adams-Moulton are adaptive stepsize, not fixed stepsize as is usual for these multi-step integrators. This is due to the fact the implementation relies on the Nordsieck vector representation of the state.

13.5 Derivatives

If in addition to state y(t) the user needs to compute the sensitivity of the state to the initial state or some parameter of the ODE, he will use the classes and interfaces from the org.apache.commons.ode.jacobians package instead of the top level ode package. These classes compute the jacobians dy(t)/dy₀ and dy(t)/dp where y₀ is the initial state and p is some ODE parameter.

The classes and interfaces in this package mimic the behavior of the classes and interfaces of the top level ode package, only adding parameters arrays for the jacobians. The behavior of these classes is to create a compound state vector z containing both the state y(t) and its derivatives dy(t)/dy₀ and dy(t)/dp and to set up an extended problem by adding the equations for the jacobians automatically. These extended state and problems are then provided to a classical underlying integrator chosen by user.

This behavior imply there will be a top level integrator knowing about state and jacobians and a low level integrator knowing only about compound state (which may be big). If the user wants to deal with the top level only, he will use the specialized step handler and event handler classes registered at top level. He can also register classical step handlers and event handlers, but in this case will see the big compound state. This state is guaranteed to contain the original state in the first elements, followed by the jacobian with respect to initial state (in row order), followed by the jacobian with respect to parameters (in row order). If for example the original state dimension is 6 and there are 3 parameters, the compound state will be a 60 elements array. The first 6 elements will be the original state, the next 36 elements will be the jacobian with respect to initial state, and the remaining 18 will be the jacobian with respect to parameters. Dealing with low level step handlers and event handlers is cumbersome if one really needs the jacobians in these methods, but it also prevents many data being copied back and forth between state and jacobians on one side and compound state on the other side.

In order to compute dy(t)/dy₀ and dy(t/dp for any t, the algorithm needs not only the ODE function f such that y'=f(t,y) but also its local jacobians df(t, y, p)/dy and df(t, y, p)/dp.

If the function f is too complex, the user can simply rely on internal differentiation using finite differences to compute these local jacobians. So rather than the FirstOrderDifferentialEquations interface he will implement the ParameterizedODE interface. Considering again our example where only ω is considered a parameter, we get:

public class BasicCircleODE implements ParameterizedODE {

    private double[] c;
    private double omega;

    public BasicCircleODE(double[] c, double omega) {
        this.c     = c;
        this.omega = omega;
    }

    public int getDimension() {
        return 2;
    }

    public void computeDerivatives(double t, double[] y, double[] yDot) {
        yDot[0] = omega * (c[1] - y[1]);
        yDot[1] = omega * (y[0] - c[0]);
    }

    public int getParametersDimension() {
        // we are only interested in the omega parameter
        return 1;
    }

    public void setParameter(int i, double value) {
        omega = value;
    }

}
        

This ODE is provided to the specialized integrator with two arrays specifying the step sizes to use for finite differences (one array for derivation with respect to state y, one array for derivation with respect to parameters p):

double[] hY = new double[] { 0.001, 0.001 };
double[] hP = new double[] { 1.0e-6 };
FirstOrderIntegratorWithJacobians integrator = new FirstOrderIntegratorWithJacobians(dp853, ode, hY, hP);
integrator.integrate(t0, y0, dy0dp, t, y, dydy0, dydp);
        

If the function f is simple, the user can simply provide the local jacobians by himself. So rather than the FirstOrderDifferentialEquations interface he will implement the ODEWithJacobians interface. Considering again our example where only ω is considered a parameter, we get:

public class EnhancedCircleODE implements ODEWithJacobians {

    private double[] c;
    private double omega;

    public EnhancedCircleODE(double[] c, double omega) {
        this.c     = c;
        this.omega = omega;
    }

    public int getDimension() {
        return 2;
    }

    public void computeDerivatives(double t, double[] y, double[] yDot) {
        yDot[0] = omega * (c[1] - y[1]);
        yDot[1] = omega * (y[0] - c[0]);
    }

    public int getParametersDimension() {
        // we are only interested in the omega parameter
        return 1;
    }

    public void computeJacobians(double t, double[] y, double[] yDot, double[][] dFdY, double[][] dFdP) {

        dFdY[0][0] = 0;
        dFdY[0][1] = -omega;
        dFdY[1][0] = omega;
        dFdY[1][1] = 0;

        dFdP[0][0] = 0;
        dFdP[0][1] = omega;
        dFdP[0][2] = c[1] - y[1];
        dFdP[1][0] = -omega;
        dFdP[1][1] = 0;
        dFdP[1][2] = y[0] - c[0];
 
    }

}
        

This ODE is provided to the specialized integrator as is:

FirstOrderIntegratorWithJacobians integrator = new FirstOrderIntegratorWithJacobians(dp853, ode);
integrator.integrate(t0, y0, dy0dp, t, y, dydy0, dydp);