org.apache.commons.math3.ode.nonstiff.AdamsNordsieckFieldTransformer Maven / Gradle / Ivy
Go to download
Show more of this group Show more artifacts with this name
Show all versions of virtdata-lib-realer Show documentation
Show all versions of virtdata-lib-realer Show documentation
With inspiration from other libraries
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math3.ode.nonstiff;
import java.util.Arrays;
import java.util.HashMap;
import java.util.Map;
import org.apache.commons.math3.Field;
import org.apache.commons.math3.RealFieldElement;
import org.apache.commons.math3.linear.Array2DRowFieldMatrix;
import org.apache.commons.math3.linear.ArrayFieldVector;
import org.apache.commons.math3.linear.FieldDecompositionSolver;
import org.apache.commons.math3.linear.FieldLUDecomposition;
import org.apache.commons.math3.linear.FieldMatrix;
import org.apache.commons.math3.util.MathArrays;
/** Transformer to Nordsieck vectors for Adams integrators.
* 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.
*
* We define scaled derivatives si(n) at step n as:
*
* s1(n) = h y'n for first derivative
* s2(n) = h2/2 y''n for second derivative
* s3(n) = h3/6 y'''n for third derivative
* ...
* sk(n) = hk/k! y(k)n for kth derivative
*
*
* With the previous definition, the classical representation of multistep methods
* uses first derivatives only, i.e. it handles yn, s1(n) and
* qn where qn is defined as:
*
* qn = [ s1(n-1) s1(n-2) ... s1(n-(k-1)) ]T
*
* (we omit the k index in the notation for clarity).
*
* Another possible representation uses the Nordsieck vector with
* higher degrees scaled derivatives all taken at the same step, i.e it handles yn,
* s1(n) and rn) where rn is defined as:
*
* rn = [ s2(n), s3(n) ... sk(n) ]T
*
* (here again we omit the k index in the notation for clarity)
*
*
* Taylor series formulas show that for any index offset i, s1(n-i) can be
* computed from s1(n), s2(n) ... sk(n), the formula being exact
* for degree k polynomials.
*
* s1(n-i) = s1(n) + ∑j>0 (j+1) (-i)j sj+1(n)
*
* 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 rn
* and qn resulting from the Taylor series formulas above is:
*
* qn = s1(n) u + P rn
*
* where u is the [ 1 1 ... 1 ]T vector and P is the (k-1)×(k-1) matrix built
* with the (j+1) (-i)j terms with i being the row number starting from 1 and j being
* the column number starting from 1:
*
* [ -2 3 -4 5 ... ]
* [ -4 12 -32 80 ... ]
* P = [ -6 27 -108 405 ... ]
* [ -8 48 -256 1280 ... ]
* [ ... ]
*
*
* 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.
*
* For {@link AdamsBashforthIntegrator Adams-Bashforth} method, the Nordsieck vector
* at step n+1 is computed from the Nordsieck vector at step n as follows:
*
* - yn+1 = yn + s1(n) + uT rn
* - s1(n+1) = h f(tn+1, yn+1)
* - rn+1 = (s1(n) - s1(n+1)) P-1 u + P-1 A P rn
*
* where A is a rows shifting matrix (the lower left part is an identity matrix):
*
* [ 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 ]
*
*
* 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:
*
* - Yn+1 = yn + s1(n) + uT rn
* - S1(n+1) = h f(tn+1, Yn+1)
* - Rn+1 = (s1(n) - s1(n+1)) P-1 u + P-1 A P rn
*
* From this predicted vector, the corrected vector is computed as follows:
*
* - yn+1 = yn + S1(n+1) + [ -1 +1 -1 +1 ... ±1 ] rn+1
* - s1(n+1) = h f(tn+1, yn+1)
* - rn+1 = Rn+1 + (s1(n+1) - S1(n+1)) P-1 u
*
* where the upper case Yn+1, S1(n+1) and Rn+1 represent the
* predicted states whereas the lower case yn+1, sn+1 and rn+1
* represent the corrected states.
*
* We observe that both methods use similar update formulas. In both cases a P-1u
* vector and a P-1 A P matrix are used that do not depend on the state,
* they only depend on k. This class handles these transformations.
*
* @param the type of the field elements
* @since 3.6
*/
public class AdamsNordsieckFieldTransformer> {
/** Cache for already computed coefficients. */
private static final Map>,
AdamsNordsieckFieldTransformer>>> CACHE =
new HashMap>,
AdamsNordsieckFieldTransformer>>>();
/** Field to which the time and state vector elements belong. */
private final Field field;
/** Update matrix for the higher order derivatives h2/2 y'', h3/6 y''' ... */
private final Array2DRowFieldMatrix 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 field, final int n) {
this.field = field;
final int rows = n - 1;
// compute coefficients
FieldMatrix bigP = buildP(rows);
FieldDecompositionSolver pSolver =
new FieldLUDecomposition(bigP).getSolver();
T[] u = MathArrays.buildArray(field, rows);
Arrays.fill(u, field.getOne());
c1 = pSolver.solve(new ArrayFieldVector(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(pSolver.solve(new Array2DRowFieldMatrix(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 the type of the field elements
*/
@SuppressWarnings("unchecked")
public static > AdamsNordsieckFieldTransformer
getInstance(final Field field, final int nSteps) {
synchronized(CACHE) {
Map>,
AdamsNordsieckFieldTransformer>> map = CACHE.get(nSteps);
if (map == null) {
map = new HashMap>,
AdamsNordsieckFieldTransformer>>();
CACHE.put(nSteps, map);
}
@SuppressWarnings("rawtypes") // use rawtype to avoid compilation problems with java 1.5
AdamsNordsieckFieldTransformer t = map.get(field);
if (t == null) {
t = new AdamsNordsieckFieldTransformer(field, nSteps);
map.put(field, (AdamsNordsieckFieldTransformer) t);
}
return (AdamsNordsieckFieldTransformer) t;
}
}
/** Build the P matrix.
* The P matrix general terms are shifted (j+1) (-i)j terms
* with i being the row number starting from 1 and j being the column
* number starting from 1:
*
* [ -2 3 -4 5 ... ]
* [ -4 12 -32 80 ... ]
* P = [ -6 27 -108 405 ... ]
* [ -8 48 -256 1280 ... ]
* [ ... ]
*
* @param rows number of rows of the matrix
* @return P matrix
*/
private FieldMatrix 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(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 (h2/2 y''n,
* h3/6 y'''n ... hk/k! y(k)n)
*/
public Array2DRowFieldMatrix 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 decomposition = new FieldLUDecomposition(new Array2DRowFieldMatrix(a, false));
final FieldMatrix x = decomposition.getSolver().solve(new Array2DRowFieldMatrix(b, false));
// extract just the Nordsieck vector [s2 ... sk]
final Array2DRowFieldMatrix truncatedX =
new Array2DRowFieldMatrix(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).
* The complete update of high order derivatives has a form similar to:
*
* rn+1 = (s1(n) - s1(n+1)) P-1 u + P-1 A P rn
*
* this method computes the P-1 A P rn part.
* @param highOrder high order scaled derivatives
* (h2/2 y'', ... hk/k! y(k))
* @return updated high order derivatives
* @see #updateHighOrderDerivativesPhase2(RealFieldElement[], RealFieldElement[], Array2DRowFieldMatrix)
*/
public Array2DRowFieldMatrix updateHighOrderDerivativesPhase1(final Array2DRowFieldMatrix highOrder) {
return update.multiply(highOrder);
}
/** Update the high order scaled derivatives Adams integrators (phase 2).
* The complete update of high order derivatives has a form similar to:
*
* rn+1 = (s1(n) - s1(n+1)) P-1 u + P-1 A P rn
*
* this method computes the (s1(n) - s1(n+1)) P-1 u part.
* 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
* (h2/2 y'', ... hk/k! y(k))
* @see #updateHighOrderDerivativesPhase1(Array2DRowFieldMatrix)
*/
public void updateHighOrderDerivativesPhase2(final T[] start,
final T[] end,
final Array2DRowFieldMatrix 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])));
}
}
}
}