Apache
/
commons-math
mirror of https://github.com/apache/commons-math.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

changed Nordsieck transformer to an Adams-specific Nordsieck transformer 

the transformer associated with BDF integrator will be quite different

git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@789158 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe 2009-06-28 21:54:33 +00:00
parent 43bc08d7c3
commit 7614049449
1 changed files with 50 additions and 61 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

111
src/java/org/apache/commons/math/ode/NordsieckTransformer.java → src/java/org/apache/commons/math/ode/nonstiff/AdamsNordsieckTransformer.java
View File

@ -15,7 +15,7 @@
* limitations under the License. * limitations under the License.
*/ */
package org.apache.commons.math.ode; package org.apache.commons.math.ode.nonstiff;
import java.util.Arrays; import java.util.Arrays;
import java.util.HashMap; import java.util.HashMap;
@ -32,13 +32,12 @@ import org.apache.commons.math.linear.FieldMatrix;
import org.apache.commons.math.linear.MatrixUtils; import org.apache.commons.math.linear.MatrixUtils;
import org.apache.commons.math.linear.MatrixVisitorException; import org.apache.commons.math.linear.MatrixVisitorException;
import org.apache.commons.math.linear.RealMatrix; import org.apache.commons.math.linear.RealMatrix;
import org.apache.commons.math.ode.nonstiff.AdamsBashforthIntegrator;
import org.apache.commons.math.ode.nonstiff.AdamsMoultonIntegrator;
/** Transformer for Nordsieck vectors. /** Transformer to Nordsieck vectors for Adams integrators.
* <p>This class i used by {@link MultistepIntegrator multistep integrators} * <p>This class i used by {@link AdamsBashforthIntegrator Adams-Bashforth} and
* to convert between classical representation with several previous first * {@link AdamsMoultonIntegrator Adams-Moulton} integrators to convert between
* derivatives and Nordsieck representation with higher order scaled derivatives.</p> * 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: * <p>We define scaled derivatives s<sub>i</sub>(n) at step n as:
* <pre> * <pre>
@ -89,18 +88,9 @@ import org.apache.commons.math.ode.nonstiff.AdamsMoultonIntegrator;
* </pre></p> * </pre></p>
* *
* <p>Changing -i into +i in the formula above can be used to compute a similar transform between * <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 Q matrix is simply * classical representation and Nordsieck vector at step start. The resulting matrix is simply
* the absolute value of matrix P.</p> * the absolute value of matrix P.</p>
* *
* <p>Using the Nordsieck vector has several advantages:
* <ul>
* <li>it greatly simplifies step interpolation as the interpolator mainly applies
* Taylor series formulas,</li>
* <li>it simplifies step changes that occur when discrete events that truncate
* the step are triggered,</li>
* <li>it allows to extend the methods in order to support adaptive stepsize (not implemented yet).</li>
* </ul></p>
*
* <p>For {@link AdamsBashforthIntegrator Adams-Bashforth} method, the Nordsieck vector * <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: * at step n+1 is computed from the Nordsieck vector at step n as follows:
* <ul> * <ul>
@ -126,16 +116,6 @@ import org.apache.commons.math.ode.nonstiff.AdamsMoultonIntegrator;
* <li>S<sub>1</sub>(n+1) = h f(t<sub>n+1</sub>, Y<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> = (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> * </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>
* From this predicted vector, the corrected vector is computed as follows: * From this predicted vector, the corrected vector is computed as follows:
* <ul> * <ul>
* <li>y<sub>n+1</sub> = y<sub>n</sub> + S<sub>1</sub>(n+1) + [ -1 +1 -1 +1 ... &plusmn;1 ] r<sub>n+1</sub></li> * <li>y<sub>n+1</sub> = y<sub>n</sub> + S<sub>1</sub>(n+1) + [ -1 +1 -1 +1 ... &plusmn;1 ] r<sub>n+1</sub></li>
@ -153,32 +133,33 @@ import org.apache.commons.math.ode.nonstiff.AdamsMoultonIntegrator;
* @version $Revision$ $Date$ * @version $Revision$ $Date$
* @since 2.0 * @since 2.0
*/ */
public class NordsieckTransformer { public class AdamsNordsieckTransformer {
/** Cache for already computed coefficients. */ /** Cache for already computed coefficients. */
private static final Map<Integer, NordsieckTransformer> cache = private static final Map<Integer, AdamsNordsieckTransformer> cache =
new HashMap<Integer, NordsieckTransformer>(); new HashMap<Integer, AdamsNordsieckTransformer>();
/** Initialization matrix for the higher order derivatives wrt y'', y''' ... */ /** Initialization matrix for the higher order derivatives wrt y'', y''' ... */
private final RealMatrix initialization; private final Array2DRowRealMatrix initialization;
/** Update matrix for the higher order derivatives h<sup>2</sup>/2y'', h<sup>3</sup>/6 y''' ... */ /** Update matrix for the higher order derivatives h<sup>2</sup>/2y'', h<sup>3</sup>/6 y''' ... */
private final RealMatrix update; private final Array2DRowRealMatrix update;
/** Update coefficients of the higher order derivatives wrt y'. */ /** Update coefficients of the higher order derivatives wrt y'. */
private final double[] c1; private final double[] c1;
/** Simple constructor. /** Simple constructor.
* @param nSteps number of steps of the multistep method * @param nSteps number of steps of the multistep method
* (including the one being computed) * (excluding the one being computed)
*/ */
private NordsieckTransformer(final int nSteps) { private AdamsNordsieckTransformer(final int nSteps) {
// compute exact coefficients // compute exact coefficients
FieldMatrix<BigFraction> bigP = buildP(nSteps); FieldMatrix<BigFraction> bigP = buildP(nSteps);
FieldDecompositionSolver<BigFraction> pSolver = FieldDecompositionSolver<BigFraction> pSolver =
new FieldLUDecompositionImpl<BigFraction>(bigP).getSolver(); new FieldLUDecompositionImpl<BigFraction>(bigP).getSolver();
BigFraction[] u = new BigFraction[nSteps - 1];
BigFraction[] u = new BigFraction[nSteps];
Arrays.fill(u, BigFraction.ONE); Arrays.fill(u, BigFraction.ONE);
BigFraction[] bigC1 = pSolver.solve(u); BigFraction[] bigC1 = pSolver.solve(u);
@ -190,12 +171,12 @@ public class NordsieckTransformer {
// shift rows // shift rows
shiftedP[i] = shiftedP[i - 1]; shiftedP[i] = shiftedP[i - 1];
} }
shiftedP[0] = new BigFraction[nSteps - 1]; shiftedP[0] = new BigFraction[nSteps];
Arrays.fill(shiftedP[0], BigFraction.ZERO); Arrays.fill(shiftedP[0], BigFraction.ZERO);
FieldMatrix<BigFraction> bigMSupdate = FieldMatrix<BigFraction> bigMSupdate =
pSolver.solve(new Array2DRowFieldMatrix<BigFraction>(shiftedP, false)); pSolver.solve(new Array2DRowFieldMatrix<BigFraction>(shiftedP, false));
// initialization coefficients, computed from a Q matrix = abs(P) // initialization coefficients, computed from a R matrix = abs(P)
bigP.walkInOptimizedOrder(new DefaultFieldMatrixChangingVisitor<BigFraction>(BigFraction.ZERO) { bigP.walkInOptimizedOrder(new DefaultFieldMatrixChangingVisitor<BigFraction>(BigFraction.ZERO) {
/** {@inheritDoc} */ /** {@inheritDoc} */
@Override @Override
@ -203,14 +184,14 @@ public class NordsieckTransformer {
return ((column & 0x1) == 0x1) ? value : value.negate(); return ((column & 0x1) == 0x1) ? value : value.negate();
} }
}); });
FieldMatrix<BigFraction> bigQInverse = FieldMatrix<BigFraction> bigRInverse =
new FieldLUDecompositionImpl<BigFraction>(bigP).getSolver().getInverse(); new FieldLUDecompositionImpl<BigFraction>(bigP).getSolver().getInverse();
// convert coefficients to double // convert coefficients to double
initialization = MatrixUtils.bigFractionMatrixToRealMatrix(bigQInverse); initialization = MatrixUtils.bigFractionMatrixToRealMatrix(bigRInverse);
update = MatrixUtils.bigFractionMatrixToRealMatrix(bigMSupdate); update = MatrixUtils.bigFractionMatrixToRealMatrix(bigMSupdate);
c1 = new double[nSteps - 1]; c1 = new double[nSteps];
for (int i = 0; i < nSteps - 1; ++i) { for (int i = 0; i < nSteps; ++i) {
c1[i] = bigC1[i].doubleValue(); c1[i] = bigC1[i].doubleValue();
} }
@ -218,37 +199,45 @@ public class NordsieckTransformer {
/** Get the Nordsieck transformer for a given number of steps. /** Get the Nordsieck transformer for a given number of steps.
* @param nSteps number of steps of the multistep method * @param nSteps number of steps of the multistep method
* (including the one being computed) * (excluding the one being computed)
* @return Nordsieck transformer for the specified number of steps * @return Nordsieck transformer for the specified number of steps
*/ */
public static NordsieckTransformer getInstance(final int nSteps) { public static AdamsNordsieckTransformer getInstance(final int nSteps) {
synchronized(cache) { synchronized(cache) {
NordsieckTransformer t = cache.get(nSteps); AdamsNordsieckTransformer t = cache.get(nSteps);
if (t == null) { if (t == null) {
t = new NordsieckTransformer(nSteps); t = new AdamsNordsieckTransformer(nSteps);
cache.put(nSteps, t); cache.put(nSteps, t);
} }
return t; return t;
} }
} }
/** Build the P matrix transforming multistep to Nordsieck. /** Get the number of steps of the method
* <p> * (excluding the one being computed).
* Multistep representation uses y(k), s<sub>1</sub>(k), s<sub>1</sub>(k-1) ... s<sub>1</sub>(k-(n-1)). * @return number of steps of the method
* Nordsieck representation uses y(k), s<sub>1</sub>(k), s<sub>2</sub>(k) ... s<sub>n</sub>(k). * (excluding the one being computed)
* The two representations share their two first components y(k) and */
* s<sub>1</sub>(k). The P matrix is used to transform the remaining ones: public int getNSteps() {
return c1.length;
}
/** Build the P matrix.
* <p>The P matrix general terms are shifted j (-i)<sup>j-1</sup> terms:
* <pre> * <pre>
* [ s<sub>1</sub>(k-1) ... s<sub>1</sub>(k-(n-1)]<sup>T</sup> = s<sub>1</sub>(k) [1 ... 1]<sup>T</sup> + P [s<sub>2</sub>(k) ... s<sub>n</sub>(k)]<sup>T</sup> * [ -2 3 -4 5 ... ]
* </pre> * [ -4 12 -32 80 ... ]
* </p> * P = [ -6 27 -108 405 ... ]
* [ -8 48 -256 1280 ... ]
* [ ... ]
* </pre></p>
* @param nSteps number of steps of the multistep method * @param nSteps number of steps of the multistep method
* (including the one being computed) * (excluding the one being computed)
* @return P matrix * @return P matrix
*/ */
private FieldMatrix<BigFraction> buildP(final int nSteps) { private FieldMatrix<BigFraction> buildP(final int nSteps) {
final BigFraction[][] pData = new BigFraction[nSteps - 1][nSteps - 1]; final BigFraction[][] pData = new BigFraction[nSteps][nSteps];
for (int i = 0; i < pData.length; ++i) { for (int i = 0; i < pData.length; ++i) {
// build the P matrix elements from Taylor series formulas // build the P matrix elements from Taylor series formulas
@ -271,7 +260,7 @@ public class NordsieckTransformer {
* will be modified * will be modified
* @return high order derivatives at step start * @return high order derivatives at step start
*/ */
public RealMatrix initializeHighOrderDerivatives(final double[] first, public Array2DRowRealMatrix initializeHighOrderDerivatives(final double[] first,
final double[][] multistep) { final double[][] multistep) {
for (int i = 0; i < multistep.length; ++i) { for (int i = 0; i < multistep.length; ++i) {
final double[] msI = multistep[i]; final double[] msI = multistep[i];
@ -282,7 +271,7 @@ public class NordsieckTransformer {
return initialization.multiply(new Array2DRowRealMatrix(multistep, false)); return initialization.multiply(new Array2DRowRealMatrix(multistep, false));
} }
/** Update the high order scaled derivatives (phase 1). /** Update the high order scaled derivatives for Adams integrators (phase 1).
* <p>The complete update of high order derivatives has a form similar to: * <p>The complete update of high order derivatives has a form similar to:
* <pre> * <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> * 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>
@ -293,11 +282,11 @@ public class NordsieckTransformer {
* @return updated high order derivatives * @return updated high order derivatives
* @see #updateHighOrderDerivativesPhase2(double[], double[], RealMatrix) * @see #updateHighOrderDerivativesPhase2(double[], double[], RealMatrix)
*/ */
public RealMatrix updateHighOrderDerivativesPhase1(final RealMatrix highOrder) { public Array2DRowRealMatrix updateHighOrderDerivativesPhase1(final Array2DRowRealMatrix highOrder) {
return update.multiply(highOrder); return update.multiply(highOrder);
} }
/** Update the high order scaled derivatives (phase 2). /** Update the high order scaled derivatives Adams integrators (phase 2).
* <p>The complete update of high order derivatives has a form similar to: * <p>The complete update of high order derivatives has a form similar to:
* <pre> * <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> * 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>