-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 -
-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];
-
- }
-
-}
-
- ODEWithJacobians ode = new EnhancedCircleODE(new double[] { cx, cy }, omega);
- FirstOrderIntegratorWithJacobians integrator =
- new FirstOrderIntegratorWithJacobians(lowIntegrator, ode);
-