001 /*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements. See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License. You may obtain a copy of the License at
008 *
009 * http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017 package org.apache.commons.math3.ode;
018
019 import java.util.ArrayList;
020 import java.util.List;
021
022 import org.apache.commons.math3.exception.DimensionMismatchException;
023 import org.apache.commons.math3.exception.MaxCountExceededException;
024
025
026 /**
027 * This class represents a combined set of first order differential equations,
028 * with at least a primary set of equations expandable by some sets of secondary
029 * equations.
030 * <p>
031 * One typical use case is the computation of the Jacobian matrix for some ODE.
032 * In this case, the primary set of equations corresponds to the raw ODE, and we
033 * add to this set another bunch of secondary equations which represent the Jacobian
034 * matrix of the primary set.
035 * </p>
036 * <p>
037 * We want the integrator to use <em>only</em> the primary set to estimate the
038 * errors and hence the step sizes. It should <em>not</em> use the secondary
039 * equations in this computation. The {@link AbstractIntegrator integrator} will
040 * be able to know where the primary set ends and so where the secondary sets begin.
041 * </p>
042 *
043 * @see FirstOrderDifferentialEquations
044 * @see JacobianMatrices
045 *
046 * @version $Id: ExpandableStatefulODE.java 1416643 2012-12-03 19:37:14Z tn $
047 * @since 3.0
048 */
049
050 public class ExpandableStatefulODE {
051
052 /** Primary differential equation. */
053 private final FirstOrderDifferentialEquations primary;
054
055 /** Mapper for primary equation. */
056 private final EquationsMapper primaryMapper;
057
058 /** Time. */
059 private double time;
060
061 /** State. */
062 private final double[] primaryState;
063
064 /** State derivative. */
065 private final double[] primaryStateDot;
066
067 /** Components of the expandable ODE. */
068 private List<SecondaryComponent> components;
069
070 /** Build an expandable set from its primary ODE set.
071 * @param primary the primary set of differential equations to be integrated.
072 */
073 public ExpandableStatefulODE(final FirstOrderDifferentialEquations primary) {
074 final int n = primary.getDimension();
075 this.primary = primary;
076 this.primaryMapper = new EquationsMapper(0, n);
077 this.time = Double.NaN;
078 this.primaryState = new double[n];
079 this.primaryStateDot = new double[n];
080 this.components = new ArrayList<ExpandableStatefulODE.SecondaryComponent>();
081 }
082
083 /** Get the primary set of differential equations.
084 * @return primary set of differential equations
085 */
086 public FirstOrderDifferentialEquations getPrimary() {
087 return primary;
088 }
089
090 /** Return the dimension of the complete set of equations.
091 * <p>
092 * The complete set of equations correspond to the primary set plus all secondary sets.
093 * </p>
094 * @return dimension of the complete set of equations
095 */
096 public int getTotalDimension() {
097 if (components.isEmpty()) {
098 // there are no secondary equations, the complete set is limited to the primary set
099 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 primaryMapper.insertEquationData(primaryStateDot, yDot);
121
122 // Add contribution for secondary equations
123 for (final SecondaryComponent component : components) {
124 component.mapper.extractEquationData(y, component.state);
125 component.equation.computeDerivatives(t, primaryState, primaryStateDot,
126 component.state, component.stateDot);
127 component.mapper.insertEquationData(component.stateDot, yDot);
128 }
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 }