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/dd9dc945 Tree: http://git-wip-us.apache.org/repos/asf/commons-math/tree/dd9dc945 Diff: http://git-wip-us.apache.org/repos/asf/commons-math/diff/dd9dc945 Branch: refs/heads/master Commit: dd9dc9457722787992347ae54d9a8fe8d4b91d2e Parents: 2a49849 Author: Luc Maisonobe <[email protected]> Authored: Wed Jan 6 14:18:26 2016 +0100 Committer: Luc Maisonobe <[email protected]> Committed: Wed Jan 6 14:18:26 2016 +0100 ---------------------------------------------------------------------- .../math4/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/dd9dc945/src/main/java/org/apache/commons/math4/ode/MultistepFieldIntegrator.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math4/ode/MultistepFieldIntegrator.java b/src/main/java/org/apache/commons/math4/ode/MultistepFieldIntegrator.java new file mode 100644 index 0000000..feec974 --- /dev/null +++ b/src/main/java/org/apache/commons/math4/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.math4.ode; + +import org.apache.commons.math4.Field; +import org.apache.commons.math4.RealFieldElement; +import org.apache.commons.math4.exception.DimensionMismatchException; +import org.apache.commons.math4.exception.MathIllegalStateException; +import org.apache.commons.math4.exception.MaxCountExceededException; +import org.apache.commons.math4.exception.NoBracketingException; +import org.apache.commons.math4.exception.NumberIsTooSmallException; +import org.apache.commons.math4.exception.util.LocalizedFormats; +import org.apache.commons.math4.linear.Array2DRowFieldMatrix; +import org.apache.commons.math4.ode.nonstiff.AdaptiveStepsizeFieldIntegrator; +import org.apache.commons.math4.ode.nonstiff.DormandPrince853FieldIntegrator; +import org.apache.commons.math4.ode.sampling.FieldStepHandler; +import org.apache.commons.math4.ode.sampling.FieldStepInterpolator; +import org.apache.commons.math4.util.FastMath; +import org.apache.commons.math4.util.MathArrays; +import org.apache.commons.math4.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.math4.ode.nonstiff.AdamsBashforthFieldIntegrator + * @see org.apache.commons.math4.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/dd9dc945/src/main/java/org/apache/commons/math4/ode/nonstiff/AdamsFieldIntegrator.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math4/ode/nonstiff/AdamsFieldIntegrator.java b/src/main/java/org/apache/commons/math4/ode/nonstiff/AdamsFieldIntegrator.java new file mode 100644 index 0000000..9039ca1 --- /dev/null +++ b/src/main/java/org/apache/commons/math4/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.math4.ode.nonstiff; + +import org.apache.commons.math4.Field; +import org.apache.commons.math4.RealFieldElement; +import org.apache.commons.math4.exception.DimensionMismatchException; +import org.apache.commons.math4.exception.MaxCountExceededException; +import org.apache.commons.math4.exception.NoBracketingException; +import org.apache.commons.math4.exception.NumberIsTooSmallException; +import org.apache.commons.math4.linear.Array2DRowFieldMatrix; +import org.apache.commons.math4.ode.FieldExpandableODE; +import org.apache.commons.math4.ode.FieldODEState; +import org.apache.commons.math4.ode.FieldODEStateAndDerivative; +import org.apache.commons.math4.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/dd9dc945/src/main/java/org/apache/commons/math4/ode/nonstiff/AdamsFieldStepInterpolator.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math4/ode/nonstiff/AdamsFieldStepInterpolator.java b/src/main/java/org/apache/commons/math4/ode/nonstiff/AdamsFieldStepInterpolator.java new file mode 100644 index 0000000..78c3c8e --- /dev/null +++ b/src/main/java/org/apache/commons/math4/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.math4.ode.nonstiff; + +import java.util.Arrays; + +import org.apache.commons.math4.RealFieldElement; +import org.apache.commons.math4.linear.Array2DRowFieldMatrix; +import org.apache.commons.math4.ode.FieldEquationsMapper; +import org.apache.commons.math4.ode.FieldODEStateAndDerivative; +import org.apache.commons.math4.ode.sampling.AbstractFieldStepInterpolator; +import org.apache.commons.math4.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/dd9dc945/src/main/java/org/apache/commons/math4/ode/nonstiff/AdamsNordsieckFieldTransformer.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math4/ode/nonstiff/AdamsNordsieckFieldTransformer.java b/src/main/java/org/apache/commons/math4/ode/nonstiff/AdamsNordsieckFieldTransformer.java new file mode 100644 index 0000000..93ec85e --- /dev/null +++ b/src/main/java/org/apache/commons/math4/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.math4.ode.nonstiff; + +import java.util.Arrays; +import java.util.HashMap; +import java.util.Map; + +import org.apache.commons.math4.Field; +import org.apache.commons.math4.RealFieldElement; +import org.apache.commons.math4.linear.Array2DRowFieldMatrix; +import org.apache.commons.math4.linear.ArrayFieldVector; +import org.apache.commons.math4.linear.FieldDecompositionSolver; +import org.apache.commons.math4.linear.FieldLUDecomposition; +import org.apache.commons.math4.linear.FieldMatrix; +import org.apache.commons.math4.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/dd9dc945/src/test/java/org/apache/commons/math4/ode/nonstiff/AbstractAdamsFieldIntegratorTest.java ---------------------------------------------------------------------- diff --git a/src/test/java/org/apache/commons/math4/ode/nonstiff/AbstractAdamsFieldIntegratorTest.java b/src/test/java/org/apache/commons/math4/ode/nonstiff/AbstractAdamsFieldIntegratorTest.java new file mode 100644 index 0000000..0d22831 --- /dev/null +++ b/src/test/java/org/apache/commons/math4/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.math4.ode.nonstiff; + + +import org.apache.commons.math4.Field; +import org.apache.commons.math4.RealFieldElement; +import org.apache.commons.math4.exception.MathIllegalStateException; +import org.apache.commons.math4.exception.MaxCountExceededException; +import org.apache.commons.math4.exception.NumberIsTooSmallException; +import org.apache.commons.math4.ode.AbstractFieldIntegrator; +import org.apache.commons.math4.ode.FieldExpandableODE; +import org.apache.commons.math4.ode.FieldODEState; +import org.apache.commons.math4.ode.FieldODEStateAndDerivative; +import org.apache.commons.math4.ode.FirstOrderFieldIntegrator; +import org.apache.commons.math4.ode.MultistepFieldIntegrator; +import org.apache.commons.math4.ode.TestFieldProblem1; +import org.apache.commons.math4.ode.TestFieldProblem5; +import org.apache.commons.math4.ode.TestFieldProblem6; +import org.apache.commons.math4.ode.TestFieldProblemAbstract; +import org.apache.commons.math4.ode.TestFieldProblemHandler; +import org.apache.commons.math4.ode.sampling.FieldStepHandler; +import org.apache.commons.math4.ode.sampling.FieldStepInterpolator; +import org.apache.commons.math4.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); + } + + } + +}
