13 Ordinary Differential Equations Integration
+
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 t0
+ and y(t0)=y0 known. The provided integrators compute an estimate
+ of y(t) from t=t0 to t=t1.
+
+
+ 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. 231-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'0(t) = ? × (c1 - y1(t))
+- y'1(t) = ? × (y0(t) - c0)
+
+
+ with some initial state y(t
0) = (y
0(t
0), y
1(t
0)).
+ In fact, the exact solution of this problem is that y(t) moves along a circle
+ centered at c = (c
0, c
1) with constant angular rate ?.
+
+
+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 reset() {}
+
+ public boolean requiresDenseOutput() { return false; }
+
+ public void handleStep(StepInterpolator interpolator, boolean isLast) throws DerivativeException {
+ double t = interpolator.getCurrentTime();
+ double[] y = interpolator.getInterpolatedY();
+ 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,
+ it should however be continuous (but not necessarily smooth) for the sake of
+ root finding. 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.
+
+
| Fixed Step Integrators |
+
+| Name |
+Order |
+
+| Euler |
+1 |
+
+| Midpoint |
+2 |
+
+| Classical Runge-Kutta |
+4 |
+
+| Gill |
+4 |
+
+| 3/8 |
+4 |
+
+
+
+
| Adaptive Stepsize Integrators |
+
+| Name |
+Integration Order |
+Error Estimation Order |
+
+| Higham and Hall |
+5 |
+4 |
+
+| Dormand-Prince 5(4) |
+5 |
+4 |
+
+| Dormand-Prince 8(5,3) |
+8 |
+5 and 3 |
+
+| Gragg-Bulirsch-Stoer |
+variable (up to 18 by default) |
+variable |
+
+| Adams-Bashforth |
+variable |
+variable |
+
+| Adams-Moulton |
+variable |
+variable |
+
+
+
+
+
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)/dy0 and dy(t0)/dp where y0 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)/dy0 and dy(t0)/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)/dy0 and dy(t0)/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);
+
+
+
+
+
+
+
+
+
+
+
+
+
+ 