AdamsFieldIntegrator.java
/*
* 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
*
* 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.
*/
/*
* This is not the original file distributed by the Apache Software Foundation
* It has been modified by the Hipparchus project
*/
package org.hipparchus.ode.nonstiff;
import org.hipparchus.Field;
import org.hipparchus.CalculusFieldElement;
import org.hipparchus.exception.MathIllegalArgumentException;
import org.hipparchus.exception.MathIllegalStateException;
import org.hipparchus.linear.Array2DRowFieldMatrix;
import org.hipparchus.linear.FieldMatrix;
import org.hipparchus.ode.FieldEquationsMapper;
import org.hipparchus.ode.FieldExpandableODE;
import org.hipparchus.ode.FieldODEState;
import org.hipparchus.ode.FieldODEStateAndDerivative;
import org.hipparchus.ode.LocalizedODEFormats;
import org.hipparchus.ode.MultistepFieldIntegrator;
import org.hipparchus.util.MathArrays;
/** Base class for {@link AdamsBashforthFieldIntegrator Adams-Bashforth} and
* {@link AdamsMoultonFieldIntegrator Adams-Moulton} integrators.
* @param <T> the type of the field elements
*/
public abstract class AdamsFieldIntegrator<T extends CalculusFieldElement<T>> extends MultistepFieldIntegrator<T> {
/** Transformer. */
private final AdamsNordsieckFieldTransformer<T> transformer;
/**
* Build an Adams integrator with the given order and step control parameters.
* @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 method excluding the one being computed
* @param order order of the method
* @param minStep minimal step (sign is irrelevant, regardless of
* integration direction, forward or backward), the last step can
* be smaller than this
* @param maxStep maximal step (sign is irrelevant, regardless of
* integration direction, forward or backward), the last step can
* be smaller than this
* @param scalAbsoluteTolerance allowed absolute error
* @param scalRelativeTolerance allowed relative error
* @exception MathIllegalArgumentException if order is 1 or less
*/
public AdamsFieldIntegrator(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 MathIllegalArgumentException {
super(field, name, nSteps, order, minStep, maxStep,
scalAbsoluteTolerance, scalRelativeTolerance);
transformer = AdamsNordsieckFieldTransformer.getInstance(field, nSteps);
}
/**
* Build an Adams integrator with the given order and step control parameters.
* @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 method excluding the one being computed
* @param order order of the method
* @param minStep minimal step (sign is irrelevant, regardless of
* integration direction, forward or backward), the last step can
* be smaller than this
* @param maxStep maximal step (sign is irrelevant, regardless of
* integration direction, forward or backward), the last step can
* be smaller than this
* @param vecAbsoluteTolerance allowed absolute error
* @param vecRelativeTolerance allowed relative error
* @exception IllegalArgumentException if order is 1 or less
*/
public AdamsFieldIntegrator(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)
throws IllegalArgumentException {
super(field, name, nSteps, order, minStep, maxStep,
vecAbsoluteTolerance, vecRelativeTolerance);
transformer = AdamsNordsieckFieldTransformer.getInstance(field, nSteps);
}
/** {@inheritDoc} */
@Override
public FieldODEStateAndDerivative<T> integrate(final FieldExpandableODE<T> equations,
final FieldODEState<T> initialState,
final T finalTime)
throws MathIllegalArgumentException, MathIllegalStateException {
sanityChecks(initialState, finalTime);
setStepStart(initIntegration(equations, initialState, finalTime));
final boolean forward = finalTime.subtract(initialState.getTime()).getReal() > 0;
// compute the initial Nordsieck vector using the configured starter integrator
start(equations, getStepStart(), finalTime);
// reuse the step that was chosen by the starter integrator
FieldODEStateAndDerivative<T> stepStart = getStepStart();
FieldODEStateAndDerivative<T> stepEnd =
AdamsFieldStateInterpolator.taylor(equations.getMapper(), stepStart,
stepStart.getTime().add(getStepSize()),
getStepSize(), scaled, nordsieck);
// main integration loop
setIsLastStep(false);
final T[] y = stepStart.getCompleteState();
do {
T[] predictedY = null;
final T[] predictedScaled = MathArrays.buildArray(getField(), y.length);
Array2DRowFieldMatrix<T> predictedNordsieck = null;
double error = 10;
while (error >= 1.0) {
// predict a first estimate of the state at step end
predictedY = stepEnd.getCompleteState();
// evaluate the derivative
final T[] yDot = computeDerivatives(stepEnd.getTime(), predictedY);
// predict Nordsieck vector at step end
for (int j = 0; j < predictedScaled.length; ++j) {
predictedScaled[j] = getStepSize().multiply(yDot[j]);
}
predictedNordsieck = updateHighOrderDerivativesPhase1(nordsieck);
updateHighOrderDerivativesPhase2(scaled, predictedScaled, predictedNordsieck);
// evaluate error
error = errorEstimation(y, stepEnd.getTime(), predictedY, predictedScaled, predictedNordsieck);
if (Double.isNaN(error)) {
throw new MathIllegalStateException(LocalizedODEFormats.NAN_APPEARING_DURING_INTEGRATION,
stepEnd.getTime().getReal());
}
if (error >= 1.0) {
// reject the step and attempt to reduce error by stepsize control
final double factor = computeStepGrowShrinkFactor(error);
rescale(getStepSizeHelper().filterStep(getStepSize().multiply(factor), forward, false));
stepEnd = AdamsFieldStateInterpolator.taylor(equations.getMapper(), getStepStart(),
getStepStart().getTime().add(getStepSize()),
getStepSize(),
scaled,
nordsieck);
}
}
final AdamsFieldStateInterpolator<T> interpolator =
finalizeStep(getStepSize(), predictedY, predictedScaled, predictedNordsieck,
forward, getStepStart(), stepEnd, equations.getMapper());
// discrete events handling
setStepStart(acceptStep(interpolator, finalTime));
scaled = interpolator.getScaled();
nordsieck = interpolator.getNordsieck();
if (!isLastStep()) {
if (resetOccurred()) {
// some events handler has triggered changes that
// invalidate the derivatives, we need to restart from scratch
start(equations, getStepStart(), finalTime);
final T nextT = getStepStart().getTime().add(getStepSize());
final boolean nextIsLast = forward ?
nextT.subtract(finalTime).getReal() >= 0 :
nextT.subtract(finalTime).getReal() <= 0;
final T hNew = nextIsLast ? finalTime.subtract(getStepStart().getTime()) : getStepSize();
rescale(hNew);
System.arraycopy(getStepStart().getCompleteState(), 0, y, 0, y.length);
} else {
// stepsize control for next step
final double factor = computeStepGrowShrinkFactor(error);
final T scaledH = getStepSize().multiply(factor);
final T nextT = getStepStart().getTime().add(scaledH);
final boolean nextIsLast = forward ?
nextT.subtract(finalTime).getReal() >= 0 :
nextT.subtract(finalTime).getReal() <= 0;
T hNew = getStepSizeHelper().filterStep(scaledH, forward, nextIsLast);
final T filteredNextT = getStepStart().getTime().add(hNew);
final boolean filteredNextIsLast = forward ?
filteredNextT.subtract(finalTime).getReal() >= 0 :
filteredNextT.subtract(finalTime).getReal() <= 0;
if (filteredNextIsLast) {
hNew = finalTime.subtract(getStepStart().getTime());
}
rescale(hNew);
System.arraycopy(predictedY, 0, y, 0, y.length);
}
stepEnd = AdamsFieldStateInterpolator.taylor(equations.getMapper(), getStepStart(),
getStepStart().getTime().add(getStepSize()),
getStepSize(), scaled, nordsieck);
}
} while (!isLastStep());
final FieldODEStateAndDerivative<T> finalState = getStepStart();
setStepStart(null);
setStepSize(null);
return finalState;
}
/** {@inheritDoc} */
@Override
protected Array2DRowFieldMatrix<T> initializeHighOrderDerivatives(final T h, final T[] t,
final T[][] y,
final T[][] yDot) {
return transformer.initializeHighOrderDerivatives(h, t, y, yDot);
}
/** Update the high order scaled derivatives for Adams integrators (phase 1).
* <p>The complete update of high order derivatives has a form similar to:
* \[
* r_{n+1} = (s_1(n) - s_1(n+1)) P^{-1} u + P^{-1} A P r_n
* \]
* this method computes the P<sup>-1</sup> A P r<sub>n</sub> part.</p>
* @param highOrder high order scaled derivatives
* (h<sup>2</sup>/2 y'', ... h<sup>k</sup>/k! y(k))
* @return updated high order derivatives
* @see #updateHighOrderDerivativesPhase2(CalculusFieldElement[], CalculusFieldElement[], Array2DRowFieldMatrix)
*/
public Array2DRowFieldMatrix<T> updateHighOrderDerivativesPhase1(final Array2DRowFieldMatrix<T> highOrder) {
return transformer.updateHighOrderDerivativesPhase1(highOrder);
}
/** Update the high order scaled derivatives Adams integrators (phase 2).
* <p>The complete update of high order derivatives has a form similar to:
* \[
* r_{n+1} = (s_1(n) - s_1(n+1)) P^{-1} u + P^{-1} A P r_n
* \]
* this method computes the (s<sub>1</sub>(n) - s<sub>1</sub>(n+1)) P<sup>-1</sup> u part.</p>
* <p>Phase 1 of the update must already have been performed.</p>
* @param start first order scaled derivatives at step start
* @param end first order scaled derivatives at step end
* @param highOrder high order scaled derivatives, will be modified
* (h<sup>2</sup>/2 y'', ... h<sup>k</sup>/k! y(k))
* @see #updateHighOrderDerivativesPhase1(Array2DRowFieldMatrix)
*/
public void updateHighOrderDerivativesPhase2(final T[] start, final T[] end,
final Array2DRowFieldMatrix<T> highOrder) {
transformer.updateHighOrderDerivativesPhase2(start, end, highOrder);
}
/** Estimate error.
* @param previousState state vector at step start
* @param predictedTime time at step end
* @param predictedState predicted state vector at step end
* @param predictedScaled predicted value of the scaled derivatives at step end
* @param predictedNordsieck predicted value of the Nordsieck vector at step end
* @return estimated normalized local discretization error
* @since 2.0
*/
protected abstract double errorEstimation(T[] previousState, T predictedTime,
T[] predictedState, T[] predictedScaled,
FieldMatrix<T> predictedNordsieck);
/** Finalize the step.
* @param stepSize step size used in the scaled and Nordsieck arrays
* @param predictedState predicted state at end of step
* @param predictedScaled predicted first scaled derivative
* @param predictedNordsieck predicted Nordsieck vector
* @param isForward integration direction indicator
* @param globalPreviousState start of the global step
* @param globalCurrentState end of the global step
* @param equationsMapper mapper for ODE equations primary and secondary components
* @return step interpolator
* @since 2.0
*/
protected abstract AdamsFieldStateInterpolator<T> finalizeStep(T stepSize, T[] predictedState,
T[] predictedScaled, Array2DRowFieldMatrix<T> predictedNordsieck,
boolean isForward,
FieldODEStateAndDerivative<T> globalPreviousState,
FieldODEStateAndDerivative<T> globalCurrentState,
FieldEquationsMapper<T> equationsMapper);
}