simplified and reorganized slightly the Adams integrators class hierarchy
this will allow a future BDF integrator development for stiff problems git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@789159 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe
parent
7614049449
commit
3cb90658a4
|
|
@ -17,13 +17,8 @@
|
|||
|
||||
package org.apache.commons.math.ode;
|
||||
|
||||
|
||||
import java.io.IOException;
|
||||
import java.io.ObjectInputStream;
|
||||
import java.io.ObjectOutputStream;
|
||||
import java.lang.reflect.Field;
|
||||
|
||||
import org.apache.commons.math.MathRuntimeException;
|
||||
import org.apache.commons.math.linear.Array2DRowRealMatrix;
|
||||
import org.apache.commons.math.linear.RealMatrix;
|
||||
import org.apache.commons.math.ode.nonstiff.AdaptiveStepsizeIntegrator;
|
||||
import org.apache.commons.math.ode.nonstiff.DormandPrince853Integrator;
|
||||
|
|
@ -41,9 +36,6 @@ import org.apache.commons.math.ode.sampling.StepInterpolator;
|
|||
*/
|
||||
public abstract class MultistepIntegrator extends AdaptiveStepsizeIntegrator {
|
||||
|
||||
/** Transformer. */
|
||||
protected final transient NordsieckTransformer transformer;
|
||||
|
||||
/** Starter integrator. */
|
||||
private FirstOrderIntegrator starter;
|
||||
|
||||
|
|
@ -56,7 +48,7 @@ public abstract class MultistepIntegrator extends AdaptiveStepsizeIntegrator {
|
|||
/** 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(k))</p>
|
||||
*/
|
||||
protected RealMatrix nordsieck;
|
||||
protected Array2DRowRealMatrix nordsieck;
|
||||
|
||||
/** Stepsize control exponent. */
|
||||
private double exp;
|
||||
|
|
@ -104,7 +96,6 @@ public abstract class MultistepIntegrator extends AdaptiveStepsizeIntegrator {
|
|||
scalAbsoluteTolerance,
|
||||
scalRelativeTolerance);
|
||||
this.nSteps = nSteps;
|
||||
transformer = NordsieckTransformer.getInstance(nSteps + 1);
|
||||
|
||||
exp = -1.0 / order;
|
||||
|
||||
|
|
@ -141,7 +132,6 @@ public abstract class MultistepIntegrator extends AdaptiveStepsizeIntegrator {
|
|||
vecAbsoluteTolerance,
|
||||
vecRelativeTolerance);
|
||||
this.nSteps = nSteps;
|
||||
transformer = NordsieckTransformer.getInstance(nSteps + 1);
|
||||
|
||||
exp = -1.0 / order;
|
||||
|
||||
|
|
@ -216,6 +206,15 @@ public abstract class MultistepIntegrator extends AdaptiveStepsizeIntegrator {
|
|||
|
||||
}
|
||||
|
||||
/** Initialize the high order scaled derivatives at step start.
|
||||
* @param first first scaled derivative at step start
|
||||
* @param multistep scaled derivatives after step start (hy'1, ..., hy'k-1)
|
||||
* will be modified
|
||||
* @return high order derivatives at step start
|
||||
*/
|
||||
protected abstract Array2DRowRealMatrix initializeHighOrderDerivatives(final double[] first,
|
||||
final double[][] multistep);
|
||||
|
||||
/** Get the minimal reduction factor for stepsize control.
|
||||
* @return minimal reduction factor
|
||||
*/
|
||||
|
|
@ -266,43 +265,15 @@ public abstract class MultistepIntegrator extends AdaptiveStepsizeIntegrator {
|
|||
return Math.min(maxGrowth, Math.max(minReduction, safety * Math.pow(error, exp)));
|
||||
}
|
||||
|
||||
/** Serialize the instance.
|
||||
* @param oos stream where object should be written
|
||||
* @throws IOException if object cannot be written to stream
|
||||
*/
|
||||
private void writeObject(ObjectOutputStream oos)
|
||||
throws IOException {
|
||||
oos.defaultWriteObject();
|
||||
oos.writeInt(nSteps);
|
||||
}
|
||||
|
||||
/** Deserialize the instance.
|
||||
* @param ois stream from which the object should be read
|
||||
* @throws ClassNotFoundException if a class in the stream cannot be found
|
||||
* @throws IOException if object cannot be read from the stream
|
||||
*/
|
||||
private void readObject(ObjectInputStream ois)
|
||||
throws ClassNotFoundException, IOException {
|
||||
try {
|
||||
|
||||
ois.defaultReadObject();
|
||||
final int nSteps = ois.readInt();
|
||||
|
||||
final Class<MultistepIntegrator> cl = MultistepIntegrator.class;
|
||||
final Field f = cl.getDeclaredField("transformer");
|
||||
f.setAccessible(true);
|
||||
f.set(this, NordsieckTransformer.getInstance(nSteps + 1));
|
||||
|
||||
} catch (NoSuchFieldException nsfe) {
|
||||
IOException ioe = new IOException();
|
||||
ioe.initCause(nsfe);
|
||||
throw ioe;
|
||||
} catch (IllegalAccessException iae) {
|
||||
IOException ioe = new IOException();
|
||||
ioe.initCause(iae);
|
||||
throw ioe;
|
||||
}
|
||||
|
||||
/** Transformer used to convert the first step to Nordsieck representation. */
|
||||
public static interface NordsieckTransformer {
|
||||
/** Initialize the high order scaled derivatives at step start.
|
||||
* @param first first scaled derivative at step start
|
||||
* @param multistep scaled derivatives after step start (hy'1, ..., hy'k-1)
|
||||
* will be modified
|
||||
* @return high order derivatives at step start
|
||||
*/
|
||||
RealMatrix initializeHighOrderDerivatives(double[] first, double[][] multistep);
|
||||
}
|
||||
|
||||
/** Specialized step handler storing the first step. */
|
||||
|
|
@ -344,7 +315,7 @@ public abstract class MultistepIntegrator extends AdaptiveStepsizeIntegrator {
|
|||
}
|
||||
multistep[i - 1] = msI;
|
||||
}
|
||||
nordsieck = transformer.initializeHighOrderDerivatives(scaled, multistep);
|
||||
nordsieck = initializeHighOrderDerivatives(scaled, multistep);
|
||||
|
||||
// stop the integrator after the first step has been handled
|
||||
throw new InitializationCompletedMarkerException();
|
||||
|
|
|
|||
|
|
@ -17,11 +17,10 @@
|
|||
|
||||
package org.apache.commons.math.ode.nonstiff;
|
||||
|
||||
import org.apache.commons.math.linear.RealMatrix;
|
||||
import org.apache.commons.math.linear.Array2DRowRealMatrix;
|
||||
import org.apache.commons.math.ode.DerivativeException;
|
||||
import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
|
||||
import org.apache.commons.math.ode.IntegratorException;
|
||||
import org.apache.commons.math.ode.MultistepIntegrator;
|
||||
import org.apache.commons.math.ode.events.CombinedEventsManager;
|
||||
import org.apache.commons.math.ode.sampling.NordsieckStepInterpolator;
|
||||
import org.apache.commons.math.ode.sampling.StepHandler;
|
||||
|
|
@ -139,10 +138,10 @@ import org.apache.commons.math.ode.sampling.StepHandler;
|
|||
* @version $Revision$ $Date$
|
||||
* @since 2.0
|
||||
*/
|
||||
public class AdamsBashforthIntegrator extends MultistepIntegrator {
|
||||
public class AdamsBashforthIntegrator extends AdamsIntegrator {
|
||||
|
||||
/**
|
||||
* Build an Adams-Bashforth with the given order and step size.
|
||||
* Build an Adams-Bashforth integrator with the given order and step control parameters.
|
||||
* @param nSteps number of steps of the method excluding the one being computed
|
||||
* @param minStep minimal step (must be positive even for backward
|
||||
* integration), the last step can be smaller than this
|
||||
|
|
@ -162,7 +161,7 @@ public class AdamsBashforthIntegrator extends MultistepIntegrator {
|
|||
}
|
||||
|
||||
/**
|
||||
* Build an Adams-Bashforth with the given order and step size.
|
||||
* Build an Adams-Bashforth integrator with the given order and step control parameters.
|
||||
* @param nSteps number of steps of the method excluding the one being computed
|
||||
* @param minStep minimal step (must be positive even for backward
|
||||
* integration), the last step can be smaller than this
|
||||
|
|
@ -261,9 +260,8 @@ public class AdamsBashforthIntegrator extends MultistepIntegrator {
|
|||
for (int j = 0; j < y0.length; ++j) {
|
||||
predictedScaled[j] = stepSize * yDot[j];
|
||||
}
|
||||
final RealMatrix nordsieckTmp =
|
||||
transformer.updateHighOrderDerivativesPhase1(nordsieck);
|
||||
transformer.updateHighOrderDerivativesPhase2(scaled, predictedScaled, nordsieckTmp);
|
||||
final Array2DRowRealMatrix nordsieckTmp = updateHighOrderDerivativesPhase1(nordsieck);
|
||||
updateHighOrderDerivativesPhase2(scaled, predictedScaled, nordsieckTmp);
|
||||
|
||||
// discrete events handling
|
||||
interpolatorTmp.reinitialize(stepEnd, stepSize, predictedScaled, nordsieckTmp);
|
||||
|
|
|
|||
|
|
@ -0,0 +1,132 @@
|
|||
/*
|
||||
* 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.math.ode.nonstiff;
|
||||
|
||||
import org.apache.commons.math.linear.Array2DRowRealMatrix;
|
||||
import org.apache.commons.math.linear.RealMatrix;
|
||||
import org.apache.commons.math.ode.DerivativeException;
|
||||
import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
|
||||
import org.apache.commons.math.ode.IntegratorException;
|
||||
import org.apache.commons.math.ode.MultistepIntegrator;
|
||||
|
||||
|
||||
/** Base class for {@link AdamsBashforthIntegrator Adams-Bashforth} and
|
||||
* {@link AdamsMoultonIntegrator Adams-Moulton} integrators.
|
||||
* @version $Revision$ $Date$
|
||||
* @since 2.0
|
||||
*/
|
||||
public abstract class AdamsIntegrator extends MultistepIntegrator {
|
||||
|
||||
/** Transformer. */
|
||||
private final AdamsNordsieckTransformer transformer;
|
||||
|
||||
/**
|
||||
* Build an Adams integrator with the given order and step control prameters.
|
||||
* @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 (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 IllegalArgumentException if order is 1 or less
|
||||
*/
|
||||
public AdamsIntegrator(final String name, final int nSteps, final int order,
|
||||
final double minStep, final double maxStep,
|
||||
final double scalAbsoluteTolerance,
|
||||
final double scalRelativeTolerance)
|
||||
throws IllegalArgumentException {
|
||||
super(name, nSteps, order, minStep, maxStep,
|
||||
scalAbsoluteTolerance, scalRelativeTolerance);
|
||||
transformer = AdamsNordsieckTransformer.getInstance(nSteps);
|
||||
}
|
||||
|
||||
/**
|
||||
* Build an Adams integrator with the given order and step control parameters.
|
||||
* @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 (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
|
||||
* @exception IllegalArgumentException if order is 1 or less
|
||||
*/
|
||||
public AdamsIntegrator(final String name, final int nSteps, final int order,
|
||||
final double minStep, final double maxStep,
|
||||
final double[] vecAbsoluteTolerance,
|
||||
final double[] vecRelativeTolerance)
|
||||
throws IllegalArgumentException {
|
||||
super(name, nSteps, order, minStep, maxStep,
|
||||
vecAbsoluteTolerance, vecRelativeTolerance);
|
||||
transformer = AdamsNordsieckTransformer.getInstance(nSteps);
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
@Override
|
||||
public abstract double integrate(final FirstOrderDifferentialEquations equations,
|
||||
final double t0, final double[] y0,
|
||||
final double t, final double[] y)
|
||||
throws DerivativeException, IntegratorException;
|
||||
|
||||
/** {@inheritDoc} */
|
||||
@Override
|
||||
protected Array2DRowRealMatrix initializeHighOrderDerivatives(final double[] first,
|
||||
final double[][] multistep) {
|
||||
return transformer.initializeHighOrderDerivatives(first, multistep);
|
||||
}
|
||||
|
||||
/** 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[], RealMatrix)
|
||||
*/
|
||||
public Array2DRowRealMatrix updateHighOrderDerivativesPhase1(final Array2DRowRealMatrix 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(RealMatrix)
|
||||
*/
|
||||
public void updateHighOrderDerivativesPhase2(final double[] start,
|
||||
final double[] end,
|
||||
final Array2DRowRealMatrix highOrder) {
|
||||
transformer.updateHighOrderDerivativesPhase2(start, end, highOrder);
|
||||
}
|
||||
|
||||
}
|
||||
|
|
@ -19,13 +19,12 @@ package org.apache.commons.math.ode.nonstiff;
|
|||
|
||||
import java.util.Arrays;
|
||||
|
||||
import org.apache.commons.math.linear.Array2DRowRealMatrix;
|
||||
import org.apache.commons.math.linear.MatrixVisitorException;
|
||||
import org.apache.commons.math.linear.RealMatrix;
|
||||
import org.apache.commons.math.linear.RealMatrixPreservingVisitor;
|
||||
import org.apache.commons.math.ode.DerivativeException;
|
||||
import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
|
||||
import org.apache.commons.math.ode.IntegratorException;
|
||||
import org.apache.commons.math.ode.MultistepIntegrator;
|
||||
import org.apache.commons.math.ode.events.CombinedEventsManager;
|
||||
import org.apache.commons.math.ode.sampling.NordsieckStepInterpolator;
|
||||
import org.apache.commons.math.ode.sampling.StepHandler;
|
||||
|
|
@ -128,7 +127,7 @@ import org.apache.commons.math.ode.sampling.StepHandler;
|
|||
* <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>
|
||||
* <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>
|
||||
|
|
@ -156,7 +155,7 @@ import org.apache.commons.math.ode.sampling.StepHandler;
|
|||
* @version $Revision$ $Date$
|
||||
* @since 2.0
|
||||
*/
|
||||
public class AdamsMoultonIntegrator extends MultistepIntegrator {
|
||||
public class AdamsMoultonIntegrator extends AdamsIntegrator {
|
||||
|
||||
/**
|
||||
* Build an Adams-Moulton integrator with the given order and error control parameters.
|
||||
|
|
@ -179,7 +178,7 @@ public class AdamsMoultonIntegrator extends MultistepIntegrator {
|
|||
}
|
||||
|
||||
/**
|
||||
* Build an Adams-Moulton integrator with the given order and step size.
|
||||
* Build an Adams-Moulton integrator with the given order and error control parameters.
|
||||
* @param nSteps number of steps of the method excluding the one being computed
|
||||
* @param minStep minimal step (must be positive even for backward
|
||||
* integration), the last step can be smaller than this
|
||||
|
|
@ -264,9 +263,8 @@ public class AdamsMoultonIntegrator extends MultistepIntegrator {
|
|||
for (int j = 0; j < y0.length; ++j) {
|
||||
predictedScaled[j] = stepSize * yDot[j];
|
||||
}
|
||||
final RealMatrix nordsieckTmp =
|
||||
transformer.updateHighOrderDerivativesPhase1(nordsieck);
|
||||
transformer.updateHighOrderDerivativesPhase2(scaled, predictedScaled, nordsieckTmp);
|
||||
final Array2DRowRealMatrix nordsieckTmp = updateHighOrderDerivativesPhase1(nordsieck);
|
||||
updateHighOrderDerivativesPhase2(scaled, predictedScaled, nordsieckTmp);
|
||||
|
||||
// apply correction (C in the PECE sequence)
|
||||
error = nordsieckTmp.walkInOptimizedOrder(new Corrector(y, predictedScaled, yTmp));
|
||||
|
|
@ -281,7 +279,7 @@ public class AdamsMoultonIntegrator extends MultistepIntegrator {
|
|||
for (int j = 0; j < y0.length; ++j) {
|
||||
correctedScaled[j] = stepSize * yDot[j];
|
||||
}
|
||||
transformer.updateHighOrderDerivativesPhase2(predictedScaled, correctedScaled, nordsieckTmp);
|
||||
updateHighOrderDerivativesPhase2(predictedScaled, correctedScaled, nordsieckTmp);
|
||||
|
||||
// discrete events handling
|
||||
interpolatorTmp.reinitialize(stepEnd, stepSize, correctedScaled, nordsieckTmp);
|
||||
|
|
|
|||
Loading…
Reference in New Issue