FixedStepRungeKuttaFieldIntegrator.java

  1. /*
  2.  * Licensed to the Apache Software Foundation (ASF) 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 ASF 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.  */

  18. /*
  19.  * This is not the original file distributed by the Apache Software Foundation
  20.  * It has been modified by the Hipparchus project
  21.  */

  23. package org.hipparchus.ode.nonstiff;


  26. import org.hipparchus.CalculusFieldElement;
  27. import org.hipparchus.Field;
  28. import org.hipparchus.exception.MathIllegalArgumentException;
  29. import org.hipparchus.exception.MathIllegalStateException;
  30. import org.hipparchus.ode.AbstractFieldIntegrator;
  31. import org.hipparchus.ode.FieldEquationsMapper;
  32. import org.hipparchus.ode.FieldExpandableODE;
  33. import org.hipparchus.ode.FieldODEState;
  34. import org.hipparchus.ode.FieldODEStateAndDerivative;
  35. import org.hipparchus.ode.nonstiff.interpolators.RungeKuttaFieldStateInterpolator;
  36. import org.hipparchus.util.MathArrays;

  38. /**
  39.  * This class implements the common part of all fixed step Runge-Kutta
  40.  * integrators for Ordinary Differential Equations.
  41.  *
  42.  * <p>These methods are explicit Runge-Kutta methods, their Butcher
  43.  * arrays are as follows :</p>
  44.  * <pre>
  45.  *    0  |
  46.  *   c2  | a21
  47.  *   c3  | a31  a32
  48.  *   ... |        ...
  49.  *   cs  | as1  as2  ...  ass-1
  50.  *       |--------------------------
  51.  *       |  b1   b2  ...   bs-1  bs
  52.  * </pre>
  53.  *
  54.  * @see EulerFieldIntegrator
  55.  * @see ClassicalRungeKuttaFieldIntegrator
  56.  * @see GillFieldIntegrator
  57.  * @see MidpointFieldIntegrator
  58.  * @param <T> the type of the field elements
  59.  */

  61. public abstract class FixedStepRungeKuttaFieldIntegrator<T extends CalculusFieldElement<T>>
  62.     extends AbstractFieldIntegrator<T> implements FieldExplicitRungeKuttaIntegrator<T> {

  64.     /** Time steps from Butcher array (without the first zero). */
  65.     private final T[] c;

  67.     /** Internal weights from Butcher array (without the first empty row). */
  68.     private final T[][] a;

  70.     /** External weights for the high order method from Butcher array. */
  71.     private final T[] b;

  73.     /** Time steps from Butcher array (without the first zero). */
  74.     private double[] realC = new double[0];

  76.     /** Internal weights from Butcher array (without the first empty row). Real version, optional. */
  77.     private double[][] realA = new double[0][];

  79.     /** External weights for the high order method from Butcher array. Real version, optional. */
  80.     private double[] realB = new double[0];

  82.     /** Integration step. */
  83.     private final T step;

  85.     /** Flag setting whether coefficients in Butcher array are interpreted as Field or real numbers. */
  86.     private boolean usingFieldCoefficients;

  88.     /** Simple constructor.
  89.      * Build a Runge-Kutta integrator with the given
  90.      * step. The default step handler does nothing.
  91.      * @param field field to which the time and state vector elements belong
  92.      * @param name name of the method
  93.      * @param step integration step
  94.      */
  95.     protected FixedStepRungeKuttaFieldIntegrator(final Field<T> field, final String name, final T step) {
  96.         super(field, name);
  97.         this.c    = getC();
  98.         this.a    = getA();
  99.         this.b    = getB();
  100.         this.step = step.abs();
  101.         this.usingFieldCoefficients = false;
  102.     }

  104.     /** Getter for the default, positive step-size assigned at constructor level.
  105.      * @return step
  106.      */
  107.     public T getDefaultStep() {
  108.         return this.step;
  109.     }

  111.     /**
  112.      * Setter for the flag between real or Field coefficients in the Butcher array.
  113.      *
  114.      * @param usingFieldCoefficients new value for flag
  115.      */
  116.     public void setUsingFieldCoefficients(boolean usingFieldCoefficients) {
  117.         this.usingFieldCoefficients = usingFieldCoefficients;
  118.     }

  120.     /** {@inheritDoc} */
  121.     @Override
  122.     public boolean isUsingFieldCoefficients() {
  123.         return usingFieldCoefficients;
  124.     }

  126.     /** {@inheritDoc} */
  127.     @Override
  128.     public int getNumberOfStages() {
  129.         return b.length;
  130.     }

  132.     /** Create an interpolator.
  133.      * @param forward integration direction indicator
  134.      * @param yDotK slopes at the intermediate points
  135.      * @param globalPreviousState start of the global step
  136.      * @param globalCurrentState end of the global step
  137.      * @param mapper equations mapper for the all equations
  138.      * @return external weights for the high order method from Butcher array
  139.      */
  140.     protected abstract RungeKuttaFieldStateInterpolator<T> createInterpolator(boolean forward, T[][] yDotK,
  141.                                                                               FieldODEStateAndDerivative<T> globalPreviousState,
  142.                                                                               FieldODEStateAndDerivative<T> globalCurrentState,
  143.                                                                               FieldEquationsMapper<T> mapper);

  145.     /** {@inheritDoc} */
  146.     @Override
  147.     protected FieldODEStateAndDerivative<T> initIntegration(FieldExpandableODE<T> eqn, FieldODEState<T> s0, T t) {
  148.         if (!isUsingFieldCoefficients()) {
  149.             realA = getRealA();
  150.             realB = getRealB();
  151.             realC = getRealC();
  152.         }
  153.         return super.initIntegration(eqn, s0, t);
  154.     }

  156.     /** {@inheritDoc} */
  157.     @Override
  158.     public FieldODEStateAndDerivative<T> integrate(final FieldExpandableODE<T> equations,
  159.                                                    final FieldODEState<T> initialState, final T finalTime)
  160.         throws MathIllegalArgumentException, MathIllegalStateException {

  162.         sanityChecks(initialState, finalTime);
  163.         setStepStart(initIntegration(equations, initialState, finalTime));
  164.         final boolean forward = finalTime.subtract(initialState.getTime()).getReal() > 0;

  166.         // create some internal working arrays
  167.         final int   stages = getNumberOfStages();
  168.         final T[][] yDotK  = MathArrays.buildArray(getField(), stages, -1);
  169.         MathArrays.buildArray(getField(), equations.getMapper().getTotalDimension());

  171.         // set up integration control objects
  172.         if (forward) {
  173.             if (getStepStart().getTime().add(step).subtract(finalTime).getReal() >= 0) {
  174.                 setStepSize(finalTime.subtract(getStepStart().getTime()));
  175.             } else {
  176.                 setStepSize(step);
  177.             }
  178.         } else {
  179.             if (getStepStart().getTime().subtract(step).subtract(finalTime).getReal() <= 0) {
  180.                 setStepSize(finalTime.subtract(getStepStart().getTime()));
  181.             } else {
  182.                 setStepSize(step.negate());
  183.             }
  184.         }

  186.         // main integration loop
  187.         setIsLastStep(false);
  188.         do {

  190.             // first stage
  191.             final T[] y = getStepStart().getCompleteState();
  192.             yDotK[0]    = getStepStart().getCompleteDerivative();

  194.             // next stages
  195.             final T[] yTmp;
  196.             if (isUsingFieldCoefficients()) {
  197.                 FieldExplicitRungeKuttaIntegrator.applyInternalButcherWeights(getEquations(), getStepStart().getTime(),
  198.                         y, getStepSize(), a, c, yDotK);
  199.                 yTmp = FieldExplicitRungeKuttaIntegrator.applyExternalButcherWeights(y, yDotK, getStepSize(), b);
  200.             } else {
  201.                 FieldExplicitRungeKuttaIntegrator.applyInternalButcherWeights(getEquations(), getStepStart().getTime(),
  202.                         y, getStepSize(), realA, realC, yDotK);
  203.                 yTmp = FieldExplicitRungeKuttaIntegrator.applyExternalButcherWeights(y, yDotK, getStepSize(), realB);
  204.             }

  206.             incrementEvaluations(stages - 1);

  208.             final T stepEnd   = getStepStart().getTime().add(getStepSize());
  209.             final T[] yDotTmp = computeDerivatives(stepEnd, yTmp);
  210.             final FieldODEStateAndDerivative<T> stateTmp = equations.getMapper().mapStateAndDerivative(stepEnd, yTmp, yDotTmp);

  212.             // discrete events handling
  213.             setStepStart(acceptStep(createInterpolator(forward, yDotK, getStepStart(), stateTmp, equations.getMapper()),
  214.                                     finalTime));

  216.             if (!isLastStep()) {

  218.                 // stepsize control for next step
  219.                 final T  nextT      = getStepStart().getTime().add(getStepSize());
  220.                 final boolean nextIsLast = forward ?
  221.                                            (nextT.subtract(finalTime).getReal() >= 0) :
  222.                                            (nextT.subtract(finalTime).getReal() <= 0);
  223.                 if (nextIsLast) {
  224.                     setStepSize(finalTime.subtract(getStepStart().getTime()));
  225.                 }
  226.             }

  228.         } while (!isLastStep());

  230.         final FieldODEStateAndDerivative<T> finalState = getStepStart();
  231.         setStepStart(null);
  232.         setStepSize(null);
  233.         return finalState;

  235.     }

  237. }
Created with JaCoCo 0.8.12.202403310830