    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    }