Apache/commons-math
1
0
Fork 0
mirror of https://github.com/apache/commons-math.git synced 2025-02-21 01:15:28 +00:00
Code Issues Packages Projects Releases Wiki Activity

improved documentation of ODE package, including the new jacobians part

Browse Source 
git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@920131 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe 2010-03-07 22:19:18 +00:00
parent 7a7eb1dccc
commit ed40ebc5cd
6 changed files with 484 additions and 33 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

21
src/main/java/org/apache/commons/math/ode/jacobians/EventHandlerWithJacobians.java
View File

@ -45,6 +45,27 @@ import org.apache.commons.math.ode.events.EventException;
* error (this event handling feature is available for all integrators, * error (this event handling feature is available for all integrators,
* including fixed step ones).</p> * including fixed step ones).</p>
* *
* <p>Note that is is possible to register a {@link
* org.apache.commons.math.ode.events.EventHandler classical event handler}
* in the low level integrator used to build a {@link FirstOrderIntegratorWithJacobians}
* rather than implementing this class. The event handlers registered at low level
* will see the big compound state whether the event handlers defined by this interface
* see the original state, and its jacobians in separate arrays.</p>
*
* <p>The compound 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 elements
* will be the jacobian with respect to parameters.</p>
*
* <p>Dealing with low level 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. So for performance
* reasons, it is recommended to use this interface <em>only</em> if jacobians are really
* needed and to use lower level handlers if only state is needed.</p>
*
* @version $Revision$ $Date$ * @version $Revision$ $Date$
* @since 2.1 * @since 2.1
*/ */

21
src/main/java/org/apache/commons/math/ode/jacobians/StepHandlerWithJacobians.java
View File

@ -32,6 +32,27 @@ import org.apache.commons.math.ode.DerivativeException;
* last one, store the points in an ephemeris, or forward them to * last one, store the points in an ephemeris, or forward them to
* specialized processing or output methods.</p> * specialized processing or output methods.</p>
* *
* <p>Note that is is possible to register a {@link
* org.apache.commons.math.ode.sampling.StepHandler classical step handler}
* in the low level integrator used to build a {@link FirstOrderIntegratorWithJacobians}
* rather than implementing this class. The step handlers registered at low level
* will see the big compound state whether the step handlers defined by this interface
* see the original state, and its jacobians in separate arrays.</p>
*
* <p>The compound 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 elements
* will be the jacobian with respect to parameters.</p>
*
* <p>Dealing with low level step 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. So for performance
* reasons, it is recommended to use this interface <em>only</em> if jacobians are really
* needed and to use lower level handlers if only state is needed.</p>
*
* @see FirstOrderIntegratorWithJacobians * @see FirstOrderIntegratorWithJacobians
* @see StepInterpolatorWithJacobians * @see StepInterpolatorWithJacobians
* @version $Revision$ $Date$ * @version $Revision$ $Date$

202
src/main/java/org/apache/commons/math/ode/jacobians/package.html
View File

@ -19,7 +19,7 @@
<body> <body>
<p> <p>
This package provides classes to solve Ordinary Differential Equations problems This package provides classes to solve Ordinary Differential Equations problems
and compute derivatives of the solution. and also compute derivatives of the solution.
</p> </p>
<p> <p>
@ -33,5 +33,205 @@ from <code>t=t<sub>0</sub></code> to <code>t=t<sub>1</sub></code>,
where <code>y</code> is the state and <code>p</code> is a parameters where <code>y</code> is the state and <code>p</code> is a parameters
array. array.
</p> </p>
<p>
The classes in this package mimic the behavior of classes and interfaces from the
<a href="../package-summary.html">org.apache.commons.math.ode</a>,
<a href="../events/package-summary.html">org.apache.commons.math.ode.events</a>
and <a href="../sampling/package-summary.html">org.apache.commons.math.ode.sampling</a>
packages, adding the jacobians <code>dy(t)/dy<sub>0</sub></code> and
<code>dy(t)/dp</code> to the methods signatures.
</p>
<p>
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)/dy<sub>0</sub> and dy(t<sub>0</sub>)/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.
</p>
<p>
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.
</p>
<p>
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
<ul>
<li>y'<sub>0</sub>(t) = &omega; &times; (c<sub>1</sub> - y<sub>1</sub>(t))</li>
<li>y'<sub>1</sub>(t) = &omega; &times; (y<sub>0</sub>(t) - c<sub>0</sub>)</li>
</ul>
with some initial state y(t<sub>0</sub>) = (y<sub>0</sub>(t<sub>0</sub>),
y<sub>1</sub>(t<sub>0</sub>)).
</p>
<p>
The point trajectory depends on the initial state y(t<sub>0</sub>) and on the ODE
parameter &omega;. We want to compute both the final point position y(t<sub>end</sub>)
and the sensitivity of this point with respect to the initial state:
dy(t<sub>end</sub>)/dy(t<sub>0</sub>) which is a 2&times;2 matrix and its sensitivity
with respect to the parameter: dy(t<sub>end</sub>)/d&omega; which is a 2&times;1 matrix.
</p>
<p>
We consider first the simplest implementation: we define only the ODE and let
the classes compute the necessary jacobians by itself:
<code><pre>
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;
}
}
</pre></code>
</p>
<p>
We compute the results we want as follows:
<code><pre>
// 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);
</pre></code>
</p>
<p>
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:
<code><pre>
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);
</pre></code>
</p>
<p>
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):
<code><pre>
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];
}
}
</pre></code>
With this implementation, the hY and hP arrays are not needed anymore:
<code><pre>
ODEWithJacobians ode = new EnhancedCircleODE(new double[] { cx, cy }, omega);
FirstOrderIntegratorWithJacobians integrator =
new FirstOrderIntegratorWithJacobians(lowIntegrator, ode);
</pre></code>
</p>
</body> </body>
</html> </html>

16
src/main/java/org/apache/commons/math/ode/package.html
View File

@ -27,6 +27,13 @@ This package solves Initial Value Problems of the form
<code>y(t<sub>0</sub>)=y<sub>0</sub></code> known. The provided <code>y(t<sub>0</sub>)=y<sub>0</sub></code> known. The provided
integrators compute an estimate of <code>y(t)</code> from integrators compute an estimate of <code>y(t)</code> from
<code>t=t<sub>0</sub></code> to <code>t=t<sub>1</sub></code>. <code>t=t<sub>0</sub></code> to <code>t=t<sub>1</sub></code>.
If in addition to <code>y(t)</code> users need to get the
derivatives with respect to the initial state
<code>dy(t)/dy(t<sub>0</sub>)</code> or the derivatives with
respect to some ODE parameters <code>dy(t)/dp</code>, then the
classes from the <a href="./jacobians/package-summary.html">
org.apache.commons.math.ode.jacobians</a> package must be used
instead of the classes in this package.
</p> </p>
<p> <p>
@ -144,8 +151,17 @@ automatic guess is wrong.
<tr><td>{@link org.apache.commons.math.ode.nonstiff.DormandPrince54Integrator Dormand-Prince 5(4)}</td><td>5</td><td>4</td></tr> <tr><td>{@link org.apache.commons.math.ode.nonstiff.DormandPrince54Integrator Dormand-Prince 5(4)}</td><td>5</td><td>4</td></tr>
<tr><td>{@link org.apache.commons.math.ode.nonstiff.DormandPrince853Integrator Dormand-Prince 8(5,3)}</td><td>8</td><td>5 and 3</td></tr> <tr><td>{@link org.apache.commons.math.ode.nonstiff.DormandPrince853Integrator Dormand-Prince 8(5,3)}</td><td>8</td><td>5 and 3</td></tr>
<tr><td>{@link org.apache.commons.math.ode.nonstiff.GraggBulirschStoerIntegrator Gragg-Bulirsch-Stoer}</td><td>variable (up to 18 by default)</td><td>variable</td></tr> <tr><td>{@link org.apache.commons.math.ode.nonstiff.GraggBulirschStoerIntegrator Gragg-Bulirsch-Stoer}</td><td>variable (up to 18 by default)</td><td>variable</td></tr>
<tr><td>{@link org.apache.commons.math.ode.nonstiff.AdamsBashforthIntegrator Adams-Bashforth}</td><td>variable</td><td>variable</td></tr>
<tr><td>{@link org.apache.commons.math.ode.nonstiff.AdamsMoultonIntegrator Adams-Moulton}</td><td>variable</td><td>variable</td></tr>
</table> </table>
</p> </p>
<p>
In the table above, the {@link org.apache.commons.math.ode.nonstiff.AdamsBashforthIntegrator
Adams-Bashforth} and {@link org.apache.commons.math.ode.nonstiff.AdamsMoultonIntegrator
Adams-Moulton} integrators appear as variable-step ones. This is an experimental extension
to the classical algorithms using the Nordsieck vector representation.
</p>
</body> </body>
</html> </html>

7
src/site/xdoc/userguide/index.xml
View File

@ -128,9 +128,10 @@
<li><a href="ode.html">13. Ordinary Differential Equations Integration</a> <li><a href="ode.html">13. Ordinary Differential Equations Integration</a>
<ul> <ul>
<li><a href="ode.html#a13.1_Overview">13.1 Overview</a></li> <li><a href="ode.html#a13.1_Overview">13.1 Overview</a></li>
<li><a href="ode.html#a13.2_Discrete_Events_Handling">13.2 Discrete Events Handling</a></li> <li><a href="ode.html#a13.2_Continuous_Output">13.2 Continuous Output</a></li>
<li><a href="ode.html#a13.3_ODE_Problems">13.3 ODE Problems</a></li> <li><a href="ode.html#a13.3_Discrete_Events_Handling">13.3 Discrete Events Handling</a></li>
<li><a href="ode.html#a13.4_Integrators">13.4 Integrators</a></li> <li><a href="ode.html#a13.4_Available_Integrators">13.4 Available Integrators</a></li>
<li><a href="ode.html#a13.5_Derivatives">13.5 Derivatives</a></li>
</ul></li> </ul></li>
<li><a href="genetics.html">14. Genetic Algorithms</a> <li><a href="genetics.html">14. Genetic Algorithms</a>
<ul> <ul>

250
src/site/xdoc/userguide/ode.xml
View File

@ -67,21 +67,79 @@
<a href="../apidocs/org/apache/commons/math/ode/FirstOrderDifferentialEquations.html">FirstOrderDifferentialEquations</a> <a href="../apidocs/org/apache/commons/math/ode/FirstOrderDifferentialEquations.html">FirstOrderDifferentialEquations</a>
interface. Then he should pass it to the integrator he prefers among all the classes that implement interface. Then he should pass it to the integrator he prefers among all the classes that implement
the <a href="../apidocs/org/apache/commons/math/ode/FirstOrderIntegrator.html">FirstOrderIntegrator</a> the <a href="../apidocs/org/apache/commons/math/ode/FirstOrderIntegrator.html">FirstOrderIntegrator</a>
interface. interface. The following example shows how to implement the simple two-dimensional problem:
<ul>
<li>y'<sub>0</sub>(t) = &#x3c9; &#xD7; (c<sub>1</sub> - y<sub>1</sub>(t))</li>
<li>y'<sub>1</sub>(t) = &#x3c9; &#xD7; (y<sub>0</sub>(t) - c<sub>0</sub>)</li>
</ul>
with some initial state y(t<sub>0</sub>) = (y<sub>0</sub>(t<sub>0</sub>), y<sub>1</sub>(t<sub>0</sub>)).
In fact, the exact solution of this problem is that y(t<sub>0</sub>) moves along a circle
centered at c = (c<sub>0</sub>, c<sub>1</sub>) with constant angular rate &#x3c9;.
</p> </p>
<source>
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]);
}
}
</source>
<p>
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):
</p>
<source>
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
</source>
</subsection>
<subsection name="13.2 Continuous Output" href="continuous">
<p> <p>
The solution of the integration problem is provided by two means. The first one is aimed towards 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 simple use: the state vector at the end of the integration process is copied in the y array of the
<code>FirstOrderIntegrator.integrate</code> method. The second one should be used when more in-depth <code>FirstOrderIntegrator.integrate</code> method, as shown by previous example. The second one
information is needed throughout the integration process. The user can register an object implementing should be used when more in-depth information is needed throughout the integration process. The user
the <a href="../apidocs/org/apache/commons/math/ode/sampling/StepHandler.html">StepHandler</a> interface or a can register an object implementing the
<a href="../apidocs/org/apache/commons/math/ode/sampling/StepHandler.html">StepHandler</a> interface or a
<a href="../apidocs/org/apache/commons/math/ode/sampling/StepNormalizer.html">StepNormalizer</a> object wrapping <a href="../apidocs/org/apache/commons/math/ode/sampling/StepNormalizer.html">StepNormalizer</a> object wrapping
a user-specified object implementing the a user-specified object implementing the
<a href="../apidocs/org/apache/commons/math/ode/sampling/FixedStepHandler.html">FixedStepHandler</a> interface <a href="../apidocs/org/apache/commons/math/ode/sampling/FixedStepHandler.html">FixedStepHandler</a> interface
into the integrator before calling the <code>FirstOrderIntegrator.integrate</code> method. The user object into the integrator before calling the <code>FirstOrderIntegrator.integrate</code> method. The user object
will be called appropriately during the integration process, allowing the user to process intermediate will be called appropriately during the integration process, allowing the user to process intermediate
results. The default step handler does nothing. 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:
</p> </p>
<source>
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);
</source>
<p> <p>
<a href="../apidocs/org/apache/commons/math/ode/ContinuousOutputModel.html">ContinuousOutputModel</a> <a href="../apidocs/org/apache/commons/math/ode/ContinuousOutputModel.html">ContinuousOutputModel</a>
is a special-purpose step handler that is able to store all steps and to provide transparent access to is a special-purpose step handler that is able to store all steps and to provide transparent access to
@ -114,7 +172,13 @@
automatic guess is wrong. automatic guess is wrong.
</p> </p>
</subsection> </subsection>
<subsection name="13.2 Discrete Events Handling" href="events"> <subsection name="13.3 Discrete Events Handling" href="events">
<p>
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).
</p>
<p> <p>
Discrete events detection is based on switching functions. The user provides Discrete events detection is based on switching functions. The user provides
a simple <a href="../apidocs/org/apache/commons/math/ode/events/EventHandler.html">g(t, y)</a> a simple <a href="../apidocs/org/apache/commons/math/ode/events/EventHandler.html">g(t, y)</a>
@ -125,7 +189,6 @@
root finding. The steps are shortened as needed to ensure the events occur 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). at step boundaries (even if the integrator is a fixed-step integrator).
</p> </p>
<p> <p>
When an event is triggered, the event time, current state and an indicator When an event is triggered, the event time, current state and an indicator
whether the switching function was increasing or decreasing at event time whether the switching function was increasing or decreasing at event time
@ -158,7 +221,7 @@ public int eventOccurred(double t, double[] y, boolean increasing) {
The second case, change state vector or derivatives is encountered when dealing 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 with discontinuous dynamical models. A typical case would be the motion of a
spacecraft when thrusters are fired for orbital maneuvers. The acceleration is spacecraft when thrusters are fired for orbital maneuvers. The acceleration is
smooth as long as no maneuver are performed, depending only on gravity, drag, smooth as long as no maneuvers are performed, depending only on gravity, drag,
third body attraction, radiation pressure. Firing a thruster introduces a third body attraction, radiation pressure. Firing a thruster introduces a
discontinuity that must be handled appropriately by the integrator. In such a case, discontinuity that must be handled appropriately by the integrator. In such a case,
we would use a switching function setting similar to this: we would use a switching function setting similar to this:
@ -187,27 +250,7 @@ public int eventOccurred(double t, double[] y, boolean increasing) {
} }
</source> </source>
</subsection> </subsection>
<subsection name="13.3 ODE Problems" href="problems"> <subsection name="13.4 Available Integrators" href="integrators">
<p>
First order ODE problems are defined by implementing the
<a href="../apidocs/org/apache/commons/math/ode/FirstOrderDifferentialEquations.html">FirstOrderDifferentialEquations</a>
interface before they can be handled by the integrators <code>FirstOrderIntegrator.integrate</code>
method.
</p>
<p>
A first order differential equations problem, as seen by an integrator is the time
derivative <code>dY/dt</code> of a state vector <code>Y</code>, both being one
dimensional arrays. From the integrator point of view, this derivative depends
only on the current time <code>t</code> and on the state vector <code>Y</code>.
</p>
<p>
For real world problems, the derivative depends also on parameters that do not
belong to the state vector (dynamical model constants for example). These
constants are completely outside of the scope of this interface, the classes
that implement it are allowed to handle them as they want.
</p>
</subsection>
<subsection name="13.4 Integrators" href="integrators">
<p> <p>
The tables below show the various integrators available for non-stiff problems. Note that the 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 implementations of Adams-Bashforth and Adams-Moulton are adaptive stepsize, not fixed stepsize
@ -238,6 +281,155 @@ public int eventOccurred(double t, double[] y, boolean increasing) {
</table> </table>
</p> </p>
</subsection> </subsection>
<subsection name="13.5 Derivatives" href="derivatives">
<p>
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 <a
href="../apidocs/org/apache/commons/math/ode/jacobians/package-summary.html">org.apache.commons.ode.jacobians</a>
package instead of the top level ode package. These classes compute the jacobians
dy(t)/dy<sub>0</sub> and dy(t<sub>0</sub>)/dp where y<sub>0</sub> is the initial state
and p is some ODE parameter.
</p>
<p>
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)/dy<sub>0</sub> and dy(t<sub>0</sub>)/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.
</p>
<p>
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.
</p>
<p>
In order to compute dy(t)/dy<sub>0</sub> and dy(t<sub>0</sub>)/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.
</p>
<p>
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 <a
href="../apidocs/org/apache/commons/math/ode/FirstOrderDifferentialEquations.html">FirstOrderDifferentialEquations</a>
interface he will implement the <a
href="../apidocs/org/apache/commons/math/ode/jacobians/ParameterizedODE.html">ParameterizedODE</a>
interface. Considering again our example where only &#x3c9; is considered a parameter, we get:
</p>
<source>
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;
}
}
</source>
<p>
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):
</p>
<source>
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);
</source>
<p>
If the function f is simple, the user can simply provide the local jacobians
by himself. So rather than the <a
href="../apidocs/org/apache/commons/math/ode/FirstOrderDifferentialEquations.html">FirstOrderDifferentialEquations</a>
interface he will implement the <a
href="../apidocs/org/apache/commons/math/ode/jacobians/ODEWithJacobians.html">ODEWithJacobians</a>
interface. Considering again our example where only &#x3c9; is considered a parameter, we get:
</p>
<source>
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];
}
}
</source>
<p>
This ODE is provided to the specialized integrator as is:
</p>
<source>
FirstOrderIntegratorWithJacobians integrator = new FirstOrderIntegratorWithJacobians(dp853, ode);
integrator.integrate(t0, y0, dy0dp, t, y, dydy0, dydp);
</source>
</subsection>
</section> </section>
</body> </body>
</document> </document>