org.apache.commons.math.ode.jacobians.package.html Maven / Gradle / Ivy
Show all versions of commons-math Show documentation
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];
}
}
ODEWithJacobians ode = new EnhancedCircleODE(new double[] { cx, cy }, omega);
FirstOrderIntegratorWithJacobians integrator =
new FirstOrderIntegratorWithJacobians(lowIntegrator, ode);