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 }