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 */
017package org.apache.commons.math3.ode;
018
019import java.util.ArrayList;
020import java.util.List;
021
022import org.apache.commons.math3.exception.DimensionMismatchException;
023import 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 1463680 2013-04-02 19:02:55Z luc $
047 * @since 3.0
048 */
049
050public 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
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}