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.math4.legacy.ode;
018
019import java.util.ArrayList;
020import java.util.List;
021
022import org.apache.commons.math4.legacy.core.RealFieldElement;
023import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
024import org.apache.commons.math4.legacy.exception.MaxCountExceededException;
025import org.apache.commons.math4.legacy.core.MathArrays;
026
027
028/**
029 * This class represents a combined set of first order differential equations,
030 * with at least a primary set of equations expandable by some sets of secondary
031 * equations.
032 * <p>
033 * One typical use case is the computation of the Jacobian matrix for some ODE.
034 * In this case, the primary set of equations corresponds to the raw ODE, and we
035 * add to this set another bunch of secondary equations which represent the Jacobian
036 * matrix of the primary set.
037 * </p>
038 * <p>
039 * We want the integrator to use <em>only</em> the primary set to estimate the
040 * errors and hence the step sizes. It should <em>not</em> use the secondary
041 * equations in this computation. The {@link FirstOrderFieldIntegrator integrator} will
042 * be able to know where the primary set ends and so where the secondary sets begin.
043 * </p>
044 *
045 * @see FirstOrderFieldDifferentialEquations
046 * @see FieldSecondaryEquations
047 *
048 * @param <T> the type of the field elements
049 * @since 3.6
050 */
051
052public class FieldExpandableODE<T extends RealFieldElement<T>> {
053
054    /** Primary differential equation. */
055    private final FirstOrderFieldDifferentialEquations<T> primary;
056
057    /** Components of the expandable ODE. */
058    private List<FieldSecondaryEquations<T>> components;
059
060    /** Mapper for all equations. */
061    private FieldEquationsMapper<T> mapper;
062
063    /** Build an expandable set from its primary ODE set.
064     * @param primary the primary set of differential equations to be integrated.
065     */
066    public FieldExpandableODE(final FirstOrderFieldDifferentialEquations<T> primary) {
067        this.primary    = primary;
068        this.components = new ArrayList<>();
069        this.mapper     = new FieldEquationsMapper<>(null, primary.getDimension());
070    }
071
072    /** Get the mapper for the set of equations.
073     * @return mapper for the set of equations
074     */
075    public FieldEquationsMapper<T> getMapper() {
076        return mapper;
077    }
078
079    /** Add a set of secondary equations to be integrated along with the primary set.
080     * @param secondary secondary equations set
081     * @return index of the secondary equation in the expanded state, to be used
082     * as the parameter to {@link FieldODEState#getSecondaryState(int)} and
083     * {@link FieldODEStateAndDerivative#getSecondaryDerivative(int)} (beware index
084     * 0 corresponds to main state, additional states start at 1)
085     */
086    public int addSecondaryEquations(final FieldSecondaryEquations<T> secondary) {
087
088        components.add(secondary);
089        mapper = new FieldEquationsMapper<>(mapper, secondary.getDimension());
090
091        return components.size();
092    }
093
094    /** Initialize equations at the start of an ODE integration.
095     * @param t0 value of the independent <I>time</I> variable at integration start
096     * @param y0 array containing the value of the state vector at integration start
097     * @param finalTime target time for the integration
098     * @exception MaxCountExceededException if the number of functions evaluations is exceeded
099     * @exception DimensionMismatchException if arrays dimensions do not match equations settings
100     */
101    public void init(final T t0, final T[] y0, final T finalTime) {
102
103        // initialize primary equations
104        int index = 0;
105        final T[] primary0 = mapper.extractEquationData(index, y0);
106        primary.init(t0, primary0, finalTime);
107
108        // initialize secondary equations
109        while (++index < mapper.getNumberOfEquations()) {
110            final T[] secondary0 = mapper.extractEquationData(index, y0);
111            components.get(index - 1).init(t0, primary0, secondary0, finalTime);
112        }
113    }
114
115    /** Get the current time derivative of the complete state vector.
116     * @param t current value of the independent <I>time</I> variable
117     * @param y array containing the current value of the complete state vector
118     * @return time derivative of the complete state vector
119     * @exception MaxCountExceededException if the number of functions evaluations is exceeded
120     * @exception DimensionMismatchException if arrays dimensions do not match equations settings
121     */
122    public T[] computeDerivatives(final T t, final T[] y)
123        throws MaxCountExceededException, DimensionMismatchException {
124
125        final T[] yDot = MathArrays.buildArray(t.getField(), mapper.getTotalDimension());
126
127        // compute derivatives of the primary equations
128        int index = 0;
129        final T[] primaryState    = mapper.extractEquationData(index, y);
130        final T[] primaryStateDot = primary.computeDerivatives(t, primaryState);
131        mapper.insertEquationData(index, primaryStateDot, yDot);
132
133        // Add contribution for secondary equations
134        while (++index < mapper.getNumberOfEquations()) {
135            final T[] componentState    = mapper.extractEquationData(index, y);
136            final T[] componentStateDot = components.get(index - 1).computeDerivatives(t, primaryState, primaryStateDot,
137                                                                                       componentState);
138            mapper.insertEquationData(index, componentStateDot, yDot);
139        }
140
141        return yDot;
142    }
143}