VariationalEquation.java

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

  19. import java.lang.reflect.Array;

  21. import org.hipparchus.exception.LocalizedCoreFormats;
  22. import org.hipparchus.exception.MathIllegalArgumentException;
  23. import org.hipparchus.exception.MathIllegalStateException;

  25. /**
  26.  * This class defines a set of {@link SecondaryODE secondary equations} to
  27.  * compute the global Jacobian matrices with respect to the initial state
  28.  * vector and, if any, to some parameters of the primary ODE set.
  29.  * <p>
  30.  * The primary set of ODE for which Jaobian matrices are requested may be:
  31.  * </p>
  32.  * <ul>
  33.  * <li>a full-fledged {@link ODEJacobiansProvider} that computes by itself
  34.  * both the ODE and its local partial derivatives,</li>
  35.  * <li>a simple {@link OrdinaryDifferentialEquation} which must therefore
  36.  * be completed with a finite differences configuration to compute local
  37.  * partial derivatives (so-called internal differentiation).</li>
  38.  * </ul>
  39.  * <p>
  40.  * As the variational equation automatically inserts {@link
  41.  * ExpandableODE#addSecondaryEquations(SecondaryODE) secondary differential
  42.  * equations}, in the {@link ExpandableODE expandable ODE}, data for
  43.  * initial state must also be inserted before integration and matrices
  44.  * result must be extracted after integration. This implies a precise
  45.  * scheduling of the calls to the various methods of this class. The
  46.  * proper scheduling is the following one:
  47.  * </p>
  48.  * <pre>
  49.  *   // set up equations
  50.  *   ODEJacobiansProvider jode       = new MyODE(...);
  51.  *   ExpandableODE        expandable = new Expandable(jode);
  52.  *   VariationalEquation  ve         = new VariationalEquation(expandable, jode);
  53.  *
  54.  *   // set up initial state
  55.  *   ODEState initWithoutDerivatives = new ODEState(t0, y0);
  56.  *   ve.setInitialMainStateJacobian(dYdY0); // only needed if the default identity matrix is not suitable
  57.  *   ve.setInitialParameterJacobian(name, dYdP); // only needed if the default zero matrix is not suitable
  58.  *   ODEState initWithDerivatives = ve.setUpInitialState(initWithoutDerivatives);
  59.  *
  60.  *   // perform integration on the expanded equations with the expanded initial state
  61.  *   ODEStateAndDerivative finalState = integrator.integrate(expandable, initWithDerivatives, finalT);
  62.  *
  63.  *   // extract Jacobian matrices
  64.  *   dYdY0 = ve.extractMainSetJacobian(finalState);
  65.  *   dYdP  = ve.extractParameterJacobian(finalState, name);
  66.  * </pre>
  67.  * <p>
  68.  * The most important part is to not forget to call {@link #setUpInitialState(ODEState)} to add
  69.  * the secondary state with the initial matrices to the {@link ODEState} used in the
  70.  * {@link ODEIntegrator#integrate(ExpandableODE, ODEState, double) integrate} method.
  71.  * Forgetting to do this and passing only a {@link ODEState} without the secondary state
  72.  * set up will trigger an error as the state vector will not have the correct dimension.
  73.  * </p>
  74.  *
  75.  * @see ExpandableODE
  76.  * @see ODEJacobiansProvider
  77.  * @see OrdinaryDifferentialEquation
  78.  * @see NamedParameterJacobianProvider
  79.  * @see ParametersController
  80.  *
  81.  */
  82. public class VariationalEquation {

  84.     /** ODE with Jacobian computation skill. */
  85.     private final ODEJacobiansProvider jode;

  87.     /** Expandable first order differential equation. */
  88.     private final ExpandableODE expandable;

  90.     /** Index of the instance in the expandable set. */
  91.     private final int index;

  93.     /** State and parameters Jacobian matrices in a row. */
  94.     private double[] matricesData;

  96.     /** Build variational equation using finite differences for local
  97.      * partial derivatives.
  98.      * @param expandable expandable set into which variational equations should be registered
  99.      * @param ode base ordinary differential equation for which Jacobians
  100.      * matrices are requested
  101.      * @param hY step used for finite difference computation with respect to state vector
  102.      * @param controller controller to change parameters
  103.      * @param paramsAndSteps parameters and steps to compute the Jacobians df/dp
  104.      * @exception MismatchedEquations if the primary set of the expandable set does
  105.      * not match the {@code ode}
  106.      */
  107.     public VariationalEquation(final ExpandableODE expandable,
  108.                                final OrdinaryDifferentialEquation ode, final double[] hY,
  109.                                final ParametersController controller,
  110.                                final ParameterConfiguration ... paramsAndSteps)
  111.         throws MismatchedEquations {
  112.         this(expandable, new ParameterJacobianWrapper(ode, hY, controller, paramsAndSteps));
  113.     }

  115.     /** Build variational equation using analytical local partial derivatives.
  116.      * <p>
  117.      * Parameters must belong to the supported ones given by {@link
  118.      * Parameterizable#getParametersNames()}, so the primary set of differential
  119.      * equations must be {@link Parameterizable}.
  120.      * </p>
  121.      * <p>Note that each selection clears the previous selected parameters.</p>
  122.      *
  123.      * @param expandable expandable set into which variational equations should be registered
  124.      * @param jode the primary first order differential equations set to extend
  125.      * @exception MismatchedEquations if the primary set of the expandable set does
  126.      * not match the {@code ode}
  127.      */
  128.     public VariationalEquation(final ExpandableODE expandable,
  129.                                final ODEJacobiansProvider jode)
  130.         throws MismatchedEquations {

  132.         // safety checks
  133.         final OrdinaryDifferentialEquation ode;
  134.         if (jode instanceof ParameterJacobianWrapper) {
  135.             ode = ((ParameterJacobianWrapper) jode).getODE();
  136.         } else {
  137.             ode = jode;
  138.         }
  139.         if (expandable.getPrimary() != ode) {
  140.             throw new MismatchedEquations();
  141.         }

  143.         this.jode       = jode;
  144.         this.expandable = expandable;
  145.         this.index      = expandable.addSecondaryEquations(new JacobiansSecondaryODE());

  147.         // set the default initial state Jacobian to the identity
  148.         // and the default initial parameters Jacobian to the null matrix
  149.         matricesData = new double[(jode.getDimension() + jode.getParametersNames().size()) * jode.getDimension()];
  150.         for (int i = 0; i < jode.getDimension(); ++i) {
  151.             matricesData[i * (jode.getDimension() + 1)] = 1.0;
  152.         }

  154.     }

  156.     /** Set the initial value of the Jacobian matrix with respect to state.
  157.      * <p>
  158.      * If this method is not called, the initial value of the Jacobian
  159.      * matrix with respect to state is set to identity.
  160.      * </p>
  161.      * <p>
  162.      * This method must be called <em>before {@link #setUpInitialState(ODEState)}</em>
  163.      * </p>
  164.      * @param dYdY0 initial Jacobian matrix w.r.t. state
  165.      * @exception MathIllegalArgumentException if matrix dimensions are incorrect
  166.      */
  167.     public void setInitialMainStateJacobian(final double[][] dYdY0)
  168.         throws MathIllegalArgumentException {

  170.         // Check dimensions
  171.         checkDimension(jode.getDimension(), dYdY0);
  172.         checkDimension(jode.getDimension(), dYdY0[0]);

  174.         // store the matrix in row major order as a single dimension array
  175.         int i = 0;
  176.         for (final double[] row : dYdY0) {
  177.             System.arraycopy(row, 0, matricesData, i, jode.getDimension());
  178.             i += jode.getDimension();
  179.         }

  181.     }

  183.     /** Set the initial value of a column of the Jacobian matrix with respect to one parameter.
  184.      * <p>
  185.      * If this method is not called for some parameter, the initial value of
  186.      * the column of the Jacobian matrix with respect to this parameter is set to zero.
  187.      * </p>
  188.      * <p>
  189.      * This method must be called <em>before {@link #setUpInitialState(ODEState)}</em>
  190.      * </p>
  191.      * @param pName parameter name
  192.      * @param dYdP initial Jacobian column vector with respect to the parameter
  193.      * @exception MathIllegalArgumentException if a parameter is not supported
  194.      * @throws MathIllegalArgumentException if the column vector does not match state dimension
  195.      */
  196.     public void setInitialParameterJacobian(final String pName, final double[] dYdP)
  197.         throws MathIllegalArgumentException {

  199.         // Check dimensions
  200.         checkDimension(jode.getDimension(), dYdP);

  202.         // store the column in a global single dimension array
  203.         int i = jode.getDimension() * jode.getDimension();
  204.         for (final String knownParameter : jode.getParametersNames()) {
  205.             if (pName.equals(knownParameter)) {
  206.                 System.arraycopy(dYdP, 0, matricesData, i, jode.getDimension());
  207.                 return;
  208.             }
  209.             i += jode.getDimension();
  210.         }

  212.         throw new MathIllegalArgumentException(LocalizedODEFormats.UNKNOWN_PARAMETER, pName);

  214.     }

  216.     /** Set up initial state.
  217.      * <p>
  218.      * This method inserts the initial Jacobian matrices data into
  219.      * an {@link ODEState ODE state} by overriding the additional
  220.      * state components corresponding to the instance. It must be
  221.      * called prior to integrate the equations.
  222.      * </p>
  223.      * <p>
  224.      * This method must be called <em>after</em>
  225.      * {@link #setInitialMainStateJacobian(double[][])} and
  226.      * {@link #setInitialParameterJacobian(String, double[])}.
  227.      * </p>
  228.      * @param initialState initial state, without the initial Jacobians
  229.      * matrices
  230.      * @return a new instance of initial state, with the initial Jacobians
  231.      * matrices properly initialized
  232.      */
  233.     public ODEState setUpInitialState(final ODEState initialState) { // NOPMD - PMD false positive

  235.         // insert the matrices data into secondary states
  236.         final double[][] secondary = new double[expandable.getMapper().getNumberOfEquations() - 1][];
  237.         for (int i = 0; i < initialState.getNumberOfSecondaryStates(); ++i) {
  238.             if (i + 1 != index) {
  239.                 secondary[i] = initialState.getSecondaryState(i + 1);
  240.             }
  241.         }
  242.         secondary[index - 1] = matricesData;

  244.         // create an updated initial state
  245.         return new ODEState(initialState.getTime(), initialState.getPrimaryState(), secondary);

  247.     }

  249.     /** Extract the Jacobian matrix with respect to state.
  250.      * @param state state from which to extract Jacobian matrix
  251.      * @return Jacobian matrix dY/dY0 with respect to state.
  252.      */
  253.     public double[][] extractMainSetJacobian(final ODEState state) {

  255.         // get current state for this set of equations from the expandable fode
  256.         final double[] p = state.getSecondaryState(index);

  258.         final double[][] dYdY0 = new double[jode.getDimension()][jode.getDimension()];
  259.         int j = 0;
  260.         for (int i = 0; i < jode.getDimension(); i++) {
  261.             System.arraycopy(p, j, dYdY0[i], 0, jode.getDimension());
  262.             j += jode.getDimension();
  263.         }

  265.         return dYdY0;

  267.     }

  269.     /** Extract the Jacobian matrix with respect to one parameter.
  270.      * @param state state from which to extract Jacobian matrix
  271.      * @param pName name of the parameter for the computed Jacobian matrix
  272.      * @return Jacobian matrix dY/dP with respect to the named parameter
  273.      */
  274.     public double[] extractParameterJacobian(final ODEState state, final String pName) {

  276.         // get current state for this set of equations from the expandable fode
  277.         final double[] p = state.getSecondaryState(index);

  279.         final double[] dYdP = new double[jode.getDimension()];
  280.         int i = jode.getDimension() * jode.getDimension();
  281.         for (final String knownParameter : jode.getParametersNames()) {
  282.             if (pName.equals(knownParameter)) {
  283.                 System.arraycopy(p, i, dYdP, 0, jode.getDimension());
  284.                 break;
  285.             }
  286.             i += jode.getDimension();
  287.         }

  289.         return dYdP;

  291.     }

  293.     /** Check array dimensions.
  294.      * @param expected expected dimension
  295.      * @param array (may be null if expected is 0)
  296.      * @throws MathIllegalArgumentException if the array dimension does not match the expected one
  297.      */
  298.     private void checkDimension(final int expected, final Object array)
  299.         throws MathIllegalArgumentException {
  300.         int arrayDimension = (array == null) ? 0 : Array.getLength(array);
  301.         if (arrayDimension != expected) {
  302.             throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
  303.                                                    arrayDimension, expected);
  304.         }
  305.     }

  307.     /** Local implementation of secondary equations. */
  308.     private class JacobiansSecondaryODE implements SecondaryODE {

  310.         /** {@inheritDoc} */
  311.         @Override
  312.         public int getDimension() {
  313.             return jode.getDimension() * (jode.getDimension() + jode.getParametersNames().size());
  314.         }

  316.         /** {@inheritDoc} */
  317.         @Override
  318.         public double[] computeDerivatives(final double t, final double[] y, final double[] yDot,
  319.                                            final double[] z)
  320.             throws MathIllegalArgumentException, MathIllegalStateException {

  322.             final double[] zDot = new double[z.length];

  324.             // variational equations:
  325.             // from d[dy/dt]/dy0 and d[dy/dt]/dp to d[dy/dy0]/dt and d[dy/dp]/dt

  327.             // compute Jacobian matrix with respect to primary state
  328.             double[][] dFdY = jode.computeMainStateJacobian(t, y, yDot);

  330.             // Dispatch Jacobian matrix in the compound secondary state vector
  331.             for (int i = 0; i < jode.getDimension(); ++i) {
  332.                 final double[] dFdYi = dFdY[i];
  333.                 for (int j = 0; j < jode.getDimension(); ++j) {
  334.                     double s = 0;
  335.                     final int startIndex = j;
  336.                     int zIndex = startIndex;
  337.                     for (int l = 0; l < jode.getDimension(); ++l) {
  338.                         s += dFdYi[l] * z[zIndex];
  339.                         zIndex += jode.getDimension();
  340.                     }
  341.                     zDot[startIndex + i * jode.getDimension()] = s;
  342.                 }
  343.             }

  345.             // compute Jacobian matrices with respect to parameters
  346.             int startIndex = jode.getDimension() * jode.getDimension();
  347.             for (final String name : jode.getParametersNames()) {
  348.                 final double[] dFdP = jode.computeParameterJacobian(t, y, yDot, name);
  349.                 for (int i = 0; i < jode.getDimension(); ++i) {
  350.                     final double[] dFdYi = dFdY[i];
  351.                     int zIndex = startIndex;
  352.                     double s = dFdP[i];
  353.                     for (int l = 0; l < jode.getDimension(); ++l) {
  354.                         s += dFdYi[l] * z[zIndex];
  355.                         zIndex++;
  356.                     }
  357.                     zDot[startIndex + i] = s;
  358.                 }
  359.                 startIndex += jode.getDimension();
  360.             }

  362.             return zDot;

  364.         }
  365.     }

  367.     /**
  368.      * Special exception for equations mismatch.
  369.      */
  370.     public static class MismatchedEquations extends MathIllegalArgumentException {

  372.         /** Serializable UID. */
  373.         private static final long serialVersionUID = 20120902L;

  375.         /** Simple constructor. */
  376.         public MismatchedEquations() {
  377.             super(LocalizedODEFormats.UNMATCHED_ODE_IN_EXPANDED_SET);
  378.         }

  380.     }

  382. }
Created with JaCoCo 0.8.11.202310140853