|
10
|
1 <!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">
|
|
|
2
|
|
|
3
|
|
|
4
|
|
|
5
|
|
|
6
|
|
|
7
|
|
|
8
|
|
|
9
|
|
|
10
|
|
|
11
|
|
|
12
|
|
|
13 <html xmlns="http://www.w3.org/1999/xhtml">
|
|
|
14 <head>
|
|
|
15 <title>Math - The Commons Math User Guide - Ordinary Differential Equations Integration</title>
|
|
|
16 <style type="text/css" media="all">
|
|
|
17 @import url("../css/maven-base.css");
|
|
|
18 @import url("../css/maven-theme.css");
|
|
|
19 @import url("../css/site.css");
|
|
|
20 </style>
|
|
|
21 <link rel="stylesheet" href="../css/print.css" type="text/css" media="print" />
|
|
|
22 <meta http-equiv="Content-Type" content="text/html; charset=ISO-8859-1" />
|
|
|
23 </head>
|
|
|
24 <body class="composite">
|
|
|
25 <div id="banner">
|
|
|
26 <span id="bannerLeft">
|
|
|
27
|
|
|
28 Commons Math User Guide
|
|
|
29
|
|
|
30 </span>
|
|
|
31 <div class="clear">
|
|
|
32 <hr/>
|
|
|
33 </div>
|
|
|
34 </div>
|
|
|
35 <div id="breadcrumbs">
|
|
|
36
|
|
|
37
|
|
|
38
|
|
|
39
|
|
|
40
|
|
|
41
|
|
|
42
|
|
|
43
|
|
|
44 <div class="xright">
|
|
|
45
|
|
|
46
|
|
|
47
|
|
|
48
|
|
|
49
|
|
|
50
|
|
|
51
|
|
|
52 </div>
|
|
|
53 <div class="clear">
|
|
|
54 <hr/>
|
|
|
55 </div>
|
|
|
56 </div>
|
|
|
57 <div id="leftColumn">
|
|
|
58 <div id="navcolumn">
|
|
|
59
|
|
|
60
|
|
|
61
|
|
|
62
|
|
|
63
|
|
|
64
|
|
|
65
|
|
|
66
|
|
|
67 <h5>User Guide</h5>
|
|
|
68 <ul>
|
|
|
69
|
|
|
70 <li class="none">
|
|
|
71 <a href="../userguide/index.html">Contents</a>
|
|
|
72 </li>
|
|
|
73
|
|
|
74 <li class="none">
|
|
|
75 <a href="../userguide/overview.html">Overview</a>
|
|
|
76 </li>
|
|
|
77
|
|
|
78 <li class="none">
|
|
|
79 <a href="../userguide/stat.html">Statistics</a>
|
|
|
80 </li>
|
|
|
81
|
|
|
82 <li class="none">
|
|
|
83 <a href="../userguide/random.html">Data Generation</a>
|
|
|
84 </li>
|
|
|
85
|
|
|
86 <li class="none">
|
|
|
87 <a href="../userguide/linear.html">Linear Algebra</a>
|
|
|
88 </li>
|
|
|
89
|
|
|
90 <li class="none">
|
|
|
91 <a href="../userguide/analysis.html">Numerical Analysis</a>
|
|
|
92 </li>
|
|
|
93
|
|
|
94 <li class="none">
|
|
|
95 <a href="../userguide/special.html">Special Functions</a>
|
|
|
96 </li>
|
|
|
97
|
|
|
98 <li class="none">
|
|
|
99 <a href="../userguide/utilities.html">Utilities</a>
|
|
|
100 </li>
|
|
|
101
|
|
|
102 <li class="none">
|
|
|
103 <a href="../userguide/complex.html">Complex Numbers</a>
|
|
|
104 </li>
|
|
|
105
|
|
|
106 <li class="none">
|
|
|
107 <a href="../userguide/distribution.html">Distributions</a>
|
|
|
108 </li>
|
|
|
109
|
|
|
110 <li class="none">
|
|
|
111 <a href="../userguide/fraction.html">Fractions</a>
|
|
|
112 </li>
|
|
|
113
|
|
|
114 <li class="none">
|
|
|
115 <a href="../userguide/transform.html">Transform Methods</a>
|
|
|
116 </li>
|
|
|
117
|
|
|
118 <li class="none">
|
|
|
119 <a href="../userguide/geometry.html">3D Geometry</a>
|
|
|
120 </li>
|
|
|
121
|
|
|
122 <li class="none">
|
|
|
123 <a href="../userguide/optimization.html">Optimization</a>
|
|
|
124 </li>
|
|
|
125
|
|
|
126 <li class="none">
|
|
|
127 <strong>Ordinary Differential Equations</strong>
|
|
|
128 </li>
|
|
|
129
|
|
|
130 <li class="none">
|
|
|
131 <a href="../userguide/genetics.html">Genetic Algorithms</a>
|
|
|
132 </li>
|
|
|
133 </ul>
|
|
|
134 <a href="http://maven.apache.org/" title="Built by Maven" class="poweredBy">
|
|
|
135 <img alt="Built by Maven" src="../images/logos/maven-feather.png"></img>
|
|
|
136 </a>
|
|
|
137
|
|
|
138
|
|
|
139
|
|
|
140
|
|
|
141
|
|
|
142
|
|
|
143
|
|
|
144
|
|
|
145 </div>
|
|
|
146 </div>
|
|
|
147 <div id="bodyColumn">
|
|
|
148 <div id="contentBox">
|
|
|
149 <div class="section"><h2><a name="a13_Ordinary_Differential_Equations_Integration"></a>13 Ordinary Differential Equations Integration</h2>
|
|
|
150 <div class="section"><h3><a name="a13.1_Overview"></a>13.1 Overview</h3>
|
|
|
151 <p>
|
|
|
152 The ode package provides classes to solve Ordinary Differential Equations problems.
|
|
|
153 </p>
|
|
|
154 <p>
|
|
|
155 This package solves Initial Value Problems of the form y'=f(t,y) with t<sub>0</sub>
|
|
|
156 and y(t<sub>0</sub>)=y<sub>0</sub> known. The provided integrators compute an estimate
|
|
|
157 of y(t) from t=t<sub>0</sub> to t=t<sub>1</sub>.
|
|
|
158 </p>
|
|
|
159 <p>
|
|
|
160 All integrators provide dense output. This means that besides computing the state vector
|
|
|
161 at discrete times, they also provide a cheap mean to get both the state and its derivative
|
|
|
162 between the time steps. They do so through classes extending the
|
|
|
163 <a href="../apidocs/org/apache/commons/math/ode/sampling/StepInterpolator.html">StepInterpolator</a>
|
|
|
164 abstract class, which are made available to the user at the end of each step.
|
|
|
165 </p>
|
|
|
166 <p>
|
|
|
167 All integrators handle multiple discrete events detection based on switching
|
|
|
168 functions. This means that the integrator can be driven by user specified discrete events
|
|
|
169 (occurring when the sign of user-supplied <i>switching function</i> changes). The steps are
|
|
|
170 shortened as needed to ensure the events occur at step boundaries (even if the integrator
|
|
|
171 is a fixed-step integrator). When the events are triggered, integration can
|
|
|
172 be stopped (this is called a G-stop facility), the state vector can be changed, or integration
|
|
|
173 can simply go on. The latter case is useful to handle discontinuities in the differential
|
|
|
174 equations gracefully and get accurate dense output even close to the discontinuity.
|
|
|
175 </p>
|
|
|
176 <p>
|
|
|
177 All integrators support setting a maximal number of evaluations of differential
|
|
|
178 equations function. If this number is exceeded, an exception will be thrown during
|
|
|
179 integration. This can be used to prevent infinite loops if for example error control or
|
|
|
180 discrete events create a really large number of extremely small steps. By default, the
|
|
|
181 maximal number of evaluation is set to <code>Integer.MAX_VALUE</code> (i.e. 2<sup>31</sup>-1
|
|
|
182 or 2147483647). It is recommended to set this maximal number to a value suited to the ODE
|
|
|
183 problem, integration range, and step size or error control settings.
|
|
|
184 </p>
|
|
|
185 <p>
|
|
|
186 The user should describe his problem in his own classes which should implement the
|
|
|
187 <a href="../apidocs/org/apache/commons/math/ode/FirstOrderDifferentialEquations.html">FirstOrderDifferentialEquations</a>
|
|
|
188 interface. Then he should pass it to the integrator he prefers among all the classes that implement
|
|
|
189 the <a href="../apidocs/org/apache/commons/math/ode/FirstOrderIntegrator.html">FirstOrderIntegrator</a>
|
|
|
190 interface. The following example shows how to implement the simple two-dimensional problem:
|
|
|
191 <ul><li>y'<sub>0</sub>(t) = ? × (c<sub>1</sub> - y<sub>1</sub>(t))</li>
|
|
|
192 <li>y'<sub>1</sub>(t) = ? × (y<sub>0</sub>(t) - c<sub>0</sub>)</li>
|
|
|
193 </ul>
|
|
|
194
|
|
|
195 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>)).
|
|
|
196 In fact, the exact solution of this problem is that y(t) moves along a circle
|
|
|
197 centered at c = (c<sub>0</sub>, c<sub>1</sub>) with constant angular rate ?.
|
|
|
198 </p>
|
|
|
199 <div class="source"><pre>
|
|
|
200 private static class CircleODE implements FirstOrderDifferentialEquations {
|
|
|
201
|
|
|
202 private double[] c;
|
|
|
203 private double omega;
|
|
|
204
|
|
|
205 public CircleODE(double[] c, double omega) {
|
|
|
206 this.c = c;
|
|
|
207 this.omega = omega;
|
|
|
208 }
|
|
|
209
|
|
|
210 public int getDimension() {
|
|
|
211 return 2;
|
|
|
212 }
|
|
|
213
|
|
|
214 public void computeDerivatives(double t, double[] y, double[] yDot) {
|
|
|
215 yDot[0] = omega * (c[1] - y[1]);
|
|
|
216 yDot[1] = omega * (y[0] - c[0]);
|
|
|
217 }
|
|
|
218
|
|
|
219 }
|
|
|
220 </pre>
|
|
|
221 </div>
|
|
|
222 <p>
|
|
|
223 Computing the state y(16.0) starting from y(0.0) = (0.0, 1.0) and integrating the ODE
|
|
|
224 is done as follows (using Dormand-Prince 8(5,3) integrator as an example):
|
|
|
225 </p>
|
|
|
226 <div class="source"><pre>
|
|
|
227 FirstOrderIntegrator dp853 = new DormandPrince853Integrator(1.0e-8, 100.0, 1.0e-10, 1.0e-10);
|
|
|
228 FirstOrderDifferentialEquations ode = new CircleODE(new double[] { 1.0, 1.0 }, 0.1);
|
|
|
229 double[] y = new double[] { 0.0, 1.0 }; // initial state
|
|
|
230 dp853.integrate(ode, 0.0, y, 16.0, y); // now y contains final state at time t=16.0
|
|
|
231 </pre>
|
|
|
232 </div>
|
|
|
233 </div>
|
|
|
234 <div class="section"><h3><a name="a13.2_Continuous_Output"></a>13.2 Continuous Output</h3>
|
|
|
235 <p>
|
|
|
236 The solution of the integration problem is provided by two means. The first one is aimed towards
|
|
|
237 simple use: the state vector at the end of the integration process is copied in the y array of the
|
|
|
238 <code>FirstOrderIntegrator.integrate</code> method, as shown by previous example. The second one
|
|
|
239 should be used when more in-depth information is needed throughout the integration process. The user
|
|
|
240 can register an object implementing the
|
|
|
241 <a href="../apidocs/org/apache/commons/math/ode/sampling/StepHandler.html">StepHandler</a> interface or a
|
|
|
242 <a href="../apidocs/org/apache/commons/math/ode/sampling/StepNormalizer.html">StepNormalizer</a> object wrapping
|
|
|
243 a user-specified object implementing the
|
|
|
244 <a href="../apidocs/org/apache/commons/math/ode/sampling/FixedStepHandler.html">FixedStepHandler</a> interface
|
|
|
245 into the integrator before calling the <code>FirstOrderIntegrator.integrate</code> method. The user object
|
|
|
246 will be called appropriately during the integration process, allowing the user to process intermediate
|
|
|
247 results. The default step handler does nothing. Considering again the previous example, we want to print the
|
|
|
248 trajectory of the point to check it really is a circle arc. We simply add the following before the call
|
|
|
249 to integrator.integrate:
|
|
|
250 </p>
|
|
|
251 <div class="source"><pre>
|
|
|
252 StepHandler stepHandler = new StepHandler() {
|
|
|
253 public void reset() {}
|
|
|
254
|
|
|
255 public boolean requiresDenseOutput() { return false; }
|
|
|
256
|
|
|
257 public void handleStep(StepInterpolator interpolator, boolean isLast) throws DerivativeException {
|
|
|
258 double t = interpolator.getCurrentTime();
|
|
|
259 double[] y = interpolator.getInterpolatedY();
|
|
|
260 System.out.println(t + " " + y[0] + " " + y[1]);
|
|
|
261 }
|
|
|
262 };
|
|
|
263 integrator.addStepHandler(stepHandler);
|
|
|
264 </pre>
|
|
|
265 </div>
|
|
|
266 <p><a href="../apidocs/org/apache/commons/math/ode/ContinuousOutputModel.html">ContinuousOutputModel</a>
|
|
|
267 is a special-purpose step handler that is able to store all steps and to provide transparent access to
|
|
|
268 any intermediate result once the integration is over. An important feature of this class is that it
|
|
|
269 implements the <code>Serializable</code> interface. This means that a complete continuous model of the
|
|
|
270 integrated function throughout the integration range can be serialized and reused later (if stored into
|
|
|
271 a persistent medium like a file system or a database) or elsewhere (if sent to another application).
|
|
|
272 Only the result of the integration is stored, there is no reference to the integrated problem by itself.
|
|
|
273 </p>
|
|
|
274 <p>
|
|
|
275 Other default implementations of the <a href="../apidocs/org/apache/commons/math/ode/sampling/StepHandler.html">StepHandler</a>
|
|
|
276 interface are available for general needs
|
|
|
277 (<a href="../apidocs/org/apache/commons/math/ode/sampling/DummyStepHandler.html">DummyStepHandler</a>,
|
|
|
278 <a href="../apidocs/org/apache/commons/math/ode/sampling/StepNormalizer.html">StepNormalizer</a>) and custom
|
|
|
279 implementations can be developed for specific needs. As an example, if an application is to be
|
|
|
280 completely driven by the integration process, then most of the application code will be run inside a
|
|
|
281 step handler specific to this application.
|
|
|
282 </p>
|
|
|
283 <p>
|
|
|
284 Some integrators (the simple ones) use fixed steps that are set at creation time. The more efficient
|
|
|
285 integrators use variable steps that are handled internally in order to control the integration error
|
|
|
286 with respect to a specified accuracy (these integrators extend the
|
|
|
287 <a href="../apidocs/org/apache/commons/math/ode/AdaptiveStepsizeIntegrator.html">AdaptiveStepsizeIntegrator</a>
|
|
|
288 abstract class). In this case, the step handler which is called after each successful step shows up
|
|
|
289 the variable stepsize. The <a href="../apidocs/org/apache/commons/math/ode/sampling/StepNormalizer.html">StepNormalizer</a>
|
|
|
290 class can be used to convert the variable stepsize into a fixed stepsize that can be handled by classes
|
|
|
291 implementing the <a href="../apidocs/org/apache/commons/math/ode/sampling/FixedStepHandler.html">FixedStepHandler</a>
|
|
|
292 interface. Adaptive stepsize integrators can automatically compute the initial stepsize by themselves,
|
|
|
293 however the user can specify it if he prefers to retain full control over the integration or if the
|
|
|
294 automatic guess is wrong.
|
|
|
295 </p>
|
|
|
296 </div>
|
|
|
297 <div class="section"><h3><a name="a13.3_Discrete_Events_Handling"></a>13.3 Discrete Events Handling</h3>
|
|
|
298 <p>
|
|
|
299 ODE problems are continuous ones. However, sometimes discrete events must be
|
|
|
300 taken into account. The most frequent case is the stop condition of the integrator
|
|
|
301 is not defined by the time t but by a target condition on state y (say y[0] = 1.0
|
|
|
302 for example).
|
|
|
303 </p>
|
|
|
304 <p>
|
|
|
305 Discrete events detection is based on switching functions. The user provides
|
|
|
306 a simple <a href="../apidocs/org/apache/commons/math/ode/events/EventHandler.html">g(t, y)</a>
|
|
|
307 function depending on the current time and state. The integrator will monitor
|
|
|
308 the value of the function throughout integration range and will trigger the
|
|
|
309 event when its sign changes. The magnitude of the value is almost irrelevant,
|
|
|
310 it should however be continuous (but not necessarily smooth) for the sake of
|
|
|
311 root finding. The steps are shortened as needed to ensure the events occur
|
|
|
312 at step boundaries (even if the integrator is a fixed-step integrator).
|
|
|
313 </p>
|
|
|
314 <p>
|
|
|
315 When an event is triggered, the event time, current state and an indicator
|
|
|
316 whether the switching function was increasing or decreasing at event time
|
|
|
317 are provided to the user. Several different options are available to him:
|
|
|
318 </p>
|
|
|
319 <ul><li>integration can be stopped (this is called a G-stop facility),</li>
|
|
|
320 <li>the state vector or the derivatives can be changed,</li>
|
|
|
321 <li>or integration can simply go on.</li>
|
|
|
322 </ul>
|
|
|
323 <p>
|
|
|
324 The first case, G-stop, is the most common one. A typical use case is when an
|
|
|
325 ODE must be solved up to some target state is reached, with a known value of
|
|
|
326 the state but an unknown occurrence time. As an example, if we want to monitor
|
|
|
327 a chemical reaction up to some predefined concentration for the first substance,
|
|
|
328 we can use the following switching function setting:
|
|
|
329 </p>
|
|
|
330 <div class="source"><pre>
|
|
|
331 public double g(double t, double[] y) {
|
|
|
332 return y[0] - targetConcentration;
|
|
|
333 }
|
|
|
334
|
|
|
335 public int eventOccurred(double t, double[] y, boolean increasing) {
|
|
|
336 return STOP;
|
|
|
337 }
|
|
|
338 </pre>
|
|
|
339 </div>
|
|
|
340 <p>
|
|
|
341 The second case, change state vector or derivatives is encountered when dealing
|
|
|
342 with discontinuous dynamical models. A typical case would be the motion of a
|
|
|
343 spacecraft when thrusters are fired for orbital maneuvers. The acceleration is
|
|
|
344 smooth as long as no maneuvers are performed, depending only on gravity, drag,
|
|
|
345 third body attraction, radiation pressure. Firing a thruster introduces a
|
|
|
346 discontinuity that must be handled appropriately by the integrator. In such a case,
|
|
|
347 we would use a switching function setting similar to this:
|
|
|
348 </p>
|
|
|
349 <div class="source"><pre>
|
|
|
350 public double g(double t, double[] y) {
|
|
|
351 return (t - tManeuverStart) * (t - tManeuverStop);
|
|
|
352 }
|
|
|
353
|
|
|
354 public int eventOccurred(double t, double[] y, boolean increasing) {
|
|
|
355 return RESET_DERIVATIVES;
|
|
|
356 }
|
|
|
357 </pre>
|
|
|
358 </div>
|
|
|
359 <p>
|
|
|
360 The third case is useful mainly for monitoring purposes, a simple example is:
|
|
|
361 </p>
|
|
|
362 <div class="source"><pre>
|
|
|
363 public double g(double t, double[] y) {
|
|
|
364 return y[0] - y[1];
|
|
|
365 }
|
|
|
366
|
|
|
367 public int eventOccurred(double t, double[] y, boolean increasing) {
|
|
|
368 logger.log("y0(t) and y1(t) curves cross at t = " + t);
|
|
|
369 return CONTINUE;
|
|
|
370 }
|
|
|
371 </pre>
|
|
|
372 </div>
|
|
|
373 </div>
|
|
|
374 <div class="section"><h3><a name="a13.4_Available_Integrators"></a>13.4 Available Integrators</h3>
|
|
|
375 <p>
|
|
|
376 The tables below show the various integrators available for non-stiff problems. Note that the
|
|
|
377 implementations of Adams-Bashforth and Adams-Moulton are adaptive stepsize, not fixed stepsize
|
|
|
378 as is usual for these multi-step integrators. This is due to the fact the implementation relies
|
|
|
379 on the Nordsieck vector representation of the state.
|
|
|
380 </p>
|
|
|
381 <p><table class="bodyTable"><tr class="a"><td><font size="+2">Fixed Step Integrators</font></td>
|
|
|
382 </tr>
|
|
|
383 <tr class="b"><font size="+1"><td>Name</td>
|
|
|
384 <td>Order</td>
|
|
|
385 </font></tr>
|
|
|
386 <tr class="a"><td><a href="../apidocs/org/apache/commons/math/ode/nonstiff/EulerIntegrator.html">Euler</a></td>
|
|
|
387 <td>1</td>
|
|
|
388 </tr>
|
|
|
389 <tr class="b"><td><a href="../apidocs/org/apache/commons/math/ode/nonstiff/MidpointIntegrator.html">Midpoint</a></td>
|
|
|
390 <td>2</td>
|
|
|
391 </tr>
|
|
|
392 <tr class="a"><td><a href="../apidocs/org/apache/commons/math/ode/nonstiff/ClassicalRungeKuttaIntegrator.html">Classical Runge-Kutta</a></td>
|
|
|
393 <td>4</td>
|
|
|
394 </tr>
|
|
|
395 <tr class="b"><td><a href="../apidocs/org/apache/commons/math/ode/nonstiff/GillIntegrator.html">Gill</a></td>
|
|
|
396 <td>4</td>
|
|
|
397 </tr>
|
|
|
398 <tr class="a"><td><a href="../apidocs/org/apache/commons/math/ode/nonstiff/ThreeEighthesIntegrator.html">3/8</a></td>
|
|
|
399 <td>4</td>
|
|
|
400 </tr>
|
|
|
401 </table>
|
|
|
402 </p>
|
|
|
403 <p><table class="bodyTable"><tr class="b"><td><font size="+2">Adaptive Stepsize Integrators</font></td>
|
|
|
404 </tr>
|
|
|
405 <tr class="a"><font size="+1"><td>Name</td>
|
|
|
406 <td>Integration Order</td>
|
|
|
407 <td>Error Estimation Order</td>
|
|
|
408 </font></tr>
|
|
|
409 <tr class="b"><td><a href="../apidocs/org/apache/commons/math/ode/nonstiff/HighamHall54Integrator.html">Higham and Hall</a></td>
|
|
|
410 <td>5</td>
|
|
|
411 <td>4</td>
|
|
|
412 </tr>
|
|
|
413 <tr class="a"><td><a href="../apidocs/org/apache/commons/math/ode/nonstiff/DormandPrince54Integrator.html">Dormand-Prince 5(4)</a></td>
|
|
|
414 <td>5</td>
|
|
|
415 <td>4</td>
|
|
|
416 </tr>
|
|
|
417 <tr class="b"><td><a href="../apidocs/org/apache/commons/math/ode/nonstiff/DormandPrince853Integrator.html">Dormand-Prince 8(5,3)</a></td>
|
|
|
418 <td>8</td>
|
|
|
419 <td>5 and 3</td>
|
|
|
420 </tr>
|
|
|
421 <tr class="a"><td><a href="../apidocs/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegrator.html">Gragg-Bulirsch-Stoer</a></td>
|
|
|
422 <td>variable (up to 18 by default)</td>
|
|
|
423 <td>variable</td>
|
|
|
424 </tr>
|
|
|
425 <tr class="b"><td><a href="../apidocs/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.html">Adams-Bashforth</a></td>
|
|
|
426 <td>variable</td>
|
|
|
427 <td>variable</td>
|
|
|
428 </tr>
|
|
|
429 <tr class="a"><td><a href="../apidocs/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.html">Adams-Moulton</a></td>
|
|
|
430 <td>variable</td>
|
|
|
431 <td>variable</td>
|
|
|
432 </tr>
|
|
|
433 </table>
|
|
|
434 </p>
|
|
|
435 </div>
|
|
|
436 <div class="section"><h3><a name="a13.5_Derivatives"></a>13.5 Derivatives</h3>
|
|
|
437 <p>
|
|
|
438 If in addition to state y(t) the user needs to compute the sensitivity of the state to
|
|
|
439 the initial state or some parameter of the ODE, he will use the classes and interfaces
|
|
|
440 from the <a href="../apidocs/org/apache/commons/math/ode/jacobians/package-summary.html">org.apache.commons.ode.jacobians</a>
|
|
|
441 package instead of the top level ode package. These classes compute the jacobians
|
|
|
442 dy(t)/dy<sub>0</sub> and dy(t<sub>0</sub>)/dp where y<sub>0</sub> is the initial state
|
|
|
443 and p is some ODE parameter.
|
|
|
444 </p>
|
|
|
445 <p>
|
|
|
446 The classes and interfaces in this package mimic the behavior of the classes and
|
|
|
447 interfaces of the top level ode package, only adding parameters arrays for the jacobians.
|
|
|
448 The behavior of these classes is to create a compound state vector z containing both
|
|
|
449 the state y(t) and its derivatives dy(t)/dy<sub>0</sub> and dy(t<sub>0</sub>)/dp and
|
|
|
450 to set up an extended problem by adding the equations for the jacobians automatically.
|
|
|
451 These extended state and problems are then provided to a classical underlying integrator
|
|
|
452 chosen by user.
|
|
|
453 </p>
|
|
|
454 <p>
|
|
|
455 This behavior imply there will be a top level integrator knowing about state and jacobians
|
|
|
456 and a low level integrator knowing only about compound state (which may be big). If the user
|
|
|
457 wants to deal with the top level only, he will use the specialized step handler and event
|
|
|
458 handler classes registered at top level. He can also register classical step handlers and
|
|
|
459 event handlers, but in this case will see the big compound state. This state is guaranteed
|
|
|
460 to contain the original state in the first elements, followed by the jacobian with respect
|
|
|
461 to initial state (in row order), followed by the jacobian with respect to parameters (in
|
|
|
462 row order). If for example the original state dimension is 6 and there are 3 parameters,
|
|
|
463 the compound state will be a 60 elements array. The first 6 elements will be the original
|
|
|
464 state, the next 36 elements will be the jacobian with respect to initial state, and the
|
|
|
465 remaining 18 will be the jacobian with respect to parameters. Dealing with low level
|
|
|
466 step handlers and event handlers is cumbersome if one really needs the jacobians in these
|
|
|
467 methods, but it also prevents many data being copied back and forth between state and
|
|
|
468 jacobians on one side and compound state on the other side.
|
|
|
469 </p>
|
|
|
470 <p>
|
|
|
471 In order to compute dy(t)/dy<sub>0</sub> and dy(t<sub>0</sub>)/dp for any t, the algorithm
|
|
|
472 needs not only the ODE function f such that y'=f(t,y) but also its local jacobians
|
|
|
473 df(t, y, p)/dy and df(t, y, p)/dp.
|
|
|
474 </p>
|
|
|
475 <p>
|
|
|
476 If the function f is too complex, the user can simply rely on internal differentiation
|
|
|
477 using finite differences to compute these local jacobians. So rather than the <a href="../apidocs/org/apache/commons/math/ode/FirstOrderDifferentialEquations.html">FirstOrderDifferentialEquations</a>
|
|
|
478 interface he will implement the <a href="../apidocs/org/apache/commons/math/ode/jacobians/ParameterizedODE.html">ParameterizedODE</a>
|
|
|
479 interface. Considering again our example where only ? is considered a parameter, we get:
|
|
|
480 </p>
|
|
|
481 <div class="source"><pre>
|
|
|
482 public class BasicCircleODE implements ParameterizedODE {
|
|
|
483
|
|
|
484 private double[] c;
|
|
|
485 private double omega;
|
|
|
486
|
|
|
487 public BasicCircleODE(double[] c, double omega) {
|
|
|
488 this.c = c;
|
|
|
489 this.omega = omega;
|
|
|
490 }
|
|
|
491
|
|
|
492 public int getDimension() {
|
|
|
493 return 2;
|
|
|
494 }
|
|
|
495
|
|
|
496 public void computeDerivatives(double t, double[] y, double[] yDot) {
|
|
|
497 yDot[0] = omega * (c[1] - y[1]);
|
|
|
498 yDot[1] = omega * (y[0] - c[0]);
|
|
|
499 }
|
|
|
500
|
|
|
501 public int getParametersDimension() {
|
|
|
502 // we are only interested in the omega parameter
|
|
|
503 return 1;
|
|
|
504 }
|
|
|
505
|
|
|
506 public void setParameter(int i, double value) {
|
|
|
507 omega = value;
|
|
|
508 }
|
|
|
509
|
|
|
510 }
|
|
|
511 </pre>
|
|
|
512 </div>
|
|
|
513 <p>
|
|
|
514 This ODE is provided to the specialized integrator with two arrays specifying the
|
|
|
515 step sizes to use for finite differences (one array for derivation with respect to
|
|
|
516 state y, one array for derivation with respect to parameters p):
|
|
|
517 </p>
|
|
|
518 <div class="source"><pre>
|
|
|
519 double[] hY = new double[] { 0.001, 0.001 };
|
|
|
520 double[] hP = new double[] { 1.0e-6 };
|
|
|
521 FirstOrderIntegratorWithJacobians integrator = new FirstOrderIntegratorWithJacobians(dp853, ode, hY, hP);
|
|
|
522 integrator.integrate(t0, y0, dy0dp, t, y, dydy0, dydp);
|
|
|
523 </pre>
|
|
|
524 </div>
|
|
|
525 <p>
|
|
|
526 If the function f is simple, the user can simply provide the local jacobians
|
|
|
527 by himself. So rather than the <a href="../apidocs/org/apache/commons/math/ode/FirstOrderDifferentialEquations.html">FirstOrderDifferentialEquations</a>
|
|
|
528 interface he will implement the <a href="../apidocs/org/apache/commons/math/ode/jacobians/ODEWithJacobians.html">ODEWithJacobians</a>
|
|
|
529 interface. Considering again our example where only ? is considered a parameter, we get:
|
|
|
530 </p>
|
|
|
531 <div class="source"><pre>
|
|
|
532 public class EnhancedCircleODE implements ODEWithJacobians {
|
|
|
533
|
|
|
534 private double[] c;
|
|
|
535 private double omega;
|
|
|
536
|
|
|
537 public EnhancedCircleODE(double[] c, double omega) {
|
|
|
538 this.c = c;
|
|
|
539 this.omega = omega;
|
|
|
540 }
|
|
|
541
|
|
|
542 public int getDimension() {
|
|
|
543 return 2;
|
|
|
544 }
|
|
|
545
|
|
|
546 public void computeDerivatives(double t, double[] y, double[] yDot) {
|
|
|
547 yDot[0] = omega * (c[1] - y[1]);
|
|
|
548 yDot[1] = omega * (y[0] - c[0]);
|
|
|
549 }
|
|
|
550
|
|
|
551 public int getParametersDimension() {
|
|
|
552 // we are only interested in the omega parameter
|
|
|
553 return 1;
|
|
|
554 }
|
|
|
555
|
|
|
556 public void computeJacobians(double t, double[] y, double[] yDot, double[][] dFdY, double[][] dFdP) {
|
|
|
557
|
|
|
558 dFdY[0][0] = 0;
|
|
|
559 dFdY[0][1] = -omega;
|
|
|
560 dFdY[1][0] = omega;
|
|
|
561 dFdY[1][1] = 0;
|
|
|
562
|
|
|
563 dFdP[0][0] = 0;
|
|
|
564 dFdP[0][1] = omega;
|
|
|
565 dFdP[0][2] = c[1] - y[1];
|
|
|
566 dFdP[1][0] = -omega;
|
|
|
567 dFdP[1][1] = 0;
|
|
|
568 dFdP[1][2] = y[0] - c[0];
|
|
|
569
|
|
|
570 }
|
|
|
571
|
|
|
572 }
|
|
|
573 </pre>
|
|
|
574 </div>
|
|
|
575 <p>
|
|
|
576 This ODE is provided to the specialized integrator as is:
|
|
|
577 </p>
|
|
|
578 <div class="source"><pre>
|
|
|
579 FirstOrderIntegratorWithJacobians integrator = new FirstOrderIntegratorWithJacobians(dp853, ode);
|
|
|
580 integrator.integrate(t0, y0, dy0dp, t, y, dydy0, dydp);
|
|
|
581 </pre>
|
|
|
582 </div>
|
|
|
583 </div>
|
|
|
584 </div>
|
|
|
585
|
|
|
586 </div>
|
|
|
587 </div>
|
|
|
588 <div class="clear">
|
|
|
589 <hr/>
|
|
|
590 </div>
|
|
|
591 <div id="footer">
|
|
|
592 <div class="xright">©
|
|
|
593 2003-2010
|
|
|
594
|
|
|
595
|
|
|
596
|
|
|
597
|
|
|
598
|
|
|
599
|
|
|
600
|
|
|
601
|
|
|
602
|
|
|
603 </div>
|
|
|
604 <div class="clear">
|
|
|
605 <hr/>
|
|
|
606 </div>
|
|
|
607 </div>
|
|
|
608 </body>
|
|
|
609 </html>
|
