commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/ode/MultistepFieldIntegrator.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;

 import org.apache.commons.math4.legacy.core.Field;
 import org.apache.commons.math4.legacy.core.RealFieldElement;
 import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
 import org.apache.commons.math4.legacy.exception.MathIllegalStateException;
 import org.apache.commons.math4.legacy.exception.MaxCountExceededException;
 import org.apache.commons.math4.legacy.exception.NoBracketingException;
 import org.apache.commons.math4.legacy.exception.NumberIsTooSmallException;
 import org.apache.commons.math4.legacy.exception.util.LocalizedFormats;
 import org.apache.commons.math4.legacy.linear.Array2DRowFieldMatrix;
 import org.apache.commons.math4.legacy.ode.nonstiff.AdaptiveStepsizeFieldIntegrator;
 import org.apache.commons.math4.legacy.ode.nonstiff.DormandPrince853FieldIntegrator;
 import org.apache.commons.math4.legacy.ode.sampling.FieldStepHandler;
 import org.apache.commons.math4.legacy.ode.sampling.FieldStepInterpolator;
 import org.apache.commons.math4.core.jdkmath.JdkMath;
 import org.apache.commons.math4.legacy.core.MathArrays;

 /**
  * This class is the base class for multistep integrators for Ordinary
  * Differential Equations.
  * <p>We define scaled derivatives s<sub>i</sub>(n) at step n as:
  * <div style="white-space: pre"><code>
  * s<sub>1</sub>(n) = h y'<sub>n</sub> for first derivative
  * s<sub>2</sub>(n) = h<sup>2</sup>/2 y''<sub>n</sub> for second derivative
  * s<sub>3</sub>(n) = h<sup>3</sup>/6 y'''<sub>n</sub> for third derivative
  * ...
  * s<sub>k</sub>(n) = h<sup>k</sup>/k! y<sup>(k)</sup><sub>n</sub> for k<sup>th</sup> derivative
  * </code></div>
  * <p>Rather than storing several previous steps separately, this implementation uses
  * the Nordsieck vector with higher degrees scaled derivatives all taken at the same
  * step (y<sub>n</sub>, s<sub>1</sub>(n) and r<sub>n</sub>) where r<sub>n</sub> is defined as:
  * <div style="white-space: pre"><code>
  * r<sub>n</sub> = [ s<sub>2</sub>(n), s<sub>3</sub>(n) ... s<sub>k</sub>(n) ]<sup>T</sup>
  * </code></div>
  * (we omit the k index in the notation for clarity)
  * <p>
  * Multistep integrators with Nordsieck representation are highly sensitive to
  * large step changes because when the step is multiplied by factor a, the
  * k<sup>th</sup> component of the Nordsieck vector is multiplied by a<sup>k</sup>
  * and the last components are the least accurate ones. The default max growth
  * factor is therefore set to a quite low value: 2<sup>1/order</sup>.
  * </p>
  *
  * @see org.apache.commons.math4.legacy.ode.nonstiff.AdamsBashforthFieldIntegrator
  * @see org.apache.commons.math4.legacy.ode.nonstiff.AdamsMoultonFieldIntegrator
  * @param <T> the type of the field elements
  * @since 3.6
  */
 public abstract class MultistepFieldIntegrator<T extends RealFieldElement<T>>
     extends AdaptiveStepsizeFieldIntegrator<T> {

     /** First scaled derivative (h y'). */
     protected T[] scaled;

     /** Nordsieck matrix of the higher scaled derivatives.
      * <p>(h<sup>2</sup>/2 y'', h<sup>3</sup>/6 y''' ..., h<sup>k</sup>/k! y<sup>(k)</sup>)</p>
      */
     protected Array2DRowFieldMatrix<T> nordsieck;

     /** Starter integrator. */
     private FirstOrderFieldIntegrator<T> starter;

     /** Number of steps of the multistep method (excluding the one being computed). */
     private final int nSteps;

     /** Stepsize control exponent. */
     private double exp;

     /** Safety factor for stepsize control. */
     private double safety;

     /** Minimal reduction factor for stepsize control. */
     private double minReduction;

     /** Maximal growth factor for stepsize control. */
     private double maxGrowth;

     /**
      * Build a multistep integrator with the given stepsize bounds.
      * <p>The default starter integrator is set to the {@link
      * DormandPrince853FieldIntegrator Dormand-Prince 8(5,3)} integrator with
      * some defaults settings.</p>
      * <p>
      * The default max growth factor is set to a quite low value: 2<sup>1/order</sup>.
      * </p>
      * @param field field to which the time and state vector elements belong
      * @param name name of the method
      * @param nSteps number of steps of the multistep method
      * (excluding the one being computed)
      * @param order order of the method
      * @param minStep minimal step (must be positive even for backward
      * integration), the last step can be smaller than this
      * @param maxStep maximal step (must be positive even for backward
      * integration)
      * @param scalAbsoluteTolerance allowed absolute error
      * @param scalRelativeTolerance allowed relative error
      * @exception NumberIsTooSmallException if number of steps is smaller than 2
      */
     protected MultistepFieldIntegrator(final Field<T> field, final String name,
                                        final int nSteps, final int order,
                                        final double minStep, final double maxStep,
                                        final double scalAbsoluteTolerance,
                                        final double scalRelativeTolerance)
         throws NumberIsTooSmallException {

         super(field, name, minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);

         if (nSteps < 2) {
             throw new NumberIsTooSmallException(
                   LocalizedFormats.INTEGRATION_METHOD_NEEDS_AT_LEAST_TWO_PREVIOUS_POINTS,
                   nSteps, 2, true);
         }

         starter = new DormandPrince853FieldIntegrator<>(field, minStep, maxStep,
                                                          scalAbsoluteTolerance,
                                                          scalRelativeTolerance);
         this.nSteps = nSteps;

         exp = -1.0 / order;

         // set the default values of the algorithm control parameters
         setSafety(0.9);
         setMinReduction(0.2);
         setMaxGrowth(JdkMath.pow(2.0, -exp));

     }

     /**
      * Build a multistep integrator with the given stepsize bounds.
      * <p>The default starter integrator is set to the {@link
      * DormandPrince853FieldIntegrator Dormand-Prince 8(5,3)} integrator with
      * some defaults settings.</p>
      * <p>
      * The default max growth factor is set to a quite low value: 2<sup>1/order</sup>.
      * </p>
      * @param field field to which the time and state vector elements belong
      * @param name name of the method
      * @param nSteps number of steps of the multistep method
      * (excluding the one being computed)
      * @param order order of the method
      * @param minStep minimal step (must be positive even for backward
      * integration), the last step can be smaller than this
      * @param maxStep maximal step (must be positive even for backward
      * integration)
      * @param vecAbsoluteTolerance allowed absolute error
      * @param vecRelativeTolerance allowed relative error
      */
     protected MultistepFieldIntegrator(final Field<T> field, final String name, final int nSteps,
                                        final int order,
                                        final double minStep, final double maxStep,
                                        final double[] vecAbsoluteTolerance,
                                        final double[] vecRelativeTolerance) {
         super(field, name, minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
         starter = new DormandPrince853FieldIntegrator<>(field, minStep, maxStep,
                                                          vecAbsoluteTolerance,
                                                          vecRelativeTolerance);
         this.nSteps = nSteps;

         exp = -1.0 / order;

         // set the default values of the algorithm control parameters
         setSafety(0.9);
         setMinReduction(0.2);
         setMaxGrowth(JdkMath.pow(2.0, -exp));

     }

     /**
      * Get the starter integrator.
      * @return starter integrator
      */
     public FirstOrderFieldIntegrator<T> getStarterIntegrator() {
         return starter;
     }

     /**
      * Set the starter integrator.
      * <p>The various step and event handlers for this starter integrator
      * will be managed automatically by the multi-step integrator. Any
      * user configuration for these elements will be cleared before use.</p>
      * @param starterIntegrator starter integrator
      */
     public void setStarterIntegrator(FirstOrderFieldIntegrator<T> starterIntegrator) {
         this.starter = starterIntegrator;
     }

     /** Start the integration.
      * <p>This method computes one step using the underlying starter integrator,
      * and initializes the Nordsieck vector at step start. The starter integrator
      * purpose is only to establish initial conditions, it does not really change
      * time by itself. The top level multistep integrator remains in charge of
      * handling time propagation and events handling as it will starts its own
      * computation right from the beginning. In a sense, the starter integrator
      * can be seen as a dummy one and so it will never trigger any user event nor
      * call any user step handler.</p>
      * @param equations complete set of differential equations to integrate
      * @param initialState initial state (time, primary and secondary state vectors)
      * @param t target time for the integration
      * (can be set to a value smaller than <code>t0</code> for backward integration)
      * @exception DimensionMismatchException if arrays dimension do not match equations settings
      * @exception NumberIsTooSmallException if integration step is too small
      * @exception MaxCountExceededException if the number of functions evaluations is exceeded
      * @exception NoBracketingException if the location of an event cannot be bracketed
      */
     protected void start(final FieldExpandableODE<T> equations, final FieldODEState<T> initialState, final T t)
         throws DimensionMismatchException, NumberIsTooSmallException,
                MaxCountExceededException, NoBracketingException {

         // make sure NO user event nor user step handler is triggered,
         // this is the task of the top level integrator, not the task
         // of the starter integrator
         starter.clearEventHandlers();
         starter.clearStepHandlers();

         // set up one specific step handler to extract initial Nordsieck vector
         starter.addStepHandler(new FieldNordsieckInitializer(equations.getMapper(), (nSteps + 3) / 2));

         // start integration, expecting a InitializationCompletedMarkerException
         try {

             starter.integrate(equations, initialState, t);

             // we should not reach this step
             throw new MathIllegalStateException(LocalizedFormats.MULTISTEP_STARTER_STOPPED_EARLY);

         } catch (InitializationCompletedMarkerException icme) { // NOPMD
             // this is the expected nominal interruption of the start integrator

             // count the evaluations used by the starter
             getEvaluationsCounter().increment(starter.getEvaluations());

         }

         // remove the specific step handler
         starter.clearStepHandlers();

     }

     /** Initialize the high order scaled derivatives at step start.
      * @param h step size to use for scaling
      * @param t first steps times
      * @param y first steps states
      * @param yDot first steps derivatives
      * @return Nordieck vector at first step (h<sup>2</sup>/2 y''<sub>n</sub>,
      * h<sup>3</sup>/6 y'''<sub>n</sub> ... h<sup>k</sup>/k! y<sup>(k)</sup><sub>n</sub>)
      */
     protected abstract Array2DRowFieldMatrix<T> initializeHighOrderDerivatives(T h, T[] t,
                                                                                T[][] y,
                                                                                T[][] yDot);

     /** Get the minimal reduction factor for stepsize control.
      * @return minimal reduction factor
      */
     public double getMinReduction() {
         return minReduction;
     }

     /** Set the minimal reduction factor for stepsize control.
      * @param minReduction minimal reduction factor
      */
     public void setMinReduction(final double minReduction) {
         this.minReduction = minReduction;
     }

     /** Get the maximal growth factor for stepsize control.
      * @return maximal growth factor
      */
     public double getMaxGrowth() {
         return maxGrowth;
     }

     /** Set the maximal growth factor for stepsize control.
      * @param maxGrowth maximal growth factor
      */
     public void setMaxGrowth(final double maxGrowth) {
         this.maxGrowth = maxGrowth;
     }

     /** Get the safety factor for stepsize control.
      * @return safety factor
      */
     public double getSafety() {
       return safety;
     }

     /** Set the safety factor for stepsize control.
      * @param safety safety factor
      */
     public void setSafety(final double safety) {
       this.safety = safety;
     }

     /** Get the number of steps of the multistep method (excluding the one being computed).
      * @return number of steps of the multistep method (excluding the one being computed)
      */
     public int getNSteps() {
       return nSteps;
     }

     /** Rescale the instance.
      * <p>Since the scaled and Nordsieck arrays are shared with the caller,
      * this method has the side effect of rescaling this arrays in the caller too.</p>
      * @param newStepSize new step size to use in the scaled and Nordsieck arrays
      */
     protected void rescale(final T newStepSize) {

         final T ratio = newStepSize.divide(getStepSize());
         for (int i = 0; i < scaled.length; ++i) {
             scaled[i] = scaled[i].multiply(ratio);
         }

         final T[][] nData = nordsieck.getDataRef();
         T power = ratio;
         for (int i = 0; i < nData.length; ++i) {
             power = power.multiply(ratio);
             final T[] nDataI = nData[i];
             for (int j = 0; j < nDataI.length; ++j) {
                 nDataI[j] = nDataI[j].multiply(power);
             }
         }

         setStepSize(newStepSize);

     }


     /** Compute step grow/shrink factor according to normalized error.
      * @param error normalized error of the current step
      * @return grow/shrink factor for next step
      */
     protected T computeStepGrowShrinkFactor(final T error) {
         return RealFieldElement.min(error.getField().getZero().add(maxGrowth),
                              RealFieldElement.max(error.getField().getZero().add(minReduction),
                                            error.pow(exp).multiply(safety)));
     }

     /** Specialized step handler storing the first step.
      */
     private class FieldNordsieckInitializer implements FieldStepHandler<T> {

         /** Equation mapper. */
         private final FieldEquationsMapper<T> mapper;

         /** Steps counter. */
         private int count;

         /** Saved start. */
         private FieldODEStateAndDerivative<T> savedStart;

         /** First steps times. */
         private final T[] t;

         /** First steps states. */
         private final T[][] y;

         /** First steps derivatives. */
         private final T[][] yDot;

         /** Simple constructor.
          * @param mapper equation mapper
          * @param nbStartPoints number of start points (including the initial point)
          */
         FieldNordsieckInitializer(final FieldEquationsMapper<T> mapper, final int nbStartPoints) {
             this.mapper = mapper;
             this.count  = 0;
             this.t      = MathArrays.buildArray(getField(), nbStartPoints);
             this.y      = MathArrays.buildArray(getField(), nbStartPoints, -1);
             this.yDot   = MathArrays.buildArray(getField(), nbStartPoints, -1);
         }

         /** {@inheritDoc} */
         @Override
         public void handleStep(FieldStepInterpolator<T> interpolator, boolean isLast)
             throws MaxCountExceededException {


             if (count == 0) {
                 // first step, we need to store also the point at the beginning of the step
                 final FieldODEStateAndDerivative<T> prev = interpolator.getPreviousState();
                 savedStart  = prev;
                 t[count]    = prev.getTime();
                 y[count]    = mapper.mapState(prev);
                 yDot[count] = mapper.mapDerivative(prev);
             }

             // store the point at the end of the step
             ++count;
             final FieldODEStateAndDerivative<T> curr = interpolator.getCurrentState();
             t[count]    = curr.getTime();
             y[count]    = mapper.mapState(curr);
             yDot[count] = mapper.mapDerivative(curr);

             if (count == t.length - 1) {

                 // this was the last point we needed, we can compute the derivatives
                 setStepSize(t[t.length - 1].subtract(t[0]).divide(t.length - 1));

                 // first scaled derivative
                 scaled = MathArrays.buildArray(getField(), yDot[0].length);
                 for (int j = 0; j < scaled.length; ++j) {
                     scaled[j] = yDot[0][j].multiply(getStepSize());
                 }

                 // higher order derivatives
                 nordsieck = initializeHighOrderDerivatives(getStepSize(), t, y, yDot);

                 // stop the integrator now that all needed steps have been handled
                 setStepStart(savedStart);
                 throw new InitializationCompletedMarkerException();

             }

         }

         /** {@inheritDoc} */
         @Override
         public void init(final FieldODEStateAndDerivative<T> initialState, T finalTime) {
             // nothing to do
         }

     }

     /** Marker exception used ONLY to stop the starter integrator after first step. */
     private static class InitializationCompletedMarkerException
         extends RuntimeException {

         /** Serializable version identifier. */
         private static final long serialVersionUID = -1914085471038046418L;

         /** Simple constructor. */
         InitializationCompletedMarkerException() {
             super((Throwable) null);
         }

     }

 }
	/*
	* 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;

	import org.apache.commons.math4.legacy.core.Field;
	import org.apache.commons.math4.legacy.core.RealFieldElement;
	import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
	import org.apache.commons.math4.legacy.exception.MathIllegalStateException;
	import org.apache.commons.math4.legacy.exception.MaxCountExceededException;
	import org.apache.commons.math4.legacy.exception.NoBracketingException;
	import org.apache.commons.math4.legacy.exception.NumberIsTooSmallException;
	import org.apache.commons.math4.legacy.exception.util.LocalizedFormats;
	import org.apache.commons.math4.legacy.linear.Array2DRowFieldMatrix;
	import org.apache.commons.math4.legacy.ode.nonstiff.AdaptiveStepsizeFieldIntegrator;
	import org.apache.commons.math4.legacy.ode.nonstiff.DormandPrince853FieldIntegrator;
	import org.apache.commons.math4.legacy.ode.sampling.FieldStepHandler;
	import org.apache.commons.math4.legacy.ode.sampling.FieldStepInterpolator;
	import org.apache.commons.math4.core.jdkmath.JdkMath;
	import org.apache.commons.math4.legacy.core.MathArrays;

	/**
	* This class is the base class for multistep integrators for Ordinary
	* Differential Equations.
	* <p>We define scaled derivatives s<sub>i</sub>(n) at step n as:
	* <div style="white-space: pre"><code>
	* s<sub>1</sub>(n) = h y'<sub>n</sub> for first derivative
	* s<sub>2</sub>(n) = h<sup>2</sup>/2 y''<sub>n</sub> for second derivative
	* s<sub>3</sub>(n) = h<sup>3</sup>/6 y'''<sub>n</sub> for third derivative
	* ...
	* s<sub>k</sub>(n) = h<sup>k</sup>/k! y<sup>(k)</sup><sub>n</sub> for k<sup>th</sup> derivative
	* </code></div>
	* <p>Rather than storing several previous steps separately, this implementation uses
	* the Nordsieck vector with higher degrees scaled derivatives all taken at the same
	* step (y<sub>n</sub>, s<sub>1</sub>(n) and r<sub>n</sub>) where r<sub>n</sub> is defined as:
	* <div style="white-space: pre"><code>
	* r<sub>n</sub> = [ s<sub>2</sub>(n), s<sub>3</sub>(n) ... s<sub>k</sub>(n) ]<sup>T</sup>
	* </code></div>
	* (we omit the k index in the notation for clarity)
	* <p>
	* Multistep integrators with Nordsieck representation are highly sensitive to
	* large step changes because when the step is multiplied by factor a, the
	* k<sup>th</sup> component of the Nordsieck vector is multiplied by a<sup>k</sup>
	* and the last components are the least accurate ones. The default max growth
	* factor is therefore set to a quite low value: 2<sup>1/order</sup>.
	* </p>
	*
	* @see org.apache.commons.math4.legacy.ode.nonstiff.AdamsBashforthFieldIntegrator
	* @see org.apache.commons.math4.legacy.ode.nonstiff.AdamsMoultonFieldIntegrator
	* @param <T> the type of the field elements
	* @since 3.6
	*/
	public abstract class MultistepFieldIntegrator<T extends RealFieldElement<T>>
	extends AdaptiveStepsizeFieldIntegrator<T> {

	/** First scaled derivative (h y'). */
	protected T[] scaled;

	/** Nordsieck matrix of the higher scaled derivatives.
	* <p>(h<sup>2</sup>/2 y'', h<sup>3</sup>/6 y''' ..., h<sup>k</sup>/k! y<sup>(k)</sup>)</p>
	*/
	protected Array2DRowFieldMatrix<T> nordsieck;

	/** Starter integrator. */
	private FirstOrderFieldIntegrator<T> starter;

	/** Number of steps of the multistep method (excluding the one being computed). */
	private final int nSteps;

	/** Stepsize control exponent. */
	private double exp;

	/** Safety factor for stepsize control. */
	private double safety;

	/** Minimal reduction factor for stepsize control. */
	private double minReduction;

	/** Maximal growth factor for stepsize control. */
	private double maxGrowth;

	/**
	* Build a multistep integrator with the given stepsize bounds.
	* <p>The default starter integrator is set to the {@link
	* DormandPrince853FieldIntegrator Dormand-Prince 8(5,3)} integrator with
	* some defaults settings.</p>
	* <p>
	* The default max growth factor is set to a quite low value: 2<sup>1/order</sup>.
	* </p>
	* @param field field to which the time and state vector elements belong
	* @param name name of the method
	* @param nSteps number of steps of the multistep method
	* (excluding the one being computed)
	* @param order order of the method
	* @param minStep minimal step (must be positive even for backward
	* integration), the last step can be smaller than this
	* @param maxStep maximal step (must be positive even for backward
	* integration)
	* @param scalAbsoluteTolerance allowed absolute error
	* @param scalRelativeTolerance allowed relative error
	* @exception NumberIsTooSmallException if number of steps is smaller than 2
	*/
	protected MultistepFieldIntegrator(final Field<T> field, final String name,
	final int nSteps, final int order,
	final double minStep, final double maxStep,
	final double scalAbsoluteTolerance,
	final double scalRelativeTolerance)
	throws NumberIsTooSmallException {

	super(field, name, minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);

	if (nSteps < 2) {
	throw new NumberIsTooSmallException(
	LocalizedFormats.INTEGRATION_METHOD_NEEDS_AT_LEAST_TWO_PREVIOUS_POINTS,
	nSteps, 2, true);
	}

	starter = new DormandPrince853FieldIntegrator<>(field, minStep, maxStep,
	scalAbsoluteTolerance,
	scalRelativeTolerance);
	this.nSteps = nSteps;

	exp = -1.0 / order;

	// set the default values of the algorithm control parameters
	setSafety(0.9);
	setMinReduction(0.2);
	setMaxGrowth(JdkMath.pow(2.0, -exp));

	}

	/**
	* Build a multistep integrator with the given stepsize bounds.
	* <p>The default starter integrator is set to the {@link
	* DormandPrince853FieldIntegrator Dormand-Prince 8(5,3)} integrator with
	* some defaults settings.</p>
	* <p>
	* The default max growth factor is set to a quite low value: 2<sup>1/order</sup>.
	* </p>
	* @param field field to which the time and state vector elements belong
	* @param name name of the method
	* @param nSteps number of steps of the multistep method
	* (excluding the one being computed)
	* @param order order of the method
	* @param minStep minimal step (must be positive even for backward
	* integration), the last step can be smaller than this
	* @param maxStep maximal step (must be positive even for backward
	* integration)
	* @param vecAbsoluteTolerance allowed absolute error
	* @param vecRelativeTolerance allowed relative error
	*/
	protected MultistepFieldIntegrator(final Field<T> field, final String name, final int nSteps,
	final int order,
	final double minStep, final double maxStep,
	final double[] vecAbsoluteTolerance,
	final double[] vecRelativeTolerance) {
	super(field, name, minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
	starter = new DormandPrince853FieldIntegrator<>(field, minStep, maxStep,
	vecAbsoluteTolerance,
	vecRelativeTolerance);
	this.nSteps = nSteps;

	exp = -1.0 / order;

	// set the default values of the algorithm control parameters
	setSafety(0.9);
	setMinReduction(0.2);
	setMaxGrowth(JdkMath.pow(2.0, -exp));

	}

	/**
	* Get the starter integrator.
	* @return starter integrator
	*/
	public FirstOrderFieldIntegrator<T> getStarterIntegrator() {
	return starter;
	}

	/**
	* Set the starter integrator.
	* <p>The various step and event handlers for this starter integrator
	* will be managed automatically by the multi-step integrator. Any
	* user configuration for these elements will be cleared before use.</p>
	* @param starterIntegrator starter integrator
	*/
	public void setStarterIntegrator(FirstOrderFieldIntegrator<T> starterIntegrator) {
	this.starter = starterIntegrator;
	}

	/** Start the integration.
	* <p>This method computes one step using the underlying starter integrator,
	* and initializes the Nordsieck vector at step start. The starter integrator
	* purpose is only to establish initial conditions, it does not really change
	* time by itself. The top level multistep integrator remains in charge of
	* handling time propagation and events handling as it will starts its own
	* computation right from the beginning. In a sense, the starter integrator
	* can be seen as a dummy one and so it will never trigger any user event nor
	* call any user step handler.</p>
	* @param equations complete set of differential equations to integrate
	* @param initialState initial state (time, primary and secondary state vectors)
	* @param t target time for the integration
	* (can be set to a value smaller than <code>t0</code> for backward integration)
	* @exception DimensionMismatchException if arrays dimension do not match equations settings
	* @exception NumberIsTooSmallException if integration step is too small
	* @exception MaxCountExceededException if the number of functions evaluations is exceeded
	* @exception NoBracketingException if the location of an event cannot be bracketed
	*/
	protected void start(final FieldExpandableODE<T> equations, final FieldODEState<T> initialState, final T t)
	throws DimensionMismatchException, NumberIsTooSmallException,
	MaxCountExceededException, NoBracketingException {

	// make sure NO user event nor user step handler is triggered,
	// this is the task of the top level integrator, not the task
	// of the starter integrator
	starter.clearEventHandlers();
	starter.clearStepHandlers();

	// set up one specific step handler to extract initial Nordsieck vector
	starter.addStepHandler(new FieldNordsieckInitializer(equations.getMapper(), (nSteps + 3) / 2));

	// start integration, expecting a InitializationCompletedMarkerException
	try {

	starter.integrate(equations, initialState, t);

	// we should not reach this step
	throw new MathIllegalStateException(LocalizedFormats.MULTISTEP_STARTER_STOPPED_EARLY);

	} catch (InitializationCompletedMarkerException icme) { // NOPMD
	// this is the expected nominal interruption of the start integrator

	// count the evaluations used by the starter
	getEvaluationsCounter().increment(starter.getEvaluations());

	}

	// remove the specific step handler
	starter.clearStepHandlers();

	}

	/** Initialize the high order scaled derivatives at step start.
	* @param h step size to use for scaling
	* @param t first steps times
	* @param y first steps states
	* @param yDot first steps derivatives
	* @return Nordieck vector at first step (h<sup>2</sup>/2 y''<sub>n</sub>,
	* h<sup>3</sup>/6 y'''<sub>n</sub> ... h<sup>k</sup>/k! y<sup>(k)</sup><sub>n</sub>)
	*/
	protected abstract Array2DRowFieldMatrix<T> initializeHighOrderDerivatives(T h, T[] t,
	T[][] y,
	T[][] yDot);

	/** Get the minimal reduction factor for stepsize control.
	* @return minimal reduction factor
	*/
	public double getMinReduction() {
	return minReduction;
	}

	/** Set the minimal reduction factor for stepsize control.
	* @param minReduction minimal reduction factor
	*/
	public void setMinReduction(final double minReduction) {
	this.minReduction = minReduction;
	}

	/** Get the maximal growth factor for stepsize control.
	* @return maximal growth factor
	*/
	public double getMaxGrowth() {
	return maxGrowth;
	}

	/** Set the maximal growth factor for stepsize control.
	* @param maxGrowth maximal growth factor
	*/
	public void setMaxGrowth(final double maxGrowth) {
	this.maxGrowth = maxGrowth;
	}

	/** Get the safety factor for stepsize control.
	* @return safety factor
	*/
	public double getSafety() {
	return safety;
	}

	/** Set the safety factor for stepsize control.
	* @param safety safety factor
	*/
	public void setSafety(final double safety) {
	this.safety = safety;
	}

	/** Get the number of steps of the multistep method (excluding the one being computed).
	* @return number of steps of the multistep method (excluding the one being computed)
	*/
	public int getNSteps() {
	return nSteps;
	}

	/** Rescale the instance.
	* <p>Since the scaled and Nordsieck arrays are shared with the caller,
	* this method has the side effect of rescaling this arrays in the caller too.</p>
	* @param newStepSize new step size to use in the scaled and Nordsieck arrays
	*/
	protected void rescale(final T newStepSize) {

	final T ratio = newStepSize.divide(getStepSize());
	for (int i = 0; i < scaled.length; ++i) {
	scaled[i] = scaled[i].multiply(ratio);
	}

	final T[][] nData = nordsieck.getDataRef();
	T power = ratio;
	for (int i = 0; i < nData.length; ++i) {
	power = power.multiply(ratio);
	final T[] nDataI = nData[i];
	for (int j = 0; j < nDataI.length; ++j) {
	nDataI[j] = nDataI[j].multiply(power);
	}
	}

	setStepSize(newStepSize);

	}


	/** Compute step grow/shrink factor according to normalized error.
	* @param error normalized error of the current step
	* @return grow/shrink factor for next step
	*/
	protected T computeStepGrowShrinkFactor(final T error) {
	return RealFieldElement.min(error.getField().getZero().add(maxGrowth),
	RealFieldElement.max(error.getField().getZero().add(minReduction),
	error.pow(exp).multiply(safety)));
	}

	/** Specialized step handler storing the first step.
	*/
	private class FieldNordsieckInitializer implements FieldStepHandler<T> {

	/** Equation mapper. */
	private final FieldEquationsMapper<T> mapper;

	/** Steps counter. */
	private int count;

	/** Saved start. */
	private FieldODEStateAndDerivative<T> savedStart;

	/** First steps times. */
	private final T[] t;

	/** First steps states. */
	private final T[][] y;

	/** First steps derivatives. */
	private final T[][] yDot;

	/** Simple constructor.
	* @param mapper equation mapper
	* @param nbStartPoints number of start points (including the initial point)
	*/
	FieldNordsieckInitializer(final FieldEquationsMapper<T> mapper, final int nbStartPoints) {
	this.mapper = mapper;
	this.count = 0;
	this.t = MathArrays.buildArray(getField(), nbStartPoints);
	this.y = MathArrays.buildArray(getField(), nbStartPoints, -1);
	this.yDot = MathArrays.buildArray(getField(), nbStartPoints, -1);
	}

	/** {@inheritDoc} */
	@Override
	public void handleStep(FieldStepInterpolator<T> interpolator, boolean isLast)
	throws MaxCountExceededException {


	if (count == 0) {
	// first step, we need to store also the point at the beginning of the step
	final FieldODEStateAndDerivative<T> prev = interpolator.getPreviousState();
	savedStart = prev;
	t[count] = prev.getTime();
	y[count] = mapper.mapState(prev);
	yDot[count] = mapper.mapDerivative(prev);
	}

	// store the point at the end of the step
	++count;
	final FieldODEStateAndDerivative<T> curr = interpolator.getCurrentState();
	t[count] = curr.getTime();
	y[count] = mapper.mapState(curr);
	yDot[count] = mapper.mapDerivative(curr);

	if (count == t.length - 1) {

	// this was the last point we needed, we can compute the derivatives
	setStepSize(t[t.length - 1].subtract(t[0]).divide(t.length - 1));

	// first scaled derivative
	scaled = MathArrays.buildArray(getField(), yDot[0].length);
	for (int j = 0; j < scaled.length; ++j) {
	scaled[j] = yDot[0][j].multiply(getStepSize());
	}

	// higher order derivatives
	nordsieck = initializeHighOrderDerivatives(getStepSize(), t, y, yDot);

	// stop the integrator now that all needed steps have been handled
	setStepStart(savedStart);
	throw new InitializationCompletedMarkerException();

	}

	}

	/** {@inheritDoc} */
	@Override
	public void init(final FieldODEStateAndDerivative<T> initialState, T finalTime) {
	// nothing to do
	}

	}

	/** Marker exception used ONLY to stop the starter integrator after first step. */
	private static class InitializationCompletedMarkerException
	extends RuntimeException {

	/** Serializable version identifier. */
	private static final long serialVersionUID = -1914085471038046418L;

	/** Simple constructor. */
	InitializationCompletedMarkerException() {
	super((Throwable) null);
	}

	}

	}