Other Commons Math examples (source code examples)
Here is a short list of links related to this Commons Math ode.xml source code file:
The search page
Other Commons Math source code examples at this package level
Click here to learn more about this project
Commons Math example source code file (ode.xml)
The Commons Math ode.xml source code<?xml version="1.0"?>
<!--
Licensed to the Apache Software Foundation (ASF) under one or more
contributor license agreements. See the NOTICE file distributed with
this work for additional information regarding copyright ownership.
The ASF licenses this file to You under the Apache License, Version 2.0
(the "License"); you may not use this file except in compliance with
the License. You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
-->
<?xml-stylesheet type="text/xsl" href="./xdoc.xsl"?>
<!-- $Revision: 920132 $ $Date: 2010-03-07 17:22:24 -0500 (Sun, 07 Mar 2010) $ -->
<document url="ode.html">
<properties>
<title>The Commons Math User Guide - Ordinary Differential Equations Integration
</properties>
<body>
<section name="13 Ordinary Differential Equations Integration">
<subsection name="13.1 Overview" href="overview">
<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
and y(t<sub>0)=y0 known. The provided integrators compute an estimate
of y(t) from t=t<sub>0 to t=t1.
</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
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 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 (i.e. 231-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
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
interface. The following example shows how to implement the simple two-dimensional problem:
<ul>
<li>y'0(t) = ω × (c1 - y1(t))
<li>y'1(t) = ω × (y0(t) - c0)
</ul>
with some initial state y(t<sub>0) = (y0(t0), y1(t0)).
In fact, the exact solution of this problem is that y(t) moves along a circle
centered at c = (c<sub>0, c1) with constant angular rate ω.
</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 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 interface or a
<a href="../apidocs/org/apache/commons/math/ode/sampling/StepNormalizer.html">StepNormalizer object wrapping
a user-specified object implementing the
<a href="../apidocs/org/apache/commons/math/ode/sampling/FixedStepHandler.html">FixedStepHandler interface
into the integrator before calling the <code>FirstOrderIntegrator.integrate 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>
<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
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 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
interface are available for general needs
(<a href="../apidocs/org/apache/commons/math/ode/sampling/DummyStepHandler.html">DummyStepHandler,
<a href="../apidocs/org/apache/commons/math/ode/sampling/StepNormalizer.html">StepNormalizer) 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
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
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
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>
</subsection>
<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)
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>the state vector or the derivatives can be changed,
<li>or integration can simply go on.
</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>
<source>
public double g(double t, double[] y) {
return y[0] - targetConcentration;
}
public int eventOccurred(double t, double[] y, boolean increasing) {
return STOP;
}
</source>
<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>
<source>
public double g(double t, double[] y) {
return (t - tManeuverStart) * (t - tManeuverStop);
}
public int eventOccurred(double t, double[] y, boolean increasing) {
return RESET_DERIVATIVES;
}
</source>
<p>
The third case is useful mainly for monitoring purposes, a simple example is:
</p>
<source>
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;
}
</source>
</subsection>
<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
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 border="1" align="center">
<tr BGCOLOR="#CCCCFF"> | Fixed Step Integrators | ||
| Name | Order | ||
| Euler | 1 | ||
| Midpoint | 2 | ||
| Classical Runge-Kutta | 4 | ||
| Gill | 4 | ||
| 3/8 | 4 | ||
| Adaptive Stepsize Integrators | |||
| Name | Integration Order | Error Estimation Order | |
| Higham and Hall | 5 | 4 | |
| Dormand-Prince 5(4) | 5 | 4 | |
| Dormand-Prince 8(5,3) | 8 | 5 and 3 | |
| Gragg-Bulirsch-Stoer | variable (up to 18 by default) | variable | |
| Adams-Bashforth | variable | variable | |
| Adams-Moulton | variable | variable | |
