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    }