|
||||||||||
PREV PACKAGE NEXT PACKAGE | FRAMES NO FRAMES |
See:
Description
Interface Summary | |
---|---|
EventHandlerWithJacobians | This interface represents a handler for discrete events triggered during ODE integration. |
ODEWithJacobians | This interface represents first order differential equations with parameters and partial derivatives. |
ParameterizedODE | This interface represents first order differential equations with parameters. |
StepHandlerWithJacobians | This interface represents a handler that should be called after each successful step. |
StepInterpolatorWithJacobians | This interface represents an interpolator over the last step during an ODE integration. |
Class Summary | |
---|---|
FirstOrderIntegratorWithJacobians | This class enhances a first order integrator for differential equations to compute also partial derivatives of the solution with respect to initial state and parameters. |
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):
With this implementation, the hY and hP arrays are not needed anymore:
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);
|
||||||||||
PREV PACKAGE NEXT PACKAGE | FRAMES NO FRAMES |