commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/ode/nonstiff/DormandPrince54FieldStepInterpolator.java - commons-math - Git at Google

 /*
  * Licensed to the Apache Software Foundation (ASF) under one or more
  * contributor license agreements.  See the NOTICE file distributed with
  * this work for additional information regarding copyright ownership.
  * The ASF 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
  *
  *      http://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.apache.commons.math4.legacy.ode.nonstiff;

 import org.apache.commons.math4.legacy.core.Field;
 import org.apache.commons.math4.legacy.core.RealFieldElement;
 import org.apache.commons.math4.legacy.ode.FieldEquationsMapper;
 import org.apache.commons.math4.legacy.ode.FieldODEStateAndDerivative;

 /**
  * This class represents an interpolator over the last step during an
  * ODE integration for the 5(4) Dormand-Prince integrator.
  *
  * @see DormandPrince54Integrator
  *
  * @param <T> the type of the field elements
  * @since 3.6
  */

 class DormandPrince54FieldStepInterpolator<T extends RealFieldElement<T>>
       extends RungeKuttaFieldStepInterpolator<T> {

     /** Last row of the Butcher-array internal weights, element 0. */
     private final T a70;

     // element 1 is zero, so it is neither stored nor used

     /** Last row of the Butcher-array internal weights, element 2. */
     private final T a72;

     /** Last row of the Butcher-array internal weights, element 3. */
     private final T a73;

     /** Last row of the Butcher-array internal weights, element 4. */
     private final T a74;

     /** Last row of the Butcher-array internal weights, element 5. */
     private final T a75;

     /** Shampine (1986) Dense output, element 0. */
     private final T d0;

     // element 1 is zero, so it is neither stored nor used

     /** Shampine (1986) Dense output, element 2. */
     private final T d2;

     /** Shampine (1986) Dense output, element 3. */
     private final T d3;

     /** Shampine (1986) Dense output, element 4. */
     private final T d4;

     /** Shampine (1986) Dense output, element 5. */
     private final T d5;

     /** Shampine (1986) Dense output, element 6. */
     private final T d6;

     /** Simple constructor.
      * @param field field to which the time and state vector elements belong
      * @param forward integration direction indicator
      * @param yDotK slopes at the intermediate points
      * @param globalPreviousState start of the global step
      * @param globalCurrentState end of the global step
      * @param softPreviousState start of the restricted step
      * @param softCurrentState end of the restricted step
      * @param mapper equations mapper for the all equations
      */
     DormandPrince54FieldStepInterpolator(final Field<T> field, final boolean forward,
                                          final T[][] yDotK,
                                          final FieldODEStateAndDerivative<T> globalPreviousState,
                                          final FieldODEStateAndDerivative<T> globalCurrentState,
                                          final FieldODEStateAndDerivative<T> softPreviousState,
                                          final FieldODEStateAndDerivative<T> softCurrentState,
                                          final FieldEquationsMapper<T> mapper) {
         super(field, forward, yDotK,
               globalPreviousState, globalCurrentState, softPreviousState, softCurrentState,
               mapper);
         final T one = field.getOne();
         a70 = one.multiply(   35.0).divide( 384.0);
         a72 = one.multiply(  500.0).divide(1113.0);
         a73 = one.multiply(  125.0).divide( 192.0);
         a74 = one.multiply(-2187.0).divide(6784.0);
         a75 = one.multiply(   11.0).divide(  84.0);
         d0  = one.multiply(-12715105075.0).divide( 11282082432.0);
         d2  = one.multiply( 87487479700.0).divide( 32700410799.0);
         d3  = one.multiply(-10690763975.0).divide(  1880347072.0);
         d4  = one.multiply(701980252875.0).divide(199316789632.0);
         d5  = one.multiply( -1453857185.0).divide(   822651844.0);
         d6  = one.multiply(    69997945.0).divide(    29380423.0);
     }

     /** {@inheritDoc} */
     @Override
     protected DormandPrince54FieldStepInterpolator<T> create(final Field<T> newField, final boolean newForward, final T[][] newYDotK,
                                                                  final FieldODEStateAndDerivative<T> newGlobalPreviousState,
                                                                  final FieldODEStateAndDerivative<T> newGlobalCurrentState,
                                                                  final FieldODEStateAndDerivative<T> newSoftPreviousState,
                                                                  final FieldODEStateAndDerivative<T> newSoftCurrentState,
                                                                  final FieldEquationsMapper<T> newMapper) {
         return new DormandPrince54FieldStepInterpolator<>(newField, newForward, newYDotK,
                                                            newGlobalPreviousState, newGlobalCurrentState,
                                                            newSoftPreviousState, newSoftCurrentState,
                                                            newMapper);
     }
     /** {@inheritDoc} */
     @SuppressWarnings("unchecked")
     @Override
     protected FieldODEStateAndDerivative<T> computeInterpolatedStateAndDerivatives(final FieldEquationsMapper<T> mapper,
                                                                                    final T time, final T theta,
                                                                                    final T thetaH, final T oneMinusThetaH) {

         // interpolate
         final T one      = time.getField().getOne();
         final T eta      = one.subtract(theta);
         final T twoTheta = theta.multiply(2);
         final T dot2     = one.subtract(twoTheta);
         final T dot3     = theta.multiply(theta.multiply(-3).add(2));
         final T dot4     = twoTheta.multiply(theta.multiply(twoTheta.subtract(3)).add(1));
         final T[] interpolatedState;
         final T[] interpolatedDerivatives;
         if (getGlobalPreviousState() != null && theta.getReal() <= 0.5) {
             final T f1        = thetaH;
             final T f2        = f1.multiply(eta);
             final T f3        = f2.multiply(theta);
             final T f4        = f3.multiply(eta);
             final T coeff0    = f1.multiply(a70).
                                 subtract(f2.multiply(a70.subtract(1))).
                                 add(f3.multiply(a70.multiply(2).subtract(1))).
                                 add(f4.multiply(d0));
             final T coeff1    = time.getField().getZero();
             final T coeff2    = f1.multiply(a72).
                                 subtract(f2.multiply(a72)).
                                 add(f3.multiply(a72.multiply(2))).
                                 add(f4.multiply(d2));
             final T coeff3    = f1.multiply(a73).
                                 subtract(f2.multiply(a73)).
                                 add(f3.multiply(a73.multiply(2))).
                                 add(f4.multiply(d3));
             final T coeff4    = f1.multiply(a74).
                                 subtract(f2.multiply(a74)).
                                 add(f3.multiply(a74.multiply(2))).
                                 add(f4.multiply(d4));
             final T coeff5    = f1.multiply(a75).
                                 subtract(f2.multiply(a75)).
                                 add(f3.multiply(a75.multiply(2))).
                                 add(f4.multiply(d5));
             final T coeff6    = f4.multiply(d6).subtract(f3);
             final T coeffDot0 = a70.
                                 subtract(dot2.multiply(a70.subtract(1))).
                                 add(dot3.multiply(a70.multiply(2).subtract(1))).
                                 add(dot4.multiply(d0));
             final T coeffDot1 = time.getField().getZero();
             final T coeffDot2 = a72.
                                 subtract(dot2.multiply(a72)).
                                 add(dot3.multiply(a72.multiply(2))).
                                 add(dot4.multiply(d2));
             final T coeffDot3 = a73.
                                 subtract(dot2.multiply(a73)).
                                 add(dot3.multiply(a73.multiply(2))).
                                 add(dot4.multiply(d3));
             final T coeffDot4 = a74.
                                 subtract(dot2.multiply(a74)).
                                 add(dot3.multiply(a74.multiply(2))).
                                 add(dot4.multiply(d4));
             final T coeffDot5 = a75.
                                 subtract(dot2.multiply(a75)).
                                 add(dot3.multiply(a75.multiply(2))).
                                 add(dot4.multiply(d5));
             final T coeffDot6 = dot4.multiply(d6).subtract(dot3);
             interpolatedState       = previousStateLinearCombination(coeff0, coeff1, coeff2, coeff3,
                                                                      coeff4, coeff5, coeff6);
             interpolatedDerivatives = derivativeLinearCombination(coeffDot0, coeffDot1, coeffDot2, coeffDot3,
                                                                   coeffDot4, coeffDot5, coeffDot6);
         } else {
             final T f1        = oneMinusThetaH.negate();
             final T f2        = oneMinusThetaH.multiply(theta);
             final T f3        = f2.multiply(theta);
             final T f4        = f3.multiply(eta);
             final T coeff0    = f1.multiply(a70).
                                 subtract(f2.multiply(a70.subtract(1))).
                                 add(f3.multiply(a70.multiply(2).subtract(1))).
                                 add(f4.multiply(d0));
             final T coeff1    = time.getField().getZero();
             final T coeff2    = f1.multiply(a72).
                                 subtract(f2.multiply(a72)).
                                 add(f3.multiply(a72.multiply(2))).
                                 add(f4.multiply(d2));
             final T coeff3    = f1.multiply(a73).
                                 subtract(f2.multiply(a73)).
                                 add(f3.multiply(a73.multiply(2))).
                                 add(f4.multiply(d3));
             final T coeff4    = f1.multiply(a74).
                                 subtract(f2.multiply(a74)).
                                 add(f3.multiply(a74.multiply(2))).
                                 add(f4.multiply(d4));
             final T coeff5    = f1.multiply(a75).
                                 subtract(f2.multiply(a75)).
                                 add(f3.multiply(a75.multiply(2))).
                                 add(f4.multiply(d5));
             final T coeff6    = f4.multiply(d6).subtract(f3);
             final T coeffDot0 = a70.
                                 subtract(dot2.multiply(a70.subtract(1))).
                                 add(dot3.multiply(a70.multiply(2).subtract(1))).
                                 add(dot4.multiply(d0));
             final T coeffDot1 = time.getField().getZero();
             final T coeffDot2 = a72.
                                 subtract(dot2.multiply(a72)).
                                 add(dot3.multiply(a72.multiply(2))).
                                 add(dot4.multiply(d2));
             final T coeffDot3 = a73.
                                 subtract(dot2.multiply(a73)).
                                 add(dot3.multiply(a73.multiply(2))).
                                 add(dot4.multiply(d3));
             final T coeffDot4 = a74.
                                 subtract(dot2.multiply(a74)).
                                 add(dot3.multiply(a74.multiply(2))).
                                 add(dot4.multiply(d4));
             final T coeffDot5 = a75.
                                 subtract(dot2.multiply(a75)).
                                 add(dot3.multiply(a75.multiply(2))).
                                 add(dot4.multiply(d5));
             final T coeffDot6 = dot4.multiply(d6).subtract(dot3);
             interpolatedState       = currentStateLinearCombination(coeff0, coeff1, coeff2, coeff3,
                                                                     coeff4, coeff5, coeff6);
             interpolatedDerivatives = derivativeLinearCombination(coeffDot0, coeffDot1, coeffDot2, coeffDot3,
                                                                   coeffDot4, coeffDot5, coeffDot6);
         }
         return new FieldODEStateAndDerivative<>(time, interpolatedState, interpolatedDerivatives);

     }

 }
	/*
	* Licensed to the Apache Software Foundation (ASF) under one or more
	* contributor license agreements. See the NOTICE file distributed with
	* this work for additional information regarding copyright ownership.
	* The ASF 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
	*
	* http://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.apache.commons.math4.legacy.ode.nonstiff;

	import org.apache.commons.math4.legacy.core.Field;
	import org.apache.commons.math4.legacy.core.RealFieldElement;
	import org.apache.commons.math4.legacy.ode.FieldEquationsMapper;
	import org.apache.commons.math4.legacy.ode.FieldODEStateAndDerivative;

	/**
	* This class represents an interpolator over the last step during an
	* ODE integration for the 5(4) Dormand-Prince integrator.
	*
	* @see DormandPrince54Integrator
	*
	* @param <T> the type of the field elements
	* @since 3.6
	*/

	class DormandPrince54FieldStepInterpolator<T extends RealFieldElement<T>>
	extends RungeKuttaFieldStepInterpolator<T> {

	/** Last row of the Butcher-array internal weights, element 0. */
	private final T a70;

	// element 1 is zero, so it is neither stored nor used

	/** Last row of the Butcher-array internal weights, element 2. */
	private final T a72;

	/** Last row of the Butcher-array internal weights, element 3. */
	private final T a73;

	/** Last row of the Butcher-array internal weights, element 4. */
	private final T a74;

	/** Last row of the Butcher-array internal weights, element 5. */
	private final T a75;

	/** Shampine (1986) Dense output, element 0. */
	private final T d0;

	// element 1 is zero, so it is neither stored nor used

	/** Shampine (1986) Dense output, element 2. */
	private final T d2;

	/** Shampine (1986) Dense output, element 3. */
	private final T d3;

	/** Shampine (1986) Dense output, element 4. */
	private final T d4;

	/** Shampine (1986) Dense output, element 5. */
	private final T d5;

	/** Shampine (1986) Dense output, element 6. */
	private final T d6;

	/** Simple constructor.
	* @param field field to which the time and state vector elements belong
	* @param forward integration direction indicator
	* @param yDotK slopes at the intermediate points
	* @param globalPreviousState start of the global step
	* @param globalCurrentState end of the global step
	* @param softPreviousState start of the restricted step
	* @param softCurrentState end of the restricted step
	* @param mapper equations mapper for the all equations
	*/
	DormandPrince54FieldStepInterpolator(final Field<T> field, final boolean forward,
	final T[][] yDotK,
	final FieldODEStateAndDerivative<T> globalPreviousState,
	final FieldODEStateAndDerivative<T> globalCurrentState,
	final FieldODEStateAndDerivative<T> softPreviousState,
	final FieldODEStateAndDerivative<T> softCurrentState,
	final FieldEquationsMapper<T> mapper) {
	super(field, forward, yDotK,
	globalPreviousState, globalCurrentState, softPreviousState, softCurrentState,
	mapper);
	final T one = field.getOne();
	a70 = one.multiply( 35.0).divide( 384.0);
	a72 = one.multiply( 500.0).divide(1113.0);
	a73 = one.multiply( 125.0).divide( 192.0);
	a74 = one.multiply(-2187.0).divide(6784.0);
	a75 = one.multiply( 11.0).divide( 84.0);
	d0 = one.multiply(-12715105075.0).divide( 11282082432.0);
	d2 = one.multiply( 87487479700.0).divide( 32700410799.0);
	d3 = one.multiply(-10690763975.0).divide( 1880347072.0);
	d4 = one.multiply(701980252875.0).divide(199316789632.0);
	d5 = one.multiply( -1453857185.0).divide( 822651844.0);
	d6 = one.multiply( 69997945.0).divide( 29380423.0);
	}

	/** {@inheritDoc} */
	@Override
	protected DormandPrince54FieldStepInterpolator<T> create(final Field<T> newField, final boolean newForward, final T[][] newYDotK,
	final FieldODEStateAndDerivative<T> newGlobalPreviousState,
	final FieldODEStateAndDerivative<T> newGlobalCurrentState,
	final FieldODEStateAndDerivative<T> newSoftPreviousState,
	final FieldODEStateAndDerivative<T> newSoftCurrentState,
	final FieldEquationsMapper<T> newMapper) {
	return new DormandPrince54FieldStepInterpolator<>(newField, newForward, newYDotK,
	newGlobalPreviousState, newGlobalCurrentState,
	newSoftPreviousState, newSoftCurrentState,
	newMapper);
	}
	/** {@inheritDoc} */
	@SuppressWarnings("unchecked")
	@Override
	protected FieldODEStateAndDerivative<T> computeInterpolatedStateAndDerivatives(final FieldEquationsMapper<T> mapper,
	final T time, final T theta,
	final T thetaH, final T oneMinusThetaH) {

	// interpolate
	final T one = time.getField().getOne();
	final T eta = one.subtract(theta);
	final T twoTheta = theta.multiply(2);
	final T dot2 = one.subtract(twoTheta);
	final T dot3 = theta.multiply(theta.multiply(-3).add(2));
	final T dot4 = twoTheta.multiply(theta.multiply(twoTheta.subtract(3)).add(1));
	final T[] interpolatedState;
	final T[] interpolatedDerivatives;
	if (getGlobalPreviousState() != null && theta.getReal() <= 0.5) {
	final T f1 = thetaH;
	final T f2 = f1.multiply(eta);
	final T f3 = f2.multiply(theta);
	final T f4 = f3.multiply(eta);
	final T coeff0 = f1.multiply(a70).
	subtract(f2.multiply(a70.subtract(1))).
	add(f3.multiply(a70.multiply(2).subtract(1))).
	add(f4.multiply(d0));
	final T coeff1 = time.getField().getZero();
	final T coeff2 = f1.multiply(a72).
	subtract(f2.multiply(a72)).
	add(f3.multiply(a72.multiply(2))).
	add(f4.multiply(d2));
	final T coeff3 = f1.multiply(a73).
	subtract(f2.multiply(a73)).
	add(f3.multiply(a73.multiply(2))).
	add(f4.multiply(d3));
	final T coeff4 = f1.multiply(a74).
	subtract(f2.multiply(a74)).
	add(f3.multiply(a74.multiply(2))).
	add(f4.multiply(d4));
	final T coeff5 = f1.multiply(a75).
	subtract(f2.multiply(a75)).
	add(f3.multiply(a75.multiply(2))).
	add(f4.multiply(d5));
	final T coeff6 = f4.multiply(d6).subtract(f3);
	final T coeffDot0 = a70.
	subtract(dot2.multiply(a70.subtract(1))).
	add(dot3.multiply(a70.multiply(2).subtract(1))).
	add(dot4.multiply(d0));
	final T coeffDot1 = time.getField().getZero();
	final T coeffDot2 = a72.
	subtract(dot2.multiply(a72)).
	add(dot3.multiply(a72.multiply(2))).
	add(dot4.multiply(d2));
	final T coeffDot3 = a73.
	subtract(dot2.multiply(a73)).
	add(dot3.multiply(a73.multiply(2))).
	add(dot4.multiply(d3));
	final T coeffDot4 = a74.
	subtract(dot2.multiply(a74)).
	add(dot3.multiply(a74.multiply(2))).
	add(dot4.multiply(d4));
	final T coeffDot5 = a75.
	subtract(dot2.multiply(a75)).
	add(dot3.multiply(a75.multiply(2))).
	add(dot4.multiply(d5));
	final T coeffDot6 = dot4.multiply(d6).subtract(dot3);
	interpolatedState = previousStateLinearCombination(coeff0, coeff1, coeff2, coeff3,
	coeff4, coeff5, coeff6);
	interpolatedDerivatives = derivativeLinearCombination(coeffDot0, coeffDot1, coeffDot2, coeffDot3,
	coeffDot4, coeffDot5, coeffDot6);
	} else {
	final T f1 = oneMinusThetaH.negate();
	final T f2 = oneMinusThetaH.multiply(theta);
	final T f3 = f2.multiply(theta);
	final T f4 = f3.multiply(eta);
	final T coeff0 = f1.multiply(a70).
	subtract(f2.multiply(a70.subtract(1))).
	add(f3.multiply(a70.multiply(2).subtract(1))).
	add(f4.multiply(d0));
	final T coeff1 = time.getField().getZero();
	final T coeff2 = f1.multiply(a72).
	subtract(f2.multiply(a72)).
	add(f3.multiply(a72.multiply(2))).
	add(f4.multiply(d2));
	final T coeff3 = f1.multiply(a73).
	subtract(f2.multiply(a73)).
	add(f3.multiply(a73.multiply(2))).
	add(f4.multiply(d3));
	final T coeff4 = f1.multiply(a74).
	subtract(f2.multiply(a74)).
	add(f3.multiply(a74.multiply(2))).
	add(f4.multiply(d4));
	final T coeff5 = f1.multiply(a75).
	subtract(f2.multiply(a75)).
	add(f3.multiply(a75.multiply(2))).
	add(f4.multiply(d5));
	final T coeff6 = f4.multiply(d6).subtract(f3);
	final T coeffDot0 = a70.
	subtract(dot2.multiply(a70.subtract(1))).
	add(dot3.multiply(a70.multiply(2).subtract(1))).
	add(dot4.multiply(d0));
	final T coeffDot1 = time.getField().getZero();
	final T coeffDot2 = a72.
	subtract(dot2.multiply(a72)).
	add(dot3.multiply(a72.multiply(2))).
	add(dot4.multiply(d2));
	final T coeffDot3 = a73.
	subtract(dot2.multiply(a73)).
	add(dot3.multiply(a73.multiply(2))).
	add(dot4.multiply(d3));
	final T coeffDot4 = a74.
	subtract(dot2.multiply(a74)).
	add(dot3.multiply(a74.multiply(2))).
	add(dot4.multiply(d4));
	final T coeffDot5 = a75.
	subtract(dot2.multiply(a75)).
	add(dot3.multiply(a75.multiply(2))).
	add(dot4.multiply(d5));
	final T coeffDot6 = dot4.multiply(d6).subtract(dot3);
	interpolatedState = currentStateLinearCombination(coeff0, coeff1, coeff2, coeff3,
	coeff4, coeff5, coeff6);
	interpolatedDerivatives = derivativeLinearCombination(coeffDot0, coeffDot1, coeffDot2, coeffDot3,
	coeffDot4, coeffDot5, coeffDot6);
	}
	return new FieldODEStateAndDerivative<>(time, interpolatedState, interpolatedDerivatives);

	}

	}