This package provides classes to solve Ordinary Differential Equations problems and also compute derivatives of the solution.

This package solves Initial Value Problems of the form `y'=f(t,y,p)` with `t0`, `y(t0)=y0` and `dy(t0)/dp` known. The provided integrator computes estimates of `y(t)`, `dy(t)/dy0` and `dy(t)/dp` from `t=t0` to `t=t1`, where `y` is the state and `p` is a parameters array.

The classes in this package mimic the behavior of classes and interfaces from the org.apache.commons.math.ode, org.apache.commons.math.ode.events and org.apache.commons.math.ode.sampling packages, adding the jacobians `dy(t)/dy0` and `dy(t)/dp` to the methods signatures.

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.

Here is a simple example of usage. We consider a two-dimensional problem where the state vector y is the solution of the ordinary differential equations

- y'0(t) = ω × (c1 - y1(t))
- y'1(t) = ω × (y0(t) - c0)

with some initial state y(t0) = (y0(t0), y1(t0)).

The point trajectory depends on the initial state y(t0) and on the ODE parameter ω. We want to compute both the final point position y(tend) and the sensitivity of this point with respect to the initial state: dy(tend)/dy(t0) which is a 2×2 matrix and its sensitivity with respect to the parameter: dy(tend)/dω which is a 2×1 matrix.

We consider first the simplest implementation: we define only the ODE and let the classes compute the necessary jacobians by itself: ``

    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;
        }

    }

We compute the results we want as follows: ``

        // low level integrator
        FirstOrderIntegrator lowIntegrator =
                new DormandPrince54Integrator(1.0e-8, 100.0, 1.0e-10, 1.0e-10);

        // set up ODE
        double cx = 1.0;
        double cy = 1.0;
        double omega = 0.1;
        ParameterizedODE  ode = new BasicCircleODE(new double[] { cx, cy }, omega);

        // set up high level integrator, using finite differences step hY and hP to compute jacobians
        double[] hY = new double[] { 1.0e-5, 1.0e-5 };
        double[] hP = new double[] { 1.0e-5 };
        FirstOrderIntegratorWithJacobians integrator =
                new FirstOrderIntegratorWithJacobians(lowIntegrator, ode, hY, hP);

        // set up initial state and derivatives
        double     t0    = 0.0;
        double[]   y0    = new double[] { 0.0, cy };
        double[][] dy0dp = new double[2][1] = { { 0.0 }, { 0.0 } }; // y0 does not depend on omega

        // solve problem
        double     t     = Math.PI / (2 * omega);
        double[]   y     = new double[2];
        double[][] dydy0 = new double[2][2];
        double[][] dydp  = new double[2][1];
        integrator.integrate(t0, y0, dy0dp, t, y, dydy0, dydp);

If in addition to getting the end state and its derivatives, we want to print the state throughout integration process, we have to register a step handler. Inserting the following before the call to integrate does the trick: ``

        StpeHandlerWithJacobians stepHandler = new StpeHandlerWithJacobians() {
                public void reset() {}

                public boolean requiresDenseOutput() { return false; }

                public void handleStep(StepInterpolatorWithJacobians interpolator, boolean isLast)
                    throws DerivativeException {
                    double   t = interpolator.getCurrentTime();
                    double[] y = interpolator.getInterpolatedY();
                    System.out.println(t + " " + y[0] + " " + y[1]);
                }
        };
        integrator.addStepHandler(stepHandler);

The implementation above relies on finite differences with small step sizes to compute the internal jacobians. Since the ODE is really simple here, a better way is to compute them exactly. So instead of implementing ParameterizedODE, we implement the ODEWithJacobians interface as follows (i.e. we replace the setParameter method by a computeJacobians method): ``

    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];

        }

    }

With this implementation, the hY and hP arrays are not needed anymore: `

        ODEWithJacobians ode = new EnhancedCircleODE(new double[] { cx, cy }, omega);
        FirstOrderIntegratorWithJacobians integrator =
                new FirstOrderIntegratorWithJacobians(lowIntegrator, ode);

`