Mercurial > hg > de.mpg.mpiwg.itgroup.digilib.core
diff libs/commons-math-2.1/docs/userguide/ode.html @ 10:5f2c5fb36e93
commons-math-2.1 added
| author | dwinter |
|---|---|
| date | Tue, 04 Jan 2011 10:00:53 +0100 |
| parents | |
| children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libs/commons-math-2.1/docs/userguide/ode.html Tue Jan 04 10:00:53 2011 +0100 @@ -0,0 +1,609 @@ +<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd"> + + + + + + + + + + + +<html xmlns="http://www.w3.org/1999/xhtml"> + <head> + <title>Math - The Commons Math User Guide - Ordinary Differential Equations Integration</title> + <style type="text/css" media="all"> + @import url("../css/maven-base.css"); + @import url("../css/maven-theme.css"); + @import url("../css/site.css"); + </style> + <link rel="stylesheet" href="../css/print.css" type="text/css" media="print" /> + <meta http-equiv="Content-Type" content="text/html; charset=ISO-8859-1" /> + </head> + <body class="composite"> + <div id="banner"> + <span id="bannerLeft"> + + Commons Math User Guide + + </span> + <div class="clear"> + <hr/> + </div> + </div> + <div id="breadcrumbs"> + + + + + + + + + <div class="xright"> + + + + + + + + </div> + <div class="clear"> + <hr/> + </div> + </div> + <div id="leftColumn"> + <div id="navcolumn"> + + + + + + + + + <h5>User Guide</h5> + <ul> + + <li class="none"> + <a href="../userguide/index.html">Contents</a> + </li> + + <li class="none"> + <a href="../userguide/overview.html">Overview</a> + </li> + + <li class="none"> + <a href="../userguide/stat.html">Statistics</a> + </li> + + <li class="none"> + <a href="../userguide/random.html">Data Generation</a> + </li> + + <li class="none"> + <a href="../userguide/linear.html">Linear Algebra</a> + </li> + + <li class="none"> + <a href="../userguide/analysis.html">Numerical Analysis</a> + </li> + + <li class="none"> + <a href="../userguide/special.html">Special Functions</a> + </li> + + <li class="none"> + <a href="../userguide/utilities.html">Utilities</a> + </li> + + <li class="none"> + <a href="../userguide/complex.html">Complex Numbers</a> + </li> + + <li class="none"> + <a href="../userguide/distribution.html">Distributions</a> + </li> + + <li class="none"> + <a href="../userguide/fraction.html">Fractions</a> + </li> + + <li class="none"> + <a href="../userguide/transform.html">Transform Methods</a> + </li> + + <li class="none"> + <a href="../userguide/geometry.html">3D Geometry</a> + </li> + + <li class="none"> + <a href="../userguide/optimization.html">Optimization</a> + </li> + + <li class="none"> + <strong>Ordinary Differential Equations</strong> + </li> + + <li class="none"> + <a href="../userguide/genetics.html">Genetic Algorithms</a> + </li> + </ul> + <a href="http://maven.apache.org/" title="Built by Maven" class="poweredBy"> + <img alt="Built by Maven" src="../images/logos/maven-feather.png"></img> + </a> + + + + + + + + + </div> + </div> + <div id="bodyColumn"> + <div id="contentBox"> + <div class="section"><h2><a name="a13_Ordinary_Differential_Equations_Integration"></a>13 Ordinary Differential Equations Integration</h2> +<div class="section"><h3><a name="a13.1_Overview"></a>13.1 Overview</h3> +<p> + The ode package provides classes to solve Ordinary Differential Equations problems. + </p> +<p> + This package solves Initial Value Problems of the form y'=f(t,y) with t<sub>0</sub> + and y(t<sub>0</sub>)=y<sub>0</sub> known. The provided integrators compute an estimate + of y(t) from t=t<sub>0</sub> to t=t<sub>1</sub>. + </p> +<p> + All integrators provide dense output. This means that besides computing the state vector + at discrete times, they also provide a cheap mean to get both the state and its derivative + between the time steps. They do so through classes extending the + <a href="../apidocs/org/apache/commons/math/ode/sampling/StepInterpolator.html">StepInterpolator</a> + abstract class, which are made available to the user at the end of each step. + </p> +<p> + All integrators handle multiple discrete events detection based on switching + functions. This means that the integrator can be driven by user specified discrete events + (occurring when the sign of user-supplied <i>switching function</i> changes). The steps are + shortened as needed to ensure the events occur at step boundaries (even if the integrator + is a fixed-step integrator). When the events are triggered, integration can + be stopped (this is called a G-stop facility), the state vector can be changed, or integration + can simply go on. The latter case is useful to handle discontinuities in the differential + equations gracefully and get accurate dense output even close to the discontinuity. + </p> +<p> + All integrators support setting a maximal number of evaluations of differential + equations function. If this number is exceeded, an exception will be thrown during + integration. This can be used to prevent infinite loops if for example error control or + discrete events create a really large number of extremely small steps. By default, the + maximal number of evaluation is set to <code>Integer.MAX_VALUE</code> (i.e. 2<sup>31</sup>-1 + or 2147483647). It is recommended to set this maximal number to a value suited to the ODE + problem, integration range, and step size or error control settings. + </p> +<p> + The user should describe his problem in his own classes which should implement the + <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. The following example shows how to implement the simple two-dimensional problem: + <ul><li>y'<sub>0</sub>(t) = ? × (c<sub>1</sub> - y<sub>1</sub>(t))</li> +<li>y'<sub>1</sub>(t) = ? × (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) moves along a circle + centered at c = (c<sub>0</sub>, c<sub>1</sub>) with constant angular rate ?. + </p> +<div class="source"><pre> +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]); + } + +} + </pre> +</div> +<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> +<div class="source"><pre> +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 + </pre> +</div> +</div> +<div class="section"><h3><a name="a13.2_Continuous_Output"></a>13.2 Continuous Output</h3> +<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, 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. 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> +<div class="source"><pre> +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); + </pre> +</div> +<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 + any intermediate result once the integration is over. An important feature of this class is that it + implements the <code>Serializable</code> interface. This means that a complete continuous model of the + integrated function throughout the integration range can be serialized and reused later (if stored into + a persistent medium like a file system or a database) or elsewhere (if sent to another application). + Only the result of the integration is stored, there is no reference to the integrated problem by itself. + </p> +<p> + Other default implementations of the <a href="../apidocs/org/apache/commons/math/ode/sampling/StepHandler.html">StepHandler</a> + interface are available for general needs + (<a href="../apidocs/org/apache/commons/math/ode/sampling/DummyStepHandler.html">DummyStepHandler</a>, + <a href="../apidocs/org/apache/commons/math/ode/sampling/StepNormalizer.html">StepNormalizer</a>) and custom + implementations can be developed for specific needs. As an example, if an application is to be + completely driven by the integration process, then most of the application code will be run inside a + step handler specific to this application. + </p> +<p> + Some integrators (the simple ones) use fixed steps that are set at creation time. The more efficient + integrators use variable steps that are handled internally in order to control the integration error + with respect to a specified accuracy (these integrators extend the + <a href="../apidocs/org/apache/commons/math/ode/AdaptiveStepsizeIntegrator.html">AdaptiveStepsizeIntegrator</a> + abstract class). In this case, the step handler which is called after each successful step shows up + the variable stepsize. The <a href="../apidocs/org/apache/commons/math/ode/sampling/StepNormalizer.html">StepNormalizer</a> + class can be used to convert the variable stepsize into a fixed stepsize that can be handled by classes + implementing the <a href="../apidocs/org/apache/commons/math/ode/sampling/FixedStepHandler.html">FixedStepHandler</a> + interface. Adaptive stepsize integrators can automatically compute the initial stepsize by themselves, + however the user can specify it if he prefers to retain full control over the integration or if the + automatic guess is wrong. + </p> +</div> +<div class="section"><h3><a name="a13.3_Discrete_Events_Handling"></a>13.3 Discrete Events Handling</h3> +<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> + function depending on the current time and state. The integrator will monitor + the value of the function throughout integration range and will trigger the + event when its sign changes. The magnitude of the value is almost irrelevant, + it should however be continuous (but not necessarily smooth) for the sake of + 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 + are provided to the user. Several different options are available to him: + </p> +<ul><li>integration can be stopped (this is called a G-stop facility),</li> +<li>the state vector or the derivatives can be changed,</li> +<li>or integration can simply go on.</li> +</ul> +<p> + The first case, G-stop, is the most common one. A typical use case is when an + ODE must be solved up to some target state is reached, with a known value of + the state but an unknown occurrence time. As an example, if we want to monitor + a chemical reaction up to some predefined concentration for the first substance, + we can use the following switching function setting: + </p> +<div class="source"><pre> +public double g(double t, double[] y) { + return y[0] - targetConcentration; +} + +public int eventOccurred(double t, double[] y, boolean increasing) { + return STOP; +} + </pre> +</div> +<p> + 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 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: + </p> +<div class="source"><pre> +public double g(double t, double[] y) { + return (t - tManeuverStart) * (t - tManeuverStop); +} + +public int eventOccurred(double t, double[] y, boolean increasing) { + return RESET_DERIVATIVES; +} + </pre> +</div> +<p> + The third case is useful mainly for monitoring purposes, a simple example is: + </p> +<div class="source"><pre> +public double g(double t, double[] y) { + return y[0] - y[1]; +} + +public int eventOccurred(double t, double[] y, boolean increasing) { + logger.log("y0(t) and y1(t) curves cross at t = " + t); + return CONTINUE; +} + </pre> +</div> +</div> +<div class="section"><h3><a name="a13.4_Available_Integrators"></a>13.4 Available Integrators</h3> +<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 + as is usual for these multi-step integrators. This is due to the fact the implementation relies + on the Nordsieck vector representation of the state. + </p> +<p><table class="bodyTable"><tr class="a"><td><font size="+2">Fixed Step Integrators</font></td> +</tr> +<tr class="b"><font size="+1"><td>Name</td> +<td>Order</td> +</font></tr> +<tr class="a"><td><a href="../apidocs/org/apache/commons/math/ode/nonstiff/EulerIntegrator.html">Euler</a></td> +<td>1</td> +</tr> +<tr class="b"><td><a href="../apidocs/org/apache/commons/math/ode/nonstiff/MidpointIntegrator.html">Midpoint</a></td> +<td>2</td> +</tr> +<tr class="a"><td><a href="../apidocs/org/apache/commons/math/ode/nonstiff/ClassicalRungeKuttaIntegrator.html">Classical Runge-Kutta</a></td> +<td>4</td> +</tr> +<tr class="b"><td><a href="../apidocs/org/apache/commons/math/ode/nonstiff/GillIntegrator.html">Gill</a></td> +<td>4</td> +</tr> +<tr class="a"><td><a href="../apidocs/org/apache/commons/math/ode/nonstiff/ThreeEighthesIntegrator.html">3/8</a></td> +<td>4</td> +</tr> +</table> +</p> +<p><table class="bodyTable"><tr class="b"><td><font size="+2">Adaptive Stepsize Integrators</font></td> +</tr> +<tr class="a"><font size="+1"><td>Name</td> +<td>Integration Order</td> +<td>Error Estimation Order</td> +</font></tr> +<tr class="b"><td><a href="../apidocs/org/apache/commons/math/ode/nonstiff/HighamHall54Integrator.html">Higham and Hall</a></td> +<td>5</td> +<td>4</td> +</tr> +<tr class="a"><td><a href="../apidocs/org/apache/commons/math/ode/nonstiff/DormandPrince54Integrator.html">Dormand-Prince 5(4)</a></td> +<td>5</td> +<td>4</td> +</tr> +<tr class="b"><td><a href="../apidocs/org/apache/commons/math/ode/nonstiff/DormandPrince853Integrator.html">Dormand-Prince 8(5,3)</a></td> +<td>8</td> +<td>5 and 3</td> +</tr> +<tr class="a"><td><a href="../apidocs/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegrator.html">Gragg-Bulirsch-Stoer</a></td> +<td>variable (up to 18 by default)</td> +<td>variable</td> +</tr> +<tr class="b"><td><a href="../apidocs/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.html">Adams-Bashforth</a></td> +<td>variable</td> +<td>variable</td> +</tr> +<tr class="a"><td><a href="../apidocs/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.html">Adams-Moulton</a></td> +<td>variable</td> +<td>variable</td> +</tr> +</table> +</p> +</div> +<div class="section"><h3><a name="a13.5_Derivatives"></a>13.5 Derivatives</h3> +<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 ? is considered a parameter, we get: + </p> +<div class="source"><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> +</div> +<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> +<div class="source"><pre> +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); + </pre> +</div> +<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 ? is considered a parameter, we get: + </p> +<div class="source"><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> +</div> +<p> + This ODE is provided to the specialized integrator as is: + </p> +<div class="source"><pre> +FirstOrderIntegratorWithJacobians integrator = new FirstOrderIntegratorWithJacobians(dp853, ode); +integrator.integrate(t0, y0, dy0dp, t, y, dydy0, dydp); + </pre> +</div> +</div> +</div> + + </div> + </div> + <div class="clear"> + <hr/> + </div> + <div id="footer"> + <div class="xright">© + 2003-2010 + + + + + + + + + + </div> + <div class="clear"> + <hr/> + </div> + </div> + </body> +</html>