/* * 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.math3.ode; import org.apache.commons.math3.Field; import org.apache.commons.math3.RealFieldElement; import org.apache.commons.math3.exception.DimensionMismatchException; import org.apache.commons.math3.exception.MathIllegalStateException; import org.apache.commons.math3.exception.MaxCountExceededException; import org.apache.commons.math3.exception.NoBracketingException; import org.apache.commons.math3.exception.NumberIsTooSmallException; import org.apache.commons.math3.exception.util.LocalizedFormats; import org.apache.commons.math3.linear.Array2DRowFieldMatrix; import org.apache.commons.math3.ode.nonstiff.AdaptiveStepsizeFieldIntegrator; import org.apache.commons.math3.ode.nonstiff.DormandPrince853FieldIntegrator; import org.apache.commons.math3.ode.sampling.FieldStepHandler; import org.apache.commons.math3.ode.sampling.FieldStepInterpolator; import org.apache.commons.math3.util.FastMath; import org.apache.commons.math3.util.MathArrays; import org.apache.commons.math3.util.MathUtils; /** * This class is the base class for multistep integrators for Ordinary Differential Equations. * *
We define scaled derivatives si(n) at step n as: * *
* s1(n) = h y'n for first derivative * s2(n) = h2/2 y''n for second derivative * s3(n) = h3/6 y'''n for third derivative * ... * sk(n) = hk/k! y(k)n for kth derivative ** *
Rather than storing several previous steps separately, this implementation uses the Nordsieck * vector with higher degrees scaled derivatives all taken at the same step (yn, * s1(n) and rn) where rn is defined as: * *
* rn = [ s2(n), s3(n) ... sk(n) ]T ** * (we omit the k index in the notation for clarity) * *
Multistep integrators with Nordsieck representation are highly sensitive to large step changes
* because when the step is multiplied by factor a, the kth component of the Nordsieck
* vector is multiplied by ak and the last components are the least accurate ones. The
* default max growth factor is therefore set to a quite low value: 21/order.
*
* @see org.apache.commons.math3.ode.nonstiff.AdamsBashforthFieldIntegrator
* @see org.apache.commons.math3.ode.nonstiff.AdamsMoultonFieldIntegrator
* @param (h2/2 y'', h3/6 y''' ..., hk/k! y(k))
*/
protected Array2DRowFieldMatrix The default starter integrator is set to the {@link DormandPrince853FieldIntegrator
* Dormand-Prince 8(5,3)} integrator with some defaults settings.
*
* The default max growth factor is set to a quite low value: 21/order.
*
* @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 The default starter integrator is set to the {@link DormandPrince853FieldIntegrator
* Dormand-Prince 8(5,3)} integrator with some defaults settings.
*
* The default max growth factor is set to a quite low value: 21/order.
*
* @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 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.
*
* @param starterIntegrator starter integrator
*/
public void setStarterIntegrator(FirstOrderFieldIntegrator 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.
*
* @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 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.
*
* @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 MathUtils.min(
error.getField().getZero().add(maxGrowth),
MathUtils.max(
error.getField().getZero().add(minReduction),
error.pow(exp).multiply(safety)));
}
/** Specialized step handler storing the first step. */
private class FieldNordsieckInitializer implements FieldStepHandlert0
* 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