View Javadoc

1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  package org.apache.commons.math3.ode;
18  
19  import java.util.ArrayList;
20  import java.util.List;
21  
22  import org.apache.commons.math3.exception.DimensionMismatchException;
23  import org.apache.commons.math3.exception.MaxCountExceededException;
24  
25  
26  /**
27   * This class represents a combined set of first order differential equations,
28   * with at least a primary set of equations expandable by some sets of secondary
29   * equations.
30   * <p>
31   * One typical use case is the computation of the Jacobian matrix for some ODE.
32   * In this case, the primary set of equations corresponds to the raw ODE, and we
33   * add to this set another bunch of secondary equations which represent the Jacobian
34   * matrix of the primary set.
35   * </p>
36   * <p>
37   * We want the integrator to use <em>only</em> the primary set to estimate the
38   * errors and hence the step sizes. It should <em>not</em> use the secondary
39   * equations in this computation. The {@link AbstractIntegrator integrator} will
40   * be able to know where the primary set ends and so where the secondary sets begin.
41   * </p>
42   *
43   * @see FirstOrderDifferentialEquations
44   * @see JacobianMatrices
45   *
46   * @version $Id: ExpandableStatefulODE.java 1463680 2013-04-02 19:02:55Z luc $
47   * @since 3.0
48   */
49  
50  public class ExpandableStatefulODE {
51  
52      /** Primary differential equation. */
53      private final FirstOrderDifferentialEquations primary;
54  
55      /** Mapper for primary equation. */
56      private final EquationsMapper primaryMapper;
57  
58      /** Time. */
59      private double time;
60  
61      /** State. */
62      private final double[] primaryState;
63  
64      /** State derivative. */
65      private final double[] primaryStateDot;
66  
67      /** Components of the expandable ODE. */
68      private List<SecondaryComponent> components;
69  
70      /** Build an expandable set from its primary ODE set.
71       * @param primary the primary set of differential equations to be integrated.
72       */
73      public ExpandableStatefulODE(final FirstOrderDifferentialEquations primary) {
74          final int n          = primary.getDimension();
75          this.primary         = primary;
76          this.primaryMapper   = new EquationsMapper(0, n);
77          this.time            = Double.NaN;
78          this.primaryState    = new double[n];
79          this.primaryStateDot = new double[n];
80          this.components      = new ArrayList<ExpandableStatefulODE.SecondaryComponent>();
81      }
82  
83      /** Get the primary set of differential equations.
84       * @return primary set of differential equations
85       */
86      public FirstOrderDifferentialEquations getPrimary() {
87          return primary;
88      }
89  
90      /** Return the dimension of the complete set of equations.
91       * <p>
92       * The complete set of equations correspond to the primary set plus all secondary sets.
93       * </p>
94       * @return dimension of the complete set of equations
95       */
96      public int getTotalDimension() {
97          if (components.isEmpty()) {
98              // there are no secondary equations, the complete set is limited to the primary set
99              return primaryMapper.getDimension();
100         } else {
101             // there are secondary equations, the complete set ends after the last set
102             final EquationsMapper lastMapper = components.get(components.size() - 1).mapper;
103             return lastMapper.getFirstIndex() + lastMapper.getDimension();
104         }
105     }
106 
107     /** Get the current time derivative of the complete state vector.
108      * @param t current value of the independent <I>time</I> variable
109      * @param y array containing the current value of the complete state vector
110      * @param yDot placeholder array where to put the time derivative of the complete state vector
111      * @exception MaxCountExceededException if the number of functions evaluations is exceeded
112      * @exception DimensionMismatchException if arrays dimensions do not match equations settings
113      */
114     public void computeDerivatives(final double t, final double[] y, final double[] yDot)
115         throws MaxCountExceededException, DimensionMismatchException {
116 
117         // compute derivatives of the primary equations
118         primaryMapper.extractEquationData(y, primaryState);
119         primary.computeDerivatives(t, primaryState, primaryStateDot);
120 
121         // Add contribution for secondary equations
122         for (final SecondaryComponent component : components) {
123             component.mapper.extractEquationData(y, component.state);
124             component.equation.computeDerivatives(t, primaryState, primaryStateDot,
125                                                   component.state, component.stateDot);
126             component.mapper.insertEquationData(component.stateDot, yDot);
127         }
128 
129         primaryMapper.insertEquationData(primaryStateDot, yDot);
130 
131     }
132 
133     /** Add a set of secondary equations to be integrated along with the primary set.
134      * @param secondary secondary equations set
135      * @return index of the secondary equation in the expanded state
136      */
137     public int addSecondaryEquations(final SecondaryEquations secondary) {
138 
139         final int firstIndex;
140         if (components.isEmpty()) {
141             // lazy creation of the components list
142             components = new ArrayList<ExpandableStatefulODE.SecondaryComponent>();
143             firstIndex = primary.getDimension();
144         } else {
145             final SecondaryComponent last = components.get(components.size() - 1);
146             firstIndex = last.mapper.getFirstIndex() + last.mapper.getDimension();
147         }
148 
149         components.add(new SecondaryComponent(secondary, firstIndex));
150 
151         return components.size() - 1;
152 
153     }
154 
155     /** Get an equations mapper for the primary equations set.
156      * @return mapper for the primary set
157      * @see #getSecondaryMappers()
158      */
159     public EquationsMapper getPrimaryMapper() {
160         return primaryMapper;
161     }
162 
163     /** Get the equations mappers for the secondary equations sets.
164      * @return equations mappers for the secondary equations sets
165      * @see #getPrimaryMapper()
166      */
167     public EquationsMapper[] getSecondaryMappers() {
168         final EquationsMapper[] mappers = new EquationsMapper[components.size()];
169         for (int i = 0; i < mappers.length; ++i) {
170             mappers[i] = components.get(i).mapper;
171         }
172         return mappers;
173     }
174 
175     /** Set current time.
176      * @param time current time
177      */
178     public void setTime(final double time) {
179         this.time = time;
180     }
181 
182     /** Get current time.
183      * @return current time
184      */
185     public double getTime() {
186         return time;
187     }
188 
189     /** Set primary part of the current state.
190      * @param primaryState primary part of the current state
191      * @throws DimensionMismatchException if the dimension of the array does not
192      * match the primary set
193      */
194     public void setPrimaryState(final double[] primaryState) throws DimensionMismatchException {
195 
196         // safety checks
197         if (primaryState.length != this.primaryState.length) {
198             throw new DimensionMismatchException(primaryState.length, this.primaryState.length);
199         }
200 
201         // set the data
202         System.arraycopy(primaryState, 0, this.primaryState, 0, primaryState.length);
203 
204     }
205 
206     /** Get primary part of the current state.
207      * @return primary part of the current state
208      */
209     public double[] getPrimaryState() {
210         return primaryState.clone();
211     }
212 
213     /** Get primary part of the current state derivative.
214      * @return primary part of the current state derivative
215      */
216     public double[] getPrimaryStateDot() {
217         return primaryStateDot.clone();
218     }
219 
220     /** Set secondary part of the current state.
221      * @param index index of the part to set as returned by {@link
222      * #addSecondaryEquations(SecondaryEquations)}
223      * @param secondaryState secondary part of the current state
224      * @throws DimensionMismatchException if the dimension of the partial state does not
225      * match the selected equations set dimension
226      */
227     public void setSecondaryState(final int index, final double[] secondaryState)
228         throws DimensionMismatchException {
229 
230         // get either the secondary state
231         double[] localArray = components.get(index).state;
232 
233         // safety checks
234         if (secondaryState.length != localArray.length) {
235             throw new DimensionMismatchException(secondaryState.length, localArray.length);
236         }
237 
238         // set the data
239         System.arraycopy(secondaryState, 0, localArray, 0, secondaryState.length);
240 
241     }
242 
243     /** Get secondary part of the current state.
244      * @param index index of the part to set as returned by {@link
245      * #addSecondaryEquations(SecondaryEquations)}
246      * @return secondary part of the current state
247      */
248     public double[] getSecondaryState(final int index) {
249         return components.get(index).state.clone();
250     }
251 
252     /** Get secondary part of the current state derivative.
253      * @param index index of the part to set as returned by {@link
254      * #addSecondaryEquations(SecondaryEquations)}
255      * @return secondary part of the current state derivative
256      */
257     public double[] getSecondaryStateDot(final int index) {
258         return components.get(index).stateDot.clone();
259     }
260 
261     /** Set the complete current state.
262      * @param completeState complete current state to copy data from
263      * @throws DimensionMismatchException if the dimension of the complete state does not
264      * match the complete equations sets dimension
265      */
266     public void setCompleteState(final double[] completeState)
267         throws DimensionMismatchException {
268 
269         // safety checks
270         if (completeState.length != getTotalDimension()) {
271             throw new DimensionMismatchException(completeState.length, getTotalDimension());
272         }
273 
274         // set the data
275         primaryMapper.extractEquationData(completeState, primaryState);
276         for (final SecondaryComponent component : components) {
277             component.mapper.extractEquationData(completeState, component.state);
278         }
279 
280     }
281 
282     /** Get the complete current state.
283      * @return complete current state
284      * @throws DimensionMismatchException if the dimension of the complete state does not
285      * match the complete equations sets dimension
286      */
287     public double[] getCompleteState() throws DimensionMismatchException {
288 
289         // allocate complete array
290         double[] completeState = new double[getTotalDimension()];
291 
292         // set the data
293         primaryMapper.insertEquationData(primaryState, completeState);
294         for (final SecondaryComponent component : components) {
295             component.mapper.insertEquationData(component.state, completeState);
296         }
297 
298         return completeState;
299 
300     }
301 
302     /** Components of the compound stateful ODE. */
303     private static class SecondaryComponent {
304 
305         /** Secondary differential equation. */
306         private final SecondaryEquations equation;
307 
308         /** Mapper between local and complete arrays. */
309         private final EquationsMapper mapper;
310 
311         /** State. */
312         private final double[] state;
313 
314         /** State derivative. */
315         private final double[] stateDot;
316 
317         /** Simple constructor.
318          * @param equation secondary differential equation
319          * @param firstIndex index to use for the first element in the complete arrays
320          */
321         public SecondaryComponent(final SecondaryEquations equation, final int firstIndex) {
322             final int n   = equation.getDimension();
323             this.equation = equation;
324             mapper        = new EquationsMapper(firstIndex, n);
325             state         = new double[n];
326             stateDot      = new double[n];
327         }
328 
329     }
330 
331 }