You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by lu...@apache.org on 2015/12/28 22:14:43 UTC
[3/4] [math] Base classes for field-based multistep integrators.
Base classes for field-based multistep integrators.
Project: http://git-wip-us.apache.org/repos/asf/commons-math/repo
Commit: http://git-wip-us.apache.org/repos/asf/commons-math/commit/ba6bd8cb
Tree: http://git-wip-us.apache.org/repos/asf/commons-math/tree/ba6bd8cb
Diff: http://git-wip-us.apache.org/repos/asf/commons-math/diff/ba6bd8cb
Branch: refs/heads/MATH_3_X
Commit: ba6bd8cbe193b86d9fa44b1c68c417ab88953e29
Parents: 4a501b6
Author: Luc Maisonobe <lu...@apache.org>
Authored: Mon Dec 28 22:13:48 2015 +0100
Committer: Luc Maisonobe <lu...@apache.org>
Committed: Mon Dec 28 22:13:48 2015 +0100
----------------------------------------------------------------------
.../math3/ode/MultistepFieldIntegrator.java | 427 +++++++++++++++++++
.../ode/nonstiff/AdamsFieldIntegrator.java | 146 +++++++
.../nonstiff/AdamsFieldStepInterpolator.java | 193 +++++++++
.../AdamsNordsieckFieldTransformer.java | 362 ++++++++++++++++
.../AbstractAdamsFieldIntegratorTest.java | 261 ++++++++++++
5 files changed, 1389 insertions(+)
----------------------------------------------------------------------
http://git-wip-us.apache.org/repos/asf/commons-math/blob/ba6bd8cb/src/main/java/org/apache/commons/math3/ode/MultistepFieldIntegrator.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/ode/MultistepFieldIntegrator.java b/src/main/java/org/apache/commons/math3/ode/MultistepFieldIntegrator.java
new file mode 100644
index 0000000..b86fc05
--- /dev/null
+++ b/src/main/java/org/apache/commons/math3/ode/MultistepFieldIntegrator.java
@@ -0,0 +1,427 @@
+/*
+ * 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.
+ * <p>We define scaled derivatives s<sub>i</sub>(n) at step n as:
+ * <pre>
+ * 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
+ * </pre></p>
+ * <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:
+ * <pre>
+ * r<sub>n</sub> = [ s<sub>2</sub>(n), s<sub>3</sub>(n) ... s<sub>k</sub>(n) ]<sup>T</sup>
+ * </pre>
+ * (we omit the k index in the notation for clarity)</p>
+ * <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.math3.ode.nonstiff.AdamsBashforthFieldIntegrator
+ * @see org.apache.commons.math3.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<T>(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(FastMath.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<T>(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(FastMath.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(final T h, final T[] t,
+ final T[][] y,
+ final 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;
+ }
+
+ /** 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 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} */
+ 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} */
+ 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);
+ }
+
+ }
+
+}
http://git-wip-us.apache.org/repos/asf/commons-math/blob/ba6bd8cb/src/main/java/org/apache/commons/math3/ode/nonstiff/AdamsFieldIntegrator.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/AdamsFieldIntegrator.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/AdamsFieldIntegrator.java
new file mode 100644
index 0000000..1657ca5
--- /dev/null
+++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/AdamsFieldIntegrator.java
@@ -0,0 +1,146 @@
+/*
+ * 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.nonstiff;
+
+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.MaxCountExceededException;
+import org.apache.commons.math3.exception.NoBracketingException;
+import org.apache.commons.math3.exception.NumberIsTooSmallException;
+import org.apache.commons.math3.linear.Array2DRowFieldMatrix;
+import org.apache.commons.math3.ode.FieldExpandableODE;
+import org.apache.commons.math3.ode.FieldODEState;
+import org.apache.commons.math3.ode.FieldODEStateAndDerivative;
+import org.apache.commons.math3.ode.MultistepFieldIntegrator;
+
+
+/** Base class for {@link AdamsBashforthFieldIntegrator Adams-Bashforth} and
+ * {@link AdamsMoultonFieldIntegrator Adams-Moulton} integrators.
+ * @param <T> the type of the field elements
+ * @since 3.6
+ */
+public abstract class AdamsFieldIntegrator<T extends RealFieldElement<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 NumberIsTooSmallException 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 NumberIsTooSmallException {
+ 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 abstract FieldODEStateAndDerivative<T> integrate(final FieldExpandableODE<T> equations,
+ final FieldODEState<T> initialState,
+ final T finalTime)
+ throws NumberIsTooSmallException, DimensionMismatchException,
+ MaxCountExceededException, NoBracketingException;
+
+ /** {@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:
+ * <pre>
+ * r<sub>n+1</sub> = (s<sub>1</sub>(n) - s<sub>1</sub>(n+1)) P<sup>-1</sup> u + P<sup>-1</sup> A P r<sub>n</sub>
+ * </pre>
+ * 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(double[], double[], 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:
+ * <pre>
+ * r<sub>n+1</sub> = (s<sub>1</sub>(n) - s<sub>1</sub>(n+1)) P<sup>-1</sup> u + P<sup>-1</sup> A P r<sub>n</sub>
+ * </pre>
+ * 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);
+ }
+
+}
http://git-wip-us.apache.org/repos/asf/commons-math/blob/ba6bd8cb/src/main/java/org/apache/commons/math3/ode/nonstiff/AdamsFieldStepInterpolator.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/AdamsFieldStepInterpolator.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/AdamsFieldStepInterpolator.java
new file mode 100644
index 0000000..9f282c0
--- /dev/null
+++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/AdamsFieldStepInterpolator.java
@@ -0,0 +1,193 @@
+/*
+ * 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.nonstiff;
+
+import java.util.Arrays;
+
+import org.apache.commons.math3.RealFieldElement;
+import org.apache.commons.math3.linear.Array2DRowFieldMatrix;
+import org.apache.commons.math3.ode.FieldEquationsMapper;
+import org.apache.commons.math3.ode.FieldODEStateAndDerivative;
+import org.apache.commons.math3.ode.sampling.AbstractFieldStepInterpolator;
+import org.apache.commons.math3.util.MathArrays;
+
+/**
+ * This class implements an interpolator for Adams integrators using Nordsieck representation.
+ *
+ * <p>This interpolator computes dense output around the current point.
+ * The interpolation equation is based on Taylor series formulas.
+ *
+ * @see AdamsBashforthFieldIntegrator
+ * @see AdamsMoultonFieldIntegrator
+ * @param <T> the type of the field elements
+ * @since 3.6
+ */
+
+class AdamsFieldStepInterpolator<T extends RealFieldElement<T>> extends AbstractFieldStepInterpolator<T> {
+
+ /** Step size used in the first scaled derivative and Nordsieck vector. */
+ private T scalingH;
+
+ /** First scaled derivative. */
+ private final T[] scaled;
+
+ /** Nordsieck vector. */
+ private final Array2DRowFieldMatrix<T> nordsieck;
+
+ /** Simple constructor.
+ * @param stepSize step size used in the scaled and Nordsieck arrays
+ * @param referenceState reference state from which Taylor expansion are estimated
+ * @param scaled first scaled derivative
+ * @param nordsieck Nordsieck vector
+ * @param isForward integration direction indicator
+ * @param equationsMapper mapper for ODE equations primary and secondary components
+ */
+ AdamsFieldStepInterpolator(final T stepSize, final FieldODEStateAndDerivative<T> referenceState,
+ final T[] scaled, final Array2DRowFieldMatrix<T> nordsieck,
+ final boolean isForward, final FieldEquationsMapper<T> equationsMapper) {
+ this(stepSize, scaled, nordsieck, isForward,
+ referenceState, taylor(referenceState, referenceState.getTime().add(stepSize), stepSize, scaled, nordsieck),
+ equationsMapper);
+ }
+
+ /** Simple constructor.
+ * @param stepSize step size used in the scaled and Nordsieck arrays
+ * @param scaled first scaled derivative
+ * @param nordsieck 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
+ */
+ private AdamsFieldStepInterpolator(final T stepSize, final T[] scaled,
+ final Array2DRowFieldMatrix<T> nordsieck,
+ final boolean isForward,
+ final FieldODEStateAndDerivative<T> globalPreviousState,
+ final FieldODEStateAndDerivative<T> globalCurrentState,
+ final FieldEquationsMapper<T> equationsMapper) {
+ this(stepSize, scaled, nordsieck,
+ isForward, globalPreviousState, globalCurrentState,
+ globalPreviousState, globalCurrentState, equationsMapper);
+ }
+
+ /** Simple constructor.
+ * @param stepSize step size used in the scaled and Nordsieck arrays
+ * @param scaled first scaled derivative
+ * @param nordsieck Nordsieck vector
+ * @param isForward integration direction indicator
+ * @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 equationsMapper mapper for ODE equations primary and secondary components
+ */
+ private AdamsFieldStepInterpolator(final T stepSize, final T[] scaled,
+ final Array2DRowFieldMatrix<T> nordsieck,
+ final boolean isForward,
+ final FieldODEStateAndDerivative<T> globalPreviousState,
+ final FieldODEStateAndDerivative<T> globalCurrentState,
+ final FieldODEStateAndDerivative<T> softPreviousState,
+ final FieldODEStateAndDerivative<T> softCurrentState,
+ final FieldEquationsMapper<T> equationsMapper) {
+ super(isForward, globalPreviousState, globalCurrentState,
+ softPreviousState, softCurrentState, equationsMapper);
+ this.scalingH = stepSize;
+ this.scaled = scaled.clone();
+ this.nordsieck = new Array2DRowFieldMatrix<T>(nordsieck.getData(), false);
+ }
+
+ /** Create a new instance.
+ * @param newForward integration direction indicator
+ * @param newGlobalPreviousState start of the global step
+ * @param newGlobalCurrentState end of the global step
+ * @param newSoftPreviousState start of the restricted step
+ * @param newSoftCurrentState end of the restricted step
+ * @param newMapper equations mapper for the all equations
+ * @return a new instance
+ */
+ protected AdamsFieldStepInterpolator<T> create(boolean newForward,
+ FieldODEStateAndDerivative<T> newGlobalPreviousState,
+ FieldODEStateAndDerivative<T> newGlobalCurrentState,
+ FieldODEStateAndDerivative<T> newSoftPreviousState,
+ FieldODEStateAndDerivative<T> newSoftCurrentState,
+ FieldEquationsMapper<T> newMapper) {
+ return new AdamsFieldStepInterpolator<T>(scalingH, scaled, nordsieck,
+ newForward,
+ newGlobalPreviousState, newGlobalCurrentState,
+ newSoftPreviousState, newSoftCurrentState,
+ newMapper);
+
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ protected FieldODEStateAndDerivative<T> computeInterpolatedStateAndDerivatives(final FieldEquationsMapper<T> equationsMapper,
+ final T time, final T theta,
+ final T thetaH, final T oneMinusThetaH) {
+ return taylor(getPreviousState(), time, scalingH, scaled, nordsieck);
+ }
+
+ /** Estimate state by applying Taylor formula.
+ * @param referenceState reference state
+ * @param time time at which state must be estimated
+ * @param stepSize step size used in the scaled and Nordsieck arrays
+ * @param scaled first scaled derivative
+ * @param nordsieck Nordsieck vector
+ * @return estimated state
+ * @param <S> the type of the field elements
+ */
+ private static <S extends RealFieldElement<S>> FieldODEStateAndDerivative<S> taylor(final FieldODEStateAndDerivative<S> referenceState,
+ final S time, final S stepSize,
+ final S[] scaled,
+ final Array2DRowFieldMatrix<S> nordsieck) {
+
+ final S x = time.subtract(referenceState.getTime());
+ final S normalizedAbscissa = x.divide(stepSize);
+
+ S[] stateVariation = MathArrays.buildArray(time.getField(), scaled.length);
+ Arrays.fill(stateVariation, time.getField().getZero());
+ S[] estimatedDerivatives = MathArrays.buildArray(time.getField(), scaled.length);
+ Arrays.fill(estimatedDerivatives, time.getField().getZero());
+
+ // apply Taylor formula from high order to low order,
+ // for the sake of numerical accuracy
+ final S[][] nData = nordsieck.getDataRef();
+ for (int i = nData.length - 1; i >= 0; --i) {
+ final int order = i + 2;
+ final S[] nDataI = nData[i];
+ final S power = normalizedAbscissa.pow(order);
+ for (int j = 0; j < nDataI.length; ++j) {
+ final S d = nDataI[j].multiply(power);
+ stateVariation[j] = stateVariation[j].add(d);
+ estimatedDerivatives[j] = estimatedDerivatives[j].add(d.multiply(order));
+ }
+ }
+
+ S[] estimatedState = referenceState.getState();
+ for (int j = 0; j < stateVariation.length; ++j) {
+ stateVariation[j] = stateVariation[j].add(scaled[j].multiply(normalizedAbscissa));
+ estimatedState[j] = estimatedState[j].add(stateVariation[j]);
+ estimatedDerivatives[j] =
+ estimatedDerivatives[j].add(scaled[j].multiply(normalizedAbscissa)).divide(x);
+ }
+
+ return new FieldODEStateAndDerivative<S>(time, estimatedState, estimatedDerivatives);
+
+ }
+
+}
http://git-wip-us.apache.org/repos/asf/commons-math/blob/ba6bd8cb/src/main/java/org/apache/commons/math3/ode/nonstiff/AdamsNordsieckFieldTransformer.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/AdamsNordsieckFieldTransformer.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/AdamsNordsieckFieldTransformer.java
new file mode 100644
index 0000000..89f8393
--- /dev/null
+++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/AdamsNordsieckFieldTransformer.java
@@ -0,0 +1,362 @@
+/*
+ * 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.nonstiff;
+
+import java.util.Arrays;
+import java.util.HashMap;
+import java.util.Map;
+
+import org.apache.commons.math3.Field;
+import org.apache.commons.math3.RealFieldElement;
+import org.apache.commons.math3.linear.Array2DRowFieldMatrix;
+import org.apache.commons.math3.linear.ArrayFieldVector;
+import org.apache.commons.math3.linear.FieldDecompositionSolver;
+import org.apache.commons.math3.linear.FieldLUDecomposition;
+import org.apache.commons.math3.linear.FieldMatrix;
+import org.apache.commons.math3.util.MathArrays;
+
+/** Transformer to Nordsieck vectors for Adams integrators.
+ * <p>This class is used by {@link AdamsBashforthIntegrator Adams-Bashforth} and
+ * {@link AdamsMoultonIntegrator Adams-Moulton} integrators to convert between
+ * classical representation with several previous first derivatives and Nordsieck
+ * representation with higher order scaled derivatives.</p>
+ *
+ * <p>We define scaled derivatives s<sub>i</sub>(n) at step n as:
+ * <pre>
+ * 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
+ * </pre></p>
+ *
+ * <p>With the previous definition, the classical representation of multistep methods
+ * uses first derivatives only, i.e. it handles y<sub>n</sub>, s<sub>1</sub>(n) and
+ * q<sub>n</sub> where q<sub>n</sub> is defined as:
+ * <pre>
+ * q<sub>n</sub> = [ s<sub>1</sub>(n-1) s<sub>1</sub>(n-2) ... s<sub>1</sub>(n-(k-1)) ]<sup>T</sup>
+ * </pre>
+ * (we omit the k index in the notation for clarity).</p>
+ *
+ * <p>Another possible representation uses the Nordsieck vector with
+ * higher degrees scaled derivatives all taken at the same step, i.e it handles y<sub>n</sub>,
+ * s<sub>1</sub>(n) and r<sub>n</sub>) where r<sub>n</sub> is defined as:
+ * <pre>
+ * r<sub>n</sub> = [ s<sub>2</sub>(n), s<sub>3</sub>(n) ... s<sub>k</sub>(n) ]<sup>T</sup>
+ * </pre>
+ * (here again we omit the k index in the notation for clarity)
+ * </p>
+ *
+ * <p>Taylor series formulas show that for any index offset i, s<sub>1</sub>(n-i) can be
+ * computed from s<sub>1</sub>(n), s<sub>2</sub>(n) ... s<sub>k</sub>(n), the formula being exact
+ * for degree k polynomials.
+ * <pre>
+ * s<sub>1</sub>(n-i) = s<sub>1</sub>(n) + ∑<sub>j>0</sub> (j+1) (-i)<sup>j</sup> s<sub>j+1</sub>(n)
+ * </pre>
+ * The previous formula can be used with several values for i to compute the transform between
+ * classical representation and Nordsieck vector at step end. The transform between r<sub>n</sub>
+ * and q<sub>n</sub> resulting from the Taylor series formulas above is:
+ * <pre>
+ * q<sub>n</sub> = s<sub>1</sub>(n) u + P r<sub>n</sub>
+ * </pre>
+ * where u is the [ 1 1 ... 1 ]<sup>T</sup> vector and P is the (k-1)×(k-1) matrix built
+ * with the (j+1) (-i)<sup>j</sup> terms with i being the row number starting from 1 and j being
+ * the column number starting from 1:
+ * <pre>
+ * [ -2 3 -4 5 ... ]
+ * [ -4 12 -32 80 ... ]
+ * P = [ -6 27 -108 405 ... ]
+ * [ -8 48 -256 1280 ... ]
+ * [ ... ]
+ * </pre></p>
+ *
+ * <p>Changing -i into +i in the formula above can be used to compute a similar transform between
+ * classical representation and Nordsieck vector at step start. The resulting matrix is simply
+ * the absolute value of matrix P.</p>
+ *
+ * <p>For {@link AdamsBashforthIntegrator Adams-Bashforth} method, the Nordsieck vector
+ * at step n+1 is computed from the Nordsieck vector at step n as follows:
+ * <ul>
+ * <li>y<sub>n+1</sub> = y<sub>n</sub> + s<sub>1</sub>(n) + u<sup>T</sup> r<sub>n</sub></li>
+ * <li>s<sub>1</sub>(n+1) = h f(t<sub>n+1</sub>, y<sub>n+1</sub>)</li>
+ * <li>r<sub>n+1</sub> = (s<sub>1</sub>(n) - s<sub>1</sub>(n+1)) P<sup>-1</sup> u + P<sup>-1</sup> A P r<sub>n</sub></li>
+ * </ul>
+ * where A is a rows shifting matrix (the lower left part is an identity matrix):
+ * <pre>
+ * [ 0 0 ... 0 0 | 0 ]
+ * [ ---------------+---]
+ * [ 1 0 ... 0 0 | 0 ]
+ * A = [ 0 1 ... 0 0 | 0 ]
+ * [ ... | 0 ]
+ * [ 0 0 ... 1 0 | 0 ]
+ * [ 0 0 ... 0 1 | 0 ]
+ * </pre></p>
+ *
+ * <p>For {@link AdamsMoultonIntegrator Adams-Moulton} method, the predicted Nordsieck vector
+ * at step n+1 is computed from the Nordsieck vector at step n as follows:
+ * <ul>
+ * <li>Y<sub>n+1</sub> = y<sub>n</sub> + s<sub>1</sub>(n) + u<sup>T</sup> r<sub>n</sub></li>
+ * <li>S<sub>1</sub>(n+1) = h f(t<sub>n+1</sub>, Y<sub>n+1</sub>)</li>
+ * <li>R<sub>n+1</sub> = (s<sub>1</sub>(n) - s<sub>1</sub>(n+1)) P<sup>-1</sup> u + P<sup>-1</sup> A P r<sub>n</sub></li>
+ * </ul>
+ * From this predicted vector, the corrected vector is computed as follows:
+ * <ul>
+ * <li>y<sub>n+1</sub> = y<sub>n</sub> + S<sub>1</sub>(n+1) + [ -1 +1 -1 +1 ... ±1 ] r<sub>n+1</sub></li>
+ * <li>s<sub>1</sub>(n+1) = h f(t<sub>n+1</sub>, y<sub>n+1</sub>)</li>
+ * <li>r<sub>n+1</sub> = R<sub>n+1</sub> + (s<sub>1</sub>(n+1) - S<sub>1</sub>(n+1)) P<sup>-1</sup> u</li>
+ * </ul>
+ * where the upper case Y<sub>n+1</sub>, S<sub>1</sub>(n+1) and R<sub>n+1</sub> represent the
+ * predicted states whereas the lower case y<sub>n+1</sub>, s<sub>n+1</sub> and r<sub>n+1</sub>
+ * represent the corrected states.</p>
+ *
+ * <p>We observe that both methods use similar update formulas. In both cases a P<sup>-1</sup>u
+ * vector and a P<sup>-1</sup> A P matrix are used that do not depend on the state,
+ * they only depend on k. This class handles these transformations.</p>
+ *
+ * @param <T> the type of the field elements
+ * @since 3.6
+ */
+public class AdamsNordsieckFieldTransformer<T extends RealFieldElement<T>> {
+
+ /** Cache for already computed coefficients. */
+ private static final Map<Integer,
+ Map<Field<? extends RealFieldElement<?>>,
+ AdamsNordsieckFieldTransformer<? extends RealFieldElement<?>>>> CACHE =
+ new HashMap<Integer, Map<Field<? extends RealFieldElement<?>>,
+ AdamsNordsieckFieldTransformer<? extends RealFieldElement<?>>>>();
+
+ /** Field to which the time and state vector elements belong. */
+ private final Field<T> field;
+
+ /** Update matrix for the higher order derivatives h<sup>2</sup>/2 y'', h<sup>3</sup>/6 y''' ... */
+ private final Array2DRowFieldMatrix<T> update;
+
+ /** Update coefficients of the higher order derivatives wrt y'. */
+ private final T[] c1;
+
+ /** Simple constructor.
+ * @param field field to which the time and state vector elements belong
+ * @param n number of steps of the multistep method
+ * (excluding the one being computed)
+ */
+ private AdamsNordsieckFieldTransformer(final Field<T> field, final int n) {
+
+ this.field = field;
+ final int rows = n - 1;
+
+ // compute coefficients
+ FieldMatrix<T> bigP = buildP(rows);
+ FieldDecompositionSolver<T> pSolver =
+ new FieldLUDecomposition<T>(bigP).getSolver();
+
+ T[] u = MathArrays.buildArray(field, rows);
+ Arrays.fill(u, field.getOne());
+ c1 = pSolver.solve(new ArrayFieldVector<T>(u, false)).toArray();
+
+ // update coefficients are computed by combining transform from
+ // Nordsieck to multistep, then shifting rows to represent step advance
+ // then applying inverse transform
+ T[][] shiftedP = bigP.getData();
+ for (int i = shiftedP.length - 1; i > 0; --i) {
+ // shift rows
+ shiftedP[i] = shiftedP[i - 1];
+ }
+ shiftedP[0] = MathArrays.buildArray(field, rows);
+ Arrays.fill(shiftedP[0], field.getZero());
+ update = new Array2DRowFieldMatrix<T>(pSolver.solve(new Array2DRowFieldMatrix<T>(shiftedP, false)).getData());
+
+ }
+
+ /** Get the Nordsieck transformer for a given field and number of steps.
+ * @param field field to which the time and state vector elements belong
+ * @param nSteps number of steps of the multistep method
+ * (excluding the one being computed)
+ * @return Nordsieck transformer for the specified field and number of steps
+ * @param <T> the type of the field elements
+ */
+ public static <T extends RealFieldElement<T>> AdamsNordsieckFieldTransformer<T>
+ getInstance(final Field<T> field, final int nSteps) {
+ synchronized(CACHE) {
+ Map<Field<? extends RealFieldElement<?>>,
+ AdamsNordsieckFieldTransformer<? extends RealFieldElement<?>>> map = CACHE.get(nSteps);
+ if (map == null) {
+ map = new HashMap<Field<? extends RealFieldElement<?>>,
+ AdamsNordsieckFieldTransformer<? extends RealFieldElement<?>>>();
+ CACHE.put(nSteps, map);
+ }
+ @SuppressWarnings("unchecked")
+ AdamsNordsieckFieldTransformer<T> t = (AdamsNordsieckFieldTransformer<T>) map.get(field);
+ if (t == null) {
+ t = new AdamsNordsieckFieldTransformer<T>(field, nSteps);
+ map.put(field, t);
+ }
+ return t;
+
+ }
+ }
+
+ /** Build the P matrix.
+ * <p>The P matrix general terms are shifted (j+1) (-i)<sup>j</sup> terms
+ * with i being the row number starting from 1 and j being the column
+ * number starting from 1:
+ * <pre>
+ * [ -2 3 -4 5 ... ]
+ * [ -4 12 -32 80 ... ]
+ * P = [ -6 27 -108 405 ... ]
+ * [ -8 48 -256 1280 ... ]
+ * [ ... ]
+ * </pre></p>
+ * @param rows number of rows of the matrix
+ * @return P matrix
+ */
+ private FieldMatrix<T> buildP(final int rows) {
+
+ final T[][] pData = MathArrays.buildArray(field, rows, rows);
+
+ for (int i = 1; i <= pData.length; ++i) {
+ // build the P matrix elements from Taylor series formulas
+ final T[] pI = pData[i - 1];
+ final int factor = -i;
+ T aj = field.getZero().add(factor);
+ for (int j = 1; j <= pI.length; ++j) {
+ pI[j - 1] = aj.multiply(j + 1);
+ aj = aj.multiply(factor);
+ }
+ }
+
+ return new Array2DRowFieldMatrix<T>(pData, false);
+
+ }
+
+ /** 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 start of 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>)
+ */
+
+ public Array2DRowFieldMatrix<T> initializeHighOrderDerivatives(final T h, final T[] t,
+ final T[][] y,
+ final T[][] yDot) {
+
+ // using Taylor series with di = ti - t0, we get:
+ // y(ti) - y(t0) - di y'(t0) = di^2 / h^2 s2 + ... + di^k / h^k sk + O(h^k)
+ // y'(ti) - y'(t0) = 2 di / h^2 s2 + ... + k di^(k-1) / h^k sk + O(h^(k-1))
+ // we write these relations for i = 1 to i= 1+n/2 as a set of n + 2 linear
+ // equations depending on the Nordsieck vector [s2 ... sk rk], so s2 to sk correspond
+ // to the appropriately truncated Taylor expansion, and rk is the Taylor remainder.
+ // The goal is to have s2 to sk as accurate as possible considering the fact the sum is
+ // truncated and we don't want the error terms to be included in s2 ... sk, so we need
+ // to solve also for the remainder
+ final T[][] a = MathArrays.buildArray(field, c1.length + 1, c1.length + 1);
+ final T[][] b = MathArrays.buildArray(field, c1.length + 1, y[0].length);
+ final T[] y0 = y[0];
+ final T[] yDot0 = yDot[0];
+ for (int i = 1; i < y.length; ++i) {
+
+ final T di = t[i].subtract(t[0]);
+ final T ratio = di.divide(h);
+ T dikM1Ohk = h.reciprocal();
+
+ // linear coefficients of equations
+ // y(ti) - y(t0) - di y'(t0) and y'(ti) - y'(t0)
+ final T[] aI = a[2 * i - 2];
+ final T[] aDotI = (2 * i - 1) < a.length ? a[2 * i - 1] : null;
+ for (int j = 0; j < aI.length; ++j) {
+ dikM1Ohk = dikM1Ohk.multiply(ratio);
+ aI[j] = di.multiply(dikM1Ohk);
+ if (aDotI != null) {
+ aDotI[j] = dikM1Ohk.multiply(j + 2);
+ }
+ }
+
+ // expected value of the previous equations
+ final T[] yI = y[i];
+ final T[] yDotI = yDot[i];
+ final T[] bI = b[2 * i - 2];
+ final T[] bDotI = (2 * i - 1) < b.length ? b[2 * i - 1] : null;
+ for (int j = 0; j < yI.length; ++j) {
+ bI[j] = yI[j].subtract(y0[j]).subtract(di.multiply(yDot0[j]));
+ if (bDotI != null) {
+ bDotI[j] = yDotI[j].subtract(yDot0[j]);
+ }
+ }
+
+ }
+
+ // solve the linear system to get the best estimate of the Nordsieck vector [s2 ... sk],
+ // with the additional terms s(k+1) and c grabbing the parts after the truncated Taylor expansion
+ final FieldLUDecomposition<T> decomposition = new FieldLUDecomposition<T>(new Array2DRowFieldMatrix<T>(a, false));
+ final FieldMatrix<T> x = decomposition.getSolver().solve(new Array2DRowFieldMatrix<T>(b, false));
+
+ // extract just the Nordsieck vector [s2 ... sk]
+ final Array2DRowFieldMatrix<T> truncatedX =
+ new Array2DRowFieldMatrix<T>(field, x.getRowDimension() - 1, x.getColumnDimension());
+ for (int i = 0; i < truncatedX.getRowDimension(); ++i) {
+ for (int j = 0; j < truncatedX.getColumnDimension(); ++j) {
+ truncatedX.setEntry(i, j, x.getEntry(i, j));
+ }
+ }
+ return truncatedX;
+
+ }
+
+ /** Update the high order scaled derivatives for Adams integrators (phase 1).
+ * <p>The complete update of high order derivatives has a form similar to:
+ * <pre>
+ * r<sub>n+1</sub> = (s<sub>1</sub>(n) - s<sub>1</sub>(n+1)) P<sup>-1</sup> u + P<sup>-1</sup> A P r<sub>n</sub>
+ * </pre>
+ * 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(double[], double[], Array2DRowFieldMatrix)
+ */
+ public Array2DRowFieldMatrix<T> updateHighOrderDerivativesPhase1(final Array2DRowFieldMatrix<T> highOrder) {
+ return update.multiply(highOrder);
+ }
+
+ /** Update the high order scaled derivatives Adams integrators (phase 2).
+ * <p>The complete update of high order derivatives has a form similar to:
+ * <pre>
+ * r<sub>n+1</sub> = (s<sub>1</sub>(n) - s<sub>1</sub>(n+1)) P<sup>-1</sup> u + P<sup>-1</sup> A P r<sub>n</sub>
+ * </pre>
+ * 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) {
+ final T[][] data = highOrder.getDataRef();
+ for (int i = 0; i < data.length; ++i) {
+ final T[] dataI = data[i];
+ final T c1I = c1[i];
+ for (int j = 0; j < dataI.length; ++j) {
+ dataI[j] = dataI[j].add(c1I.multiply(start[j].subtract(end[j])));
+ }
+ }
+ }
+
+}
http://git-wip-us.apache.org/repos/asf/commons-math/blob/ba6bd8cb/src/test/java/org/apache/commons/math3/ode/nonstiff/AbstractAdamsFieldIntegratorTest.java
----------------------------------------------------------------------
diff --git a/src/test/java/org/apache/commons/math3/ode/nonstiff/AbstractAdamsFieldIntegratorTest.java b/src/test/java/org/apache/commons/math3/ode/nonstiff/AbstractAdamsFieldIntegratorTest.java
new file mode 100644
index 0000000..63ab60d
--- /dev/null
+++ b/src/test/java/org/apache/commons/math3/ode/nonstiff/AbstractAdamsFieldIntegratorTest.java
@@ -0,0 +1,261 @@
+/*
+ * 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.nonstiff;
+
+
+import org.apache.commons.math3.Field;
+import org.apache.commons.math3.RealFieldElement;
+import org.apache.commons.math3.exception.MathIllegalStateException;
+import org.apache.commons.math3.exception.MaxCountExceededException;
+import org.apache.commons.math3.exception.NumberIsTooSmallException;
+import org.apache.commons.math3.ode.AbstractFieldIntegrator;
+import org.apache.commons.math3.ode.FieldExpandableODE;
+import org.apache.commons.math3.ode.FieldODEState;
+import org.apache.commons.math3.ode.FieldODEStateAndDerivative;
+import org.apache.commons.math3.ode.FirstOrderFieldIntegrator;
+import org.apache.commons.math3.ode.MultistepFieldIntegrator;
+import org.apache.commons.math3.ode.TestFieldProblem1;
+import org.apache.commons.math3.ode.TestFieldProblem5;
+import org.apache.commons.math3.ode.TestFieldProblem6;
+import org.apache.commons.math3.ode.TestFieldProblemAbstract;
+import org.apache.commons.math3.ode.TestFieldProblemHandler;
+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.junit.Assert;
+import org.junit.Test;
+
+public abstract class AbstractAdamsFieldIntegratorTest {
+
+ protected abstract <T extends RealFieldElement<T>> AdamsFieldIntegrator<T>
+ createIntegrator(Field<T> field, final int nSteps, final double minStep, final double maxStep,
+ final double scalAbsoluteTolerance, final double scalRelativeTolerance);
+
+ protected abstract <T extends RealFieldElement<T>> AdamsFieldIntegrator<T>
+ createIntegrator(Field<T> field, final int nSteps, final double minStep, final double maxStep,
+ final double[] vecAbsoluteTolerance, final double[] vecRelativeTolerance);
+
+ @Test(expected=NumberIsTooSmallException.class)
+ public abstract void testMinStep();
+
+ protected <T extends RealFieldElement<T>> void doDimensionCheck(final Field<T> field) {
+ TestFieldProblem1<T> pb = new TestFieldProblem1<T>(field);
+
+ double minStep = pb.getFinalTime().subtract(pb.getInitialState().getTime()).multiply(0.1).getReal();
+ double maxStep = pb.getFinalTime().subtract(pb.getInitialState().getTime()).getReal();
+ double[] vecAbsoluteTolerance = { 1.0e-15, 1.0e-16 };
+ double[] vecRelativeTolerance = { 1.0e-15, 1.0e-16 };
+
+ FirstOrderFieldIntegrator<T> integ = createIntegrator(field, 4, minStep, maxStep,
+ vecAbsoluteTolerance,
+ vecRelativeTolerance);
+ TestFieldProblemHandler<T> handler = new TestFieldProblemHandler<T>(pb, integ);
+ integ.addStepHandler(handler);
+ integ.integrate(new FieldExpandableODE<T>(pb), pb.getInitialState(), pb.getFinalTime());
+
+ }
+
+ @Test
+ public abstract void testIncreasingTolerance();
+
+ protected <T extends RealFieldElement<T>> void doTestIncreasingTolerance(final Field<T> field,
+ int ratioMin, int ratioMax) {
+
+ int previousCalls = Integer.MAX_VALUE;
+ for (int i = -12; i < -5; ++i) {
+ TestFieldProblem1<T> pb = new TestFieldProblem1<T>(field);
+ double minStep = 0;
+ double maxStep = pb.getFinalTime().subtract(pb.getInitialState().getTime()).getReal();
+ double scalAbsoluteTolerance = FastMath.pow(10.0, i);
+ double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
+
+ FirstOrderFieldIntegrator<T> integ = createIntegrator(field, 4, minStep, maxStep,
+ scalAbsoluteTolerance,
+ scalRelativeTolerance);
+ TestFieldProblemHandler<T> handler = new TestFieldProblemHandler<T>(pb, integ);
+ integ.addStepHandler(handler);
+ integ.integrate(new FieldExpandableODE<T>(pb), pb.getInitialState(), pb.getFinalTime());
+
+ Assert.assertTrue(handler.getMaximalValueError().getReal() > ratioMin * scalAbsoluteTolerance);
+ Assert.assertTrue(handler.getMaximalValueError().getReal() < ratioMax * scalAbsoluteTolerance);
+
+ int calls = pb.getCalls();
+ Assert.assertEquals(integ.getEvaluations(), calls);
+ Assert.assertTrue(calls <= previousCalls);
+ previousCalls = calls;
+
+ }
+
+ }
+
+ @Test(expected = MaxCountExceededException.class)
+ public abstract void exceedMaxEvaluations();
+
+ protected <T extends RealFieldElement<T>> void doExceedMaxEvaluations(final Field<T> field) {
+
+ TestFieldProblem1<T> pb = new TestFieldProblem1<T>(field);
+ double range = pb.getFinalTime().subtract(pb.getInitialState().getTime()).getReal();
+
+ FirstOrderFieldIntegrator<T> integ = createIntegrator(field, 2, 0, range, 1.0e-12, 1.0e-12);
+ TestFieldProblemHandler<T> handler = new TestFieldProblemHandler<T>(pb, integ);
+ integ.addStepHandler(handler);
+ integ.setMaxEvaluations(650);
+ integ.integrate(new FieldExpandableODE<T>(pb), pb.getInitialState(), pb.getFinalTime());
+
+ }
+
+ @Test
+ public abstract void backward();
+
+ protected <T extends RealFieldElement<T>> void doBackward(final Field<T> field,
+ final double epsilonLast,
+ final double epsilonMaxValue,
+ final double epsilonMaxTime,
+ final String name) {
+
+ TestFieldProblem5<T> pb = new TestFieldProblem5<T>(field);
+ double range = pb.getFinalTime().subtract(pb.getInitialState().getTime()).getReal();
+
+ AdamsFieldIntegrator<T> integ = createIntegrator(field, 4, 0, range, 1.0e-12, 1.0e-12);
+ integ.setStarterIntegrator(new PerfectStarter<T>(pb, (integ.getNSteps() + 5) / 2));
+ TestFieldProblemHandler<T> handler = new TestFieldProblemHandler<T>(pb, integ);
+ integ.addStepHandler(handler);
+ integ.integrate(new FieldExpandableODE<T>(pb), pb.getInitialState(), pb.getFinalTime());
+
+ Assert.assertEquals(0.0, handler.getLastError().getReal(), epsilonLast);
+ Assert.assertEquals(0.0, handler.getMaximalValueError().getReal(), epsilonMaxValue);
+ Assert.assertEquals(0, handler.getMaximalTimeError().getReal(), epsilonMaxTime);
+ Assert.assertEquals(name, integ.getName());
+ }
+
+ @Test
+ public abstract void polynomial();
+
+ protected <T extends RealFieldElement<T>> void doPolynomial(final Field<T> field,
+ final int nLimit,
+ final double epsilonBad,
+ final double epsilonGood) {
+ TestFieldProblem6<T> pb = new TestFieldProblem6<T>(field);
+ double range = pb.getFinalTime().subtract(pb.getInitialState().getTime()).abs().getReal();
+
+ for (int nSteps = 2; nSteps < 8; ++nSteps) {
+ AdamsFieldIntegrator<T> integ = createIntegrator(field, nSteps, 1.0e-6 * range, 0.1 * range, 1.0e-4, 1.0e-4);
+ integ.setStarterIntegrator(new PerfectStarter<T>(pb, nSteps));
+ TestFieldProblemHandler<T> handler = new TestFieldProblemHandler<T>(pb, integ);
+ integ.addStepHandler(handler);
+ integ.integrate(new FieldExpandableODE<T>(pb), pb.getInitialState(), pb.getFinalTime());
+ if (nSteps < nLimit) {
+ Assert.assertTrue(handler.getMaximalValueError().getReal() > epsilonBad);
+ } else {
+ Assert.assertTrue(handler.getMaximalValueError().getReal() < epsilonGood);
+ }
+ }
+
+ }
+
+ @Test(expected=MathIllegalStateException.class)
+ public abstract void testStartFailure();
+
+ protected <T extends RealFieldElement<T>> void doTestStartFailure(final Field<T> field) {
+ TestFieldProblem1<T> pb = new TestFieldProblem1<T>(field);
+ double minStep = pb.getFinalTime().subtract(pb.getInitialState().getTime()).multiply(0.0001).getReal();
+ double maxStep = pb.getFinalTime().subtract(pb.getInitialState().getTime()).getReal();
+ double scalAbsoluteTolerance = 1.0e-6;
+ double scalRelativeTolerance = 1.0e-7;
+
+ MultistepFieldIntegrator<T> integ = createIntegrator(field, 6, minStep, maxStep,
+ scalAbsoluteTolerance,
+ scalRelativeTolerance);
+ integ.setStarterIntegrator(new DormandPrince853FieldIntegrator<T>(field, maxStep * 0.5, maxStep, 0.1, 0.1));
+ TestFieldProblemHandler<T> handler = new TestFieldProblemHandler<T>(pb, integ);
+ integ.addStepHandler(handler);
+ integ.integrate(new FieldExpandableODE<T>(pb), pb.getInitialState(), pb.getFinalTime());
+
+ }
+
+ private static class PerfectStarter<T extends RealFieldElement<T>> extends AbstractFieldIntegrator<T> {
+
+ private final PerfectInterpolator<T> interpolator;
+ private final int nbSteps;
+
+ public PerfectStarter(final TestFieldProblemAbstract<T> problem, final int nbSteps) {
+ super(problem.getField(), "perfect-starter");
+ this.interpolator = new PerfectInterpolator<T>(problem);
+ this.nbSteps = nbSteps;
+ }
+
+ public FieldODEStateAndDerivative<T> integrate(FieldExpandableODE<T> equations,
+ FieldODEState<T> initialState, T finalTime) {
+ T tStart = initialState.getTime().add(finalTime.subtract(initialState.getTime()).multiply(0.01));
+ getEvaluationsCounter().increment(nbSteps);
+ interpolator.setCurrentTime(initialState.getTime());
+ for (int i = 0; i < nbSteps; ++i) {
+ T tK = initialState.getTime().multiply(nbSteps - 1 - (i + 1)).add(tStart.multiply(i + 1)).divide(nbSteps - 1);
+ interpolator.setPreviousTime(interpolator.getCurrentTime());
+ interpolator.setCurrentTime(tK);
+ for (FieldStepHandler<T> handler : getStepHandlers()) {
+ handler.handleStep(interpolator, i == nbSteps - 1);
+ }
+ }
+ return interpolator.getInterpolatedState(tStart);
+ }
+
+ }
+
+ private static class PerfectInterpolator<T extends RealFieldElement<T>> implements FieldStepInterpolator<T> {
+ private final TestFieldProblemAbstract<T> problem;
+ private T previousTime;
+ private T currentTime;
+
+ public PerfectInterpolator(final TestFieldProblemAbstract<T> problem) {
+ this.problem = problem;
+ }
+
+ public void setPreviousTime(T previousTime) {
+ this.previousTime = previousTime;
+ }
+
+ public void setCurrentTime(T currentTime) {
+ this.currentTime = currentTime;
+ }
+
+ public T getCurrentTime() {
+ return currentTime;
+ }
+
+ public boolean isForward() {
+ return problem.getFinalTime().subtract(problem.getInitialState().getTime()).getReal() >= 0;
+ }
+
+ public FieldODEStateAndDerivative<T> getPreviousState() {
+ return getInterpolatedState(previousTime);
+ }
+
+ public FieldODEStateAndDerivative<T> getCurrentState() {
+ return getInterpolatedState(currentTime);
+ }
+
+ public FieldODEStateAndDerivative<T> getInterpolatedState(T time) {
+ T[] y = problem.computeTheoreticalState(time);
+ T[] yDot = problem.computeDerivatives(time, y);
+ return new FieldODEStateAndDerivative<T>(time, y, yDot);
+ }
+
+ }
+
+}