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.math.ode;
18
19 import java.util.ArrayList;
20 import java.util.List;
21
22 import org.apache.commons.math.exception.DimensionMismatchException;
23
24
25 /**
26 * This class represents a combined set of first order differential equations,
27 * with at least a primary set of equations expandable by some sets of secondary
28 * equations.
29 * <p>
30 * One typical use case is the computation of the Jacobian matrix for some ODE.
31 * In this case, the primary set of equations corresponds to the raw ODE, and we
32 * add to this set another bunch of secondary equations which represent the Jacobian
33 * matrix of the primary set.
34 * </p>
35 * <p>
36 * We want the integrator to use <em>only</em> the primary set to estimate the
37 * errors and hence the step sizes. It should <em>not</em> use the secondary
38 * equations in this computation. The {@link AbstractIntegrator integrator} will
39 * be able to know where the primary set ends and so where the secondary sets begin.
40 * </p>
41 *
42 * @see FirstOrderDifferentialEquations
43 * @see JacobianMatrices
44 *
45 * @version $Id: ExpandableStatefulODE.java 1178235 2011-10-02 19:43:17Z luc $
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 */
111 public void computeDerivatives(final double t, final double[] y, final double[] yDot) {
112
113 // compute derivatives of the primary equations
114 primaryMapper.extractEquationData(y, primaryState);
115 primary.computeDerivatives(t, primaryState, primaryStateDot);
116 primaryMapper.insertEquationData(primaryStateDot, yDot);
117
118 // Add contribution for secondary equations
119 for (final SecondaryComponent component : components) {
120 component.mapper.extractEquationData(y, component.state);
121 component.equation.computeDerivatives(t, primaryState, primaryStateDot,
122 component.state, component.stateDot);
123 component.mapper.insertEquationData(component.stateDot, yDot);
124 }
125
126 }
127
128 /** Add a set of secondary equations to be integrated along with the primary set.
129 * @param secondary secondary equations set
130 * @return index of the secondary equation in the expanded state
131 */
132 public int addSecondaryEquations(final SecondaryEquations secondary) {
133
134 final int firstIndex;
135 if (components.isEmpty()) {
136 // lazy creation of the components list
137 components = new ArrayList<ExpandableStatefulODE.SecondaryComponent>();
138 firstIndex = primary.getDimension();
139 } else {
140 final SecondaryComponent last = components.get(components.size() - 1);
141 firstIndex = last.mapper.getFirstIndex() + last.mapper.getDimension();
142 }
143
144 components.add(new SecondaryComponent(secondary, firstIndex));
145
146 return components.size() - 1;
147
148 }
149
150 /** Get an equations mapper for the primary equations set.
151 * @return mapper for the primary set
152 * @see #getSecondaryMappers()
153 */
154 public EquationsMapper getPrimaryMapper() {
155 return primaryMapper;
156 }
157
158 /** Get the equations mappers for the secondary equations sets.
159 * @return equations mappers for the secondary equations sets
160 * @see #getPrimaryMapper()
161 */
162 public EquationsMapper[] getSecondaryMappers() {
163 final EquationsMapper[] mappers = new EquationsMapper[components.size()];
164 for (int i = 0; i < mappers.length; ++i) {
165 mappers[i] = components.get(i).mapper;
166 }
167 return mappers;
168 }
169
170 /** Set current time.
171 * @param time current time
172 */
173 public void setTime(final double time) {
174 this.time = time;
175 }
176
177 /** Get current time.
178 * @return current time
179 */
180 public double getTime() {
181 return time;
182 }
183
184 /** Set primary part of the current state.
185 * @param primaryState primary part of the current state
186 * @throws DimensionMismatchException if the dimension of the array does not
187 * match the primary set
188 */
189 public void setPrimaryState(final double[] primaryState) throws DimensionMismatchException {
190
191 // safety checks
192 if (primaryState.length != this.primaryState.length) {
193 throw new DimensionMismatchException(primaryState.length, this.primaryState.length);
194 }
195
196 // set the data
197 System.arraycopy(primaryState, 0, this.primaryState, 0, primaryState.length);
198
199 }
200
201 /** Get primary part of the current state.
202 * @return primary part of the current state
203 */
204 public double[] getPrimaryState() {
205 return primaryState.clone();
206 }
207
208 /** Get primary part of the current state derivative.
209 * @return primary part of the current state derivative
210 */
211 public double[] getPrimaryStateDot() {
212 return primaryStateDot.clone();
213 }
214
215 /** Set secondary part of the current state.
216 * @param index index of the part to set as returned by {@link
217 * #addSecondaryEquations(SecondaryEquations)}
218 * @param secondaryState secondary part of the current state
219 * @throws DimensionMismatchException if the dimension of the partial state does not
220 * match the selected equations set dimension
221 */
222 public void setSecondaryState(final int index, final double[] secondaryState)
223 throws DimensionMismatchException {
224
225 // get either the secondary state
226 double[] localArray = components.get(index).state;
227
228 // safety checks
229 if (secondaryState.length != localArray.length) {
230 throw new DimensionMismatchException(secondaryState.length, localArray.length);
231 }
232
233 // set the data
234 System.arraycopy(secondaryState, 0, localArray, 0, secondaryState.length);
235
236 }
237
238 /** Get secondary part of the current state.
239 * @param index index of the part to set as returned by {@link
240 * #addSecondaryEquations(SecondaryEquations)}
241 * @return secondary part of the current state
242 */
243 public double[] getSecondaryState(final int index) {
244 return components.get(index).state.clone();
245 }
246
247 /** Get secondary part of the current state derivative.
248 * @param index index of the part to set as returned by {@link
249 * #addSecondaryEquations(SecondaryEquations)}
250 * @return secondary part of the current state derivative
251 */
252 public double[] getSecondaryStateDot(final int index) {
253 return components.get(index).stateDot.clone();
254 }
255
256 /** Set the complete current state.
257 * @param completeState complete current state to copy data from
258 * @throws DimensionMismatchException if the dimension of the complete state does not
259 * match the complete equations sets dimension
260 */
261 public void setCompleteState(final double[] completeState)
262 throws DimensionMismatchException {
263
264 // safety checks
265 if (completeState.length != getTotalDimension()) {
266 throw new DimensionMismatchException(completeState.length, getTotalDimension());
267 }
268
269 // set the data
270 primaryMapper.extractEquationData(completeState, primaryState);
271 for (final SecondaryComponent component : components) {
272 component.mapper.extractEquationData(completeState, component.state);
273 }
274
275 }
276
277 /** Get the complete current state.
278 * @return complete current state
279 * @throws DimensionMismatchException if the dimension of the complete state does not
280 * match the complete equations sets dimension
281 */
282 public double[] getCompleteState() {
283
284 // allocate complete array
285 double[] completeState = new double[getTotalDimension()];
286
287 // set the data
288 primaryMapper.insertEquationData(primaryState, completeState);
289 for (final SecondaryComponent component : components) {
290 component.mapper.insertEquationData(component.state, completeState);
291 }
292
293 return completeState;
294
295 }
296
297 /** Components of the compound stateful ODE. */
298 private static class SecondaryComponent {
299
300 /** Secondary differential equation. */
301 private final SecondaryEquations equation;
302
303 /** Mapper between local and complete arrays. */
304 private final EquationsMapper mapper;
305
306 /** State. */
307 private final double[] state;
308
309 /** State derivative. */
310 private final double[] stateDot;
311
312 /** Simple constructor.
313 * @param equation secondary differential equation
314 * @param firstIndex index to use for the first element in the complete arrays
315 */
316 public SecondaryComponent(final SecondaryEquations equation, final int firstIndex) {
317 final int n = equation.getDimension();
318 this.equation = equation;
319 mapper = new EquationsMapper(firstIndex, n);
320 state = new double[n];
321 stateDot = new double[n];
322 }
323
324 }
325
326 }