commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From l..@apache.org
Subject svn commit: r920131 - in /commons/proper/math/trunk/src: main/java/org/apache/commons/math/ode/ main/java/org/apache/commons/math/ode/jacobians/ site/xdoc/userguide/
Date Sun, 07 Mar 2010 22:19:19 GMT
Author: luc
Date: Sun Mar  7 22:19:18 2010
New Revision: 920131

URL: http://svn.apache.org/viewvc?rev=920131&view=rev
Log:
improved documentation of ODE package, including the new jacobians part

Modified:
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/jacobians/EventHandlerWithJacobians.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/jacobians/StepHandlerWithJacobians.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/jacobians/package.html
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/package.html
    commons/proper/math/trunk/src/site/xdoc/userguide/index.xml
    commons/proper/math/trunk/src/site/xdoc/userguide/ode.xml

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/jacobians/EventHandlerWithJacobians.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/jacobians/EventHandlerWithJacobians.java?rev=920131&r1=920130&r2=920131&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/jacobians/EventHandlerWithJacobians.java
(original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/jacobians/EventHandlerWithJacobians.java
Sun Mar  7 22:19:18 2010
@@ -45,6 +45,27 @@
  * error (this event handling feature is available for all integrators,
  * 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$
  * @since 2.1
  */

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/jacobians/StepHandlerWithJacobians.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/jacobians/StepHandlerWithJacobians.java?rev=920131&r1=920130&r2=920131&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/jacobians/StepHandlerWithJacobians.java
(original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/jacobians/StepHandlerWithJacobians.java
Sun Mar  7 22:19:18 2010
@@ -32,6 +32,27 @@
  * last one, store the points in an ephemeris, or forward them to
  * 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 StepInterpolatorWithJacobians
  * @version $Revision$ $Date$

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/jacobians/package.html
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/jacobians/package.html?rev=920131&r1=920130&r2=920131&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/jacobians/package.html
(original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/jacobians/package.html
Sun Mar  7 22:19:18 2010
@@ -19,7 +19,7 @@
 <body>
 <p>
 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>
@@ -33,5 +33,205 @@
 where <code>y</code> is the state and <code>p</code> is a parameters
 array.
 </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>
 </html>

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/package.html
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/package.html?rev=920131&r1=920130&r2=920131&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/package.html (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/package.html Sun Mar
 7 22:19:18 2010
@@ -27,6 +27,13 @@
 <code>y(t<sub>0</sub>)=y<sub>0</sub></code> known. The
provided
 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>.
+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>
@@ -144,8 +151,17 @@
 <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.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>
 </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>
 </html>

Modified: commons/proper/math/trunk/src/site/xdoc/userguide/index.xml
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/site/xdoc/userguide/index.xml?rev=920131&r1=920130&r2=920131&view=diff
==============================================================================
--- commons/proper/math/trunk/src/site/xdoc/userguide/index.xml (original)
+++ commons/proper/math/trunk/src/site/xdoc/userguide/index.xml Sun Mar  7 22:19:18 2010
@@ -128,9 +128,10 @@
         <li><a href="ode.html">13. Ordinary Differential Equations Integration</a>
                 <ul>
                 <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.3_ODE_Problems">13.3 ODE Problems</a></li>
-                <li><a href="ode.html#a13.4_Integrators">13.4 Integrators</a></li>
+                <li><a href="ode.html#a13.2_Continuous_Output">13.2 Continuous
Output</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_Available_Integrators">13.4 Available
Integrators</a></li>
+                <li><a href="ode.html#a13.5_Derivatives">13.5 Derivatives</a></li>
                 </ul></li>
         <li><a href="genetics.html">14. Genetic Algorithms</a>   
                 <ul>

Modified: commons/proper/math/trunk/src/site/xdoc/userguide/ode.xml
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/site/xdoc/userguide/ode.xml?rev=920131&r1=920130&r2=920131&view=diff
==============================================================================
--- commons/proper/math/trunk/src/site/xdoc/userguide/ode.xml (original)
+++ commons/proper/math/trunk/src/site/xdoc/userguide/ode.xml Sun Mar  7 22:19:18 2010
@@ -67,21 +67,79 @@
           <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
           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>
+        <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>
           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
-          <code>FirstOrderIntegrator.integrate</code> method. The second one
should be used when more in-depth
-          information is needed throughout the integration process. The user can register
an object implementing
-          the <a href="../apidocs/org/apache/commons/math/ode/sampling/StepHandler.html">StepHandler</a>
interface or a
+          <code>FirstOrderIntegrator.integrate</code> method, as shown by previous
example. The second one
+          should be used when more in-depth information is needed throughout the integration
process. The user
+          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 user-specified object implementing the
           <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
           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>
+        <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>
           <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
@@ -114,7 +172,13 @@
           automatic guess is wrong.
         </p>
       </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>
           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>
@@ -125,7 +189,6 @@
           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).
         </p>
-
         <p>
           When an event is triggered, the event time, current state and an indicator
           whether the switching function was increasing or decreasing at event time
@@ -158,7 +221,7 @@
          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
          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
          discontinuity that must be handled appropriately by the integrator. In such a case,
          we would use a switching function setting similar to this:
@@ -187,27 +250,7 @@
 }
         </source>
       </subsection>
-      <subsection name="13.3 ODE Problems" href="problems">
-        <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">
+      <subsection name="13.4 Available Integrators" href="integrators">
         <p>
           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
@@ -238,6 +281,155 @@
           </table>
         </p>
       </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>
   </body>
 </document>



Mime
View raw message