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 * @since 3.0
047 */
048
049public class ExpandableStatefulODE {
050
051    /** Primary differential equation. */
052    private final FirstOrderDifferentialEquations primary;
053
054    /** Mapper for primary equation. */
055    private final EquationsMapper primaryMapper;
056
057    /** Time. */
058    private double time;
059
060    /** State. */
061    private final double[] primaryState;
062
063    /** State derivative. */
064    private final double[] primaryStateDot;
065
066    /** Components of the expandable ODE. */
067    private List<SecondaryComponent> components;
068
069    /** Build an expandable set from its primary ODE set.
070     * @param primary the primary set of differential equations to be integrated.
071     */
072    public ExpandableStatefulODE(final FirstOrderDifferentialEquations primary) {
073        final int n          = primary.getDimension();
074        this.primary         = primary;
075        this.primaryMapper   = new EquationsMapper(0, n);
076        this.time            = Double.NaN;
077        this.primaryState    = new double[n];
078        this.primaryStateDot = new double[n];
079        this.components      = new ArrayList<ExpandableStatefulODE.SecondaryComponent>();
080    }
081
082    /** Get the primary set of differential equations.
083     * @return primary set of differential equations
084     */
085    public FirstOrderDifferentialEquations getPrimary() {
086        return primary;
087    }
088
089    /** Return the dimension of the complete set of equations.
090     * <p>
091     * The complete set of equations correspond to the primary set plus all secondary sets.
092     * </p>
093     * @return dimension of the complete set of equations
094     */
095    public int getTotalDimension() {
096        if (components.isEmpty()) {
097            // there are no secondary equations, the complete set is limited to the primary set
098            return primaryMapper.getDimension();
099        } 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}