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   * @since 3.0
47   */
48  
49  public class ExpandableStatefulODE {
50  
51      /** Primary differential equation. */
52      private final FirstOrderDifferentialEquations primary;
53  
54      /** Mapper for primary equation. */
55      private final EquationsMapper primaryMapper;
56  
57      /** Time. */
58      private double time;
59  
60      /** State. */
61      private final double[] primaryState;
62  
63      /** State derivative. */
64      private final double[] primaryStateDot;
65  
66      /** Components of the expandable ODE. */
67      private List<SecondaryComponent> components;
68  
69      /** Build an expandable set from its primary ODE set.
70       * @param primary the primary set of differential equations to be integrated.
71       */
72      public ExpandableStatefulODE(final FirstOrderDifferentialEquations primary) {
73          final int n          = primary.getDimension();
74          this.primary         = primary;
75          this.primaryMapper   = new EquationsMapper(0, n);
76          this.time            = Double.NaN;
77          this.primaryState    = new double[n];
78          this.primaryStateDot = new double[n];
79          this.components      = new ArrayList<ExpandableStatefulODE.SecondaryComponent>();
80      }
81  
82      /** Get the primary set of differential equations.
83       * @return primary set of differential equations
84       */
85      public FirstOrderDifferentialEquations getPrimary() {
86          return primary;
87      }
88  
89      /** Return the dimension of the complete set of equations.
90       * <p>
91       * The complete set of equations correspond to the primary set plus all secondary sets.
92       * </p>
93       * @return dimension of the complete set of equations
94       */
95      public int getTotalDimension() {
96          if (components.isEmpty()) {
97              // there are no secondary equations, the complete set is limited to the primary set
98              return primaryMapper.getDimension();
99          } else {
100             // there are secondary equations, the complete set ends after the last set
101             final EquationsMapper lastMapper = components.get(components.size() - 1).mapper;
102             return lastMapper.getFirstIndex() + lastMapper.getDimension();
103         }
104     }
105 
106     /** Get the current time derivative of the complete state vector.
107      * @param t current value of the independent <I>time</I> variable
108      * @param y array containing the current value of the complete state vector
109      * @param yDot placeholder array where to put the time derivative of the complete state vector
110      * @exception MaxCountExceededException if the number of functions evaluations is exceeded
111      * @exception DimensionMismatchException if arrays dimensions do not match equations settings
112      */
113     public void computeDerivatives(final double t, final double[] y, final double[] yDot)
114         throws MaxCountExceededException, DimensionMismatchException {
115 
116         // compute derivatives of the primary equations
117         primaryMapper.extractEquationData(y, primaryState);
118         primary.computeDerivatives(t, primaryState, primaryStateDot);
119 
120         // Add contribution for secondary equations
121         for (final SecondaryComponent component : components) {
122             component.mapper.extractEquationData(y, component.state);
123             component.equation.computeDerivatives(t, primaryState, primaryStateDot,
124                                                   component.state, component.stateDot);
125             component.mapper.insertEquationData(component.stateDot, yDot);
126         }
127 
128         primaryMapper.insertEquationData(primaryStateDot, yDot);
129 
130     }
131 
132     /** Add a set of secondary equations to be integrated along with the primary set.
133      * @param secondary secondary equations set
134      * @return index of the secondary equation in the expanded state
135      */
136     public int addSecondaryEquations(final SecondaryEquations secondary) {
137 
138         final int firstIndex;
139         if (components.isEmpty()) {
140             // lazy creation of the components list
141             components = new ArrayList<ExpandableStatefulODE.SecondaryComponent>();
142             firstIndex = primary.getDimension();
143         } else {
144             final SecondaryComponent last = components.get(components.size() - 1);
145             firstIndex = last.mapper.getFirstIndex() + last.mapper.getDimension();
146         }
147 
148         components.add(new SecondaryComponent(secondary, firstIndex));
149 
150         return components.size() - 1;
151 
152     }
153 
154     /** Get an equations mapper for the primary equations set.
155      * @return mapper for the primary set
156      * @see #getSecondaryMappers()
157      */
158     public EquationsMapper getPrimaryMapper() {
159         return primaryMapper;
160     }
161 
162     /** Get the equations mappers for the secondary equations sets.
163      * @return equations mappers for the secondary equations sets
164      * @see #getPrimaryMapper()
165      */
166     public EquationsMapper[] getSecondaryMappers() {
167         final EquationsMapper[] mappers = new EquationsMapper[components.size()];
168         for (int i = 0; i < mappers.length; ++i) {
169             mappers[i] = components.get(i).mapper;
170         }
171         return mappers;
172     }
173 
174     /** Set current time.
175      * @param time current time
176      */
177     public void setTime(final double time) {
178         this.time = time;
179     }
180 
181     /** Get current time.
182      * @return current time
183      */
184     public double getTime() {
185         return time;
186     }
187 
188     /** Set primary part of the current state.
189      * @param primaryState primary part of the current state
190      * @throws DimensionMismatchException if the dimension of the array does not
191      * match the primary set
192      */
193     public void setPrimaryState(final double[] primaryState) throws DimensionMismatchException {
194 
195         // safety checks
196         if (primaryState.length != this.primaryState.length) {
197             throw new DimensionMismatchException(primaryState.length, this.primaryState.length);
198         }
199 
200         // set the data
201         System.arraycopy(primaryState, 0, this.primaryState, 0, primaryState.length);
202 
203     }
204 
205     /** Get primary part of the current state.
206      * @return primary part of the current state
207      */
208     public double[] getPrimaryState() {
209         return primaryState.clone();
210     }
211 
212     /** Get primary part of the current state derivative.
213      * @return primary part of the current state derivative
214      */
215     public double[] getPrimaryStateDot() {
216         return primaryStateDot.clone();
217     }
218 
219     /** Set secondary part of the current state.
220      * @param index index of the part to set as returned by {@link
221      * #addSecondaryEquations(SecondaryEquations)}
222      * @param secondaryState secondary part of the current state
223      * @throws DimensionMismatchException if the dimension of the partial state does not
224      * match the selected equations set dimension
225      */
226     public void setSecondaryState(final int index, final double[] secondaryState)
227         throws DimensionMismatchException {
228 
229         // get either the secondary state
230         double[] localArray = components.get(index).state;
231 
232         // safety checks
233         if (secondaryState.length != localArray.length) {
234             throw new DimensionMismatchException(secondaryState.length, localArray.length);
235         }
236 
237         // set the data
238         System.arraycopy(secondaryState, 0, localArray, 0, secondaryState.length);
239 
240     }
241 
242     /** Get secondary part of the current state.
243      * @param index index of the part to set as returned by {@link
244      * #addSecondaryEquations(SecondaryEquations)}
245      * @return secondary part of the current state
246      */
247     public double[] getSecondaryState(final int index) {
248         return components.get(index).state.clone();
249     }
250 
251     /** Get secondary part of the current state derivative.
252      * @param index index of the part to set as returned by {@link
253      * #addSecondaryEquations(SecondaryEquations)}
254      * @return secondary part of the current state derivative
255      */
256     public double[] getSecondaryStateDot(final int index) {
257         return components.get(index).stateDot.clone();
258     }
259 
260     /** Set the complete current state.
261      * @param completeState complete current state to copy data from
262      * @throws DimensionMismatchException if the dimension of the complete state does not
263      * match the complete equations sets dimension
264      */
265     public void setCompleteState(final double[] completeState)
266         throws DimensionMismatchException {
267 
268         // safety checks
269         if (completeState.length != getTotalDimension()) {
270             throw new DimensionMismatchException(completeState.length, getTotalDimension());
271         }
272 
273         // set the data
274         primaryMapper.extractEquationData(completeState, primaryState);
275         for (final SecondaryComponent component : components) {
276             component.mapper.extractEquationData(completeState, component.state);
277         }
278 
279     }
280 
281     /** Get the complete current state.
282      * @return complete current state
283      * @throws DimensionMismatchException if the dimension of the complete state does not
284      * match the complete equations sets dimension
285      */
286     public double[] getCompleteState() throws DimensionMismatchException {
287 
288         // allocate complete array
289         double[] completeState = new double[getTotalDimension()];
290 
291         // set the data
292         primaryMapper.insertEquationData(primaryState, completeState);
293         for (final SecondaryComponent component : components) {
294             component.mapper.insertEquationData(component.state, completeState);
295         }
296 
297         return completeState;
298 
299     }
300 
301     /** Components of the compound stateful ODE. */
302     private static class SecondaryComponent {
303 
304         /** Secondary differential equation. */
305         private final SecondaryEquations equation;
306 
307         /** Mapper between local and complete arrays. */
308         private final EquationsMapper mapper;
309 
310         /** State. */
311         private final double[] state;
312 
313         /** State derivative. */
314         private final double[] stateDot;
315 
316         /** Simple constructor.
317          * @param equation secondary differential equation
318          * @param firstIndex index to use for the first element in the complete arrays
319          */
320         SecondaryComponent(final SecondaryEquations equation, final int firstIndex) {
321             final int n   = equation.getDimension();
322             this.equation = equation;
323             mapper        = new EquationsMapper(firstIndex, n);
324             state         = new double[n];
325             stateDot      = new double[n];
326         }
327 
328     }
329 
330 }