Commons Math example source code file (AdamsNordsieckTransformer.java)
The Commons Math AdamsNordsieckTransformer.java source code/*
* 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 java.util.Arrays;
import java.util.HashMap;
import java.util.Map;
import org.apache.commons.math.fraction.BigFraction;
import org.apache.commons.math.linear.Array2DRowFieldMatrix;
import org.apache.commons.math.linear.Array2DRowRealMatrix;
import org.apache.commons.math.linear.DefaultFieldMatrixChangingVisitor;
import org.apache.commons.math.linear.FieldDecompositionSolver;
import org.apache.commons.math.linear.FieldLUDecompositionImpl;
import org.apache.commons.math.linear.FieldMatrix;
import org.apache.commons.math.linear.MatrixUtils;
/** Transformer to Nordsieck vectors for Adams integrators.
* <p>This class i 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 si(n) at step n as:
* <pre>
* s<sub>1(n) = h y'n for first derivative
* s<sub>2(n) = h2/2 y''n for second derivative
* s<sub>3(n) = h3/6 y'''n for third derivative
* ...
* s<sub>k(n) = hk/k! y(k)n for kth derivative
* </pre>
*
* <p>With the previous definition, the classical representation of multistep methods
* uses first derivatives only, i.e. it handles y<sub>n, s1(n) and
* q<sub>n where qn is defined as:
* <pre>
* q<sub>n = [ s1(n-1) s1(n-2) ... s1(n-(k-1)) ]T
* </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,
* s<sub>1(n) and rn) where rn is defined as:
* <pre>
* r<sub>n = [ s2(n), s3(n) ... sk(n) ]T
* </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, s1(n-i) can be
* computed from s<sub>1(n), s2(n) ... sk(n), the formula being exact
* for degree k polynomials.
* <pre>
* s<sub>1(n-i) = s1(n) + ?j j (-i)j-1 sj(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
* and q<sub>n resulting from the Taylor series formulas above is:
* <pre>
* q<sub>n = s1(n) u + P rn
* </pre>
* where u is the [ 1 1 ... 1 ]<sup>T vector and P is the (k-1)×(k-1) matrix built
* with the j (-i)<sup>j-1 terms:
* <pre>
* [ -2 3 -4 5 ... ]
* [ -4 12 -32 80 ... ]
* P = [ -6 27 -108 405 ... ]
* [ -8 48 -256 1280 ... ]
* [ ... ]
* </pre>
*
* <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>yn+1 = yn + s1(n) + uT rn
* <li>s1(n+1) = h f(tn+1, yn+1)
* <li>rn+1 = (s1(n) - s1(n+1)) P-1 u + P-1 A P rn
* </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>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>Yn+1 = yn + s1(n) + uT rn
* <li>S1(n+1) = h f(tn+1, Yn+1)
* <li>Rn+1 = (s1(n) - s1(n+1)) P-1 u + P-1 A P rn
* </ul>
* From this predicted vector, the corrected vector is computed as follows:
* <ul>
* <li>yn+1 = yn + S1(n+1) + [ -1 +1 -1 +1 ... ±1 ] rn+1
* <li>s1(n+1) = h f(tn+1, yn+1)
* <li>rn+1 = Rn+1 + (s1(n+1) - S1(n+1)) P-1 u
* </ul>
* where the upper case Y<sub>n+1, S1(n+1) and Rn+1 represent the
* predicted states whereas the lower case y<sub>n+1, sn+1 and rn+1
* represent the corrected states.</p>
*
* <p>We observe that both methods use similar update formulas. In both cases a P-1u
* vector and a P<sup>-1 A P matrix are used that do not depend on the state,
* they only depend on k. This class handles these transformations.</p>
*
* @version $Revision: 810196 $ $Date: 2009-09-01 15:47:46 -0400 (Tue, 01 Sep 2009) $
* @since 2.0
*/
public class AdamsNordsieckTransformer {
/** Cache for already computed coefficients. */
private static final Map<Integer, AdamsNordsieckTransformer> CACHE =
new HashMap<Integer, AdamsNordsieckTransformer>();
/** Initialization matrix for the higher order derivatives wrt y'', y''' ... */
private final Array2DRowRealMatrix initialization;
/** Update matrix for the higher order derivatives h<sup>2/2y'', h3/6 y''' ... */
private final Array2DRowRealMatrix update;
/** Update coefficients of the higher order derivatives wrt y'. */
private final double[] c1;
/** Simple constructor.
* @param nSteps number of steps of the multistep method
* (excluding the one being computed)
*/
private AdamsNordsieckTransformer(final int nSteps) {
// compute exact coefficients
FieldMatrix<BigFraction> bigP = buildP(nSteps);
FieldDecompositionSolver<BigFraction> pSolver =
new FieldLUDecompositionImpl<BigFraction>(bigP).getSolver();
BigFraction[] u = new BigFraction[nSteps];
Arrays.fill(u, BigFraction.ONE);
BigFraction[] bigC1 = pSolver.solve(u);
// update coefficients are computed by combining transform from
// Nordsieck to multistep, then shifting rows to represent step advance
// then applying inverse transform
BigFraction[][] shiftedP = bigP.getData();
for (int i = shiftedP.length - 1; i > 0; --i) {
// shift rows
shiftedP[i] = shiftedP[i - 1];
}
shiftedP[0] = new BigFraction[nSteps];
Arrays.fill(shiftedP[0], BigFraction.ZERO);
FieldMatrix<BigFraction> bigMSupdate =
pSolver.solve(new Array2DRowFieldMatrix<BigFraction>(shiftedP, false));
// initialization coefficients, computed from a R matrix = abs(P)
bigP.walkInOptimizedOrder(new DefaultFieldMatrixChangingVisitor<BigFraction>(BigFraction.ZERO) {
/** {@inheritDoc} */
@Override
public BigFraction visit(int row, int column, BigFraction value) {
return ((column & 0x1) == 0x1) ? value : value.negate();
}
});
FieldMatrix<BigFraction> bigRInverse =
new FieldLUDecompositionImpl<BigFraction>(bigP).getSolver().getInverse();
// convert coefficients to double
initialization = MatrixUtils.bigFractionMatrixToRealMatrix(bigRInverse);
update = MatrixUtils.bigFractionMatrixToRealMatrix(bigMSupdate);
c1 = new double[nSteps];
for (int i = 0; i < nSteps; ++i) {
c1[i] = bigC1[i].doubleValue();
}
}
/** Get the Nordsieck transformer for a given number of steps.
* @param nSteps number of steps of the multistep method
* (excluding the one being computed)
* @return Nordsieck transformer for the specified number of steps
*/
public static AdamsNordsieckTransformer getInstance(final int nSteps) {
synchronized(CACHE) {
AdamsNordsieckTransformer t = CACHE.get(nSteps);
if (t == null) {
t = new AdamsNordsieckTransformer(nSteps);
CACHE.put(nSteps, t);
}
return t;
}
}
/** Get the number of steps of the method
* (excluding the one being computed).
* @return number of steps of the method
* (excluding the one being computed)
*/
public int getNSteps() {
return c1.length;
}
/** Build the P matrix.
* <p>The P matrix general terms are shifted j (-i)j-1 terms:
* <pre>
* [ -2 3 -4 5 ... ]
* [ -4 12 -32 80 ... ]
* P = [ -6 27 -108 405 ... ]
* [ -8 48 -256 1280 ... ]
* [ ... ]
* </pre>
* @param nSteps number of steps of the multistep method
* (excluding the one being computed)
* @return P matrix
*/
private FieldMatrix<BigFraction> buildP(final int nSteps) {
final BigFraction[][] pData = new BigFraction[nSteps][nSteps];
for (int i = 0; i < pData.length; ++i) {
// build the P matrix elements from Taylor series formulas
final BigFraction[] pI = pData[i];
final int factor = -(i + 1);
int aj = factor;
for (int j = 0; j < pI.length; ++j) {
pI[j] = new BigFraction(aj * (j + 2));
aj *= factor;
}
}
return new Array2DRowFieldMatrix<BigFraction>(pData, false);
}
/** 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
*/
public Array2DRowRealMatrix initializeHighOrderDerivatives(final double[] first,
final double[][] multistep) {
for (int i = 0; i < multistep.length; ++i) {
final double[] msI = multistep[i];
for (int j = 0; j < first.length; ++j) {
msI[j] -= first[j];
}
}
return initialization.multiply(new Array2DRowRealMatrix(multistep, false));
}
/** 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 = (s1(n) - s1(n+1)) P-1 u + P-1 A P rn
* </pre>
* this method computes the P<sup>-1 A P rn part.
* @param highOrder high order scaled derivatives
* (h<sup>2/2 y'', ... hk/k! y(k))
* @return updated high order derivatives
* @see #updateHighOrderDerivativesPhase2(double[], double[], Array2DRowRealMatrix)
*/
public Array2DRowRealMatrix updateHighOrderDerivativesPhase1(final Array2DRowRealMatrix 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 = (s1(n) - s1(n+1)) P-1 u + P-1 A P rn
* </pre>
* this method computes the (s<sub>1(n) - s1(n+1)) P-1 u part.
* <p>Phase 1 of the update must already have been performed.
* @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/2 y'', ... hk/k! y(k))
* @see #updateHighOrderDerivativesPhase1(Array2DRowRealMatrix)
*/
public void updateHighOrderDerivativesPhase2(final double[] start,
final double[] end,
final Array2DRowRealMatrix highOrder) {
final double[][] data = highOrder.getDataRef();
for (int i = 0; i < data.length; ++i) {
final double[] dataI = data[i];
final double c1I = c1[i];
for (int j = 0; j < dataI.length; ++j) {
dataI[j] += c1I * (start[j] - end[j]);
}
}
}
}
Other Commons Math examples (source code examples)Here is a short list of links related to this Commons Math AdamsNordsieckTransformer.java source code file: |
The search page
Other Commons Math source code examples at this package level
Click here to learn more about this project
