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 public 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 }