ExpandableODE.java

/*
 * Licensed to the Hipparchus project under one or more
 * contributor license agreements.  See the NOTICE file distributed with
 * this work for additional information regarding copyright ownership.
 * The Hipparchus project licenses this file to You under the Apache License, Version 2.0
 * (the "License"); you may not use this file except in compliance with
 * the License.  You may obtain a copy of the License at
 *
 *      https://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
package org.hipparchus.ode;

import java.util.ArrayList;
import java.util.List;

import org.hipparchus.exception.MathIllegalArgumentException;
import org.hipparchus.exception.MathIllegalStateException;


/**
 * This class represents a combined set of first order differential equations,
 * with at least a primary set of equations expandable by some sets of secondary
 * equations.
 * <p>
 * One typical use case is the computation of the Jacobian matrix for some ODE.
 * In this case, the primary set of equations corresponds to the raw ODE, and we
 * add to this set another bunch of secondary equations which represent the Jacobian
 * matrix of the primary set.
 * </p>
 * <p>
 * We want the integrator to use <em>only</em> the primary set to estimate the
 * errors and hence the step sizes. It should <em>not</em> use the secondary
 * equations in this computation. The {@link AbstractIntegrator integrator} will
 * be able to know where the primary set ends and so where the secondary sets begin.
 * </p>
 *
 * @see OrdinaryDifferentialEquation
 * @see VariationalEquation
 *
 */

public class ExpandableODE {

    /** Primary differential equation. */
    private final OrdinaryDifferentialEquation primary;

    /** Components of the expandable ODE. */
    private List<SecondaryODE> components;

    /** Mapper for all equations. */
    private EquationsMapper mapper;

    /** Build an expandable set from its primary ODE set.
     * @param primary the primary set of differential equations to be integrated.
     */
    public ExpandableODE(final OrdinaryDifferentialEquation primary) {
        this.primary    = primary;
        this.components = new ArrayList<>();
        this.mapper     = new EquationsMapper(null, primary.getDimension());
    }

    /** Get the primary set of differential equations to be integrated.
     * @return primary set of differential equations to be integrated
     */
    public OrdinaryDifferentialEquation getPrimary() {
        return primary;
    }

    /** Get the mapper for the set of equations.
     * @return mapper for the set of equations
     */
    public EquationsMapper getMapper() {
        return mapper;
    }

    /** Add a set of secondary equations to be integrated along with the primary set.
     * @param secondary secondary equations set
     * @return index of the secondary equation in the expanded state, to be used
     * as the parameter to {@link FieldODEState#getSecondaryState(int)} and
     * {@link FieldODEStateAndDerivative#getSecondaryDerivative(int)} (beware index
     * 0 corresponds to primary state, secondary states start at 1)
     */
    public int addSecondaryEquations(final SecondaryODE secondary) {

        components.add(secondary);
        mapper = new EquationsMapper(mapper, secondary.getDimension());

        return components.size();

    }

    /** Initialize equations at the start of an ODE integration.
     * @param s0 state at integration start
     * @param finalTime target time for the integration
     * @exception MathIllegalStateException if the number of functions evaluations is exceeded
     * @exception MathIllegalArgumentException if arrays dimensions do not match equations settings
     */
    public void init(final ODEState s0, final double finalTime) {

        final double t0 = s0.getTime();

        // initialize primary equations
        final double[] primary0 = s0.getPrimaryState();
        primary.init(t0, primary0, finalTime);

        // initialize secondary equations
        for (int index = 1; index < mapper.getNumberOfEquations(); ++index) {
            final double[] secondary0 = s0.getSecondaryState(index);
            components.get(index - 1).init(t0, primary0, secondary0, finalTime);
        }

    }

    /** Get the current time derivative of the complete state vector.
     * @param t current value of the independent <I>time</I> variable
     * @param y array containing the current value of the complete state vector
     * @return time derivative of the complete state vector
     * @exception MathIllegalStateException if the number of functions evaluations is exceeded
     * @exception MathIllegalArgumentException if arrays dimensions do not match equations settings
     */
    public double[] computeDerivatives(final double t, final double[] y)
        throws MathIllegalArgumentException, MathIllegalStateException {

        final double[] yDot = new double[mapper.getTotalDimension()];

        // compute derivatives of the primary equations
        final double[] primaryState    = mapper.extractEquationData(0, y);
        final double[] primaryStateDot = primary.computeDerivatives(t, primaryState);

        // Add contribution for secondary equations
        for (int index = 1; index < mapper.getNumberOfEquations(); ++index) {
            final double[] componentState    = mapper.extractEquationData(index, y);
            final double[] componentStateDot = components.get(index - 1).computeDerivatives(t, primaryState, primaryStateDot,
                                                                                            componentState);
            mapper.insertEquationData(index, componentStateDot, yDot);
        }

        // we retrieve the primaryStateDot array after the secondary equations have
        // been computed in case they change the main state derivatives; this happens
        // for example in optimal control when the secondary equations handle co-state,
        // which changes control, and the control changes the primary state
        mapper.insertEquationData(0, primaryStateDot, yDot);

        return yDot;

    }

}