| /* |
| * 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.legacy.ode.nonstiff; |
| |
| import java.util.Arrays; |
| import java.util.HashMap; |
| import java.util.Map; |
| |
| import org.apache.commons.math4.legacy.core.Field; |
| import org.apache.commons.math4.legacy.core.RealFieldElement; |
| import org.apache.commons.math4.legacy.linear.Array2DRowFieldMatrix; |
| import org.apache.commons.math4.legacy.linear.ArrayFieldVector; |
| import org.apache.commons.math4.legacy.linear.FieldDecompositionSolver; |
| import org.apache.commons.math4.legacy.linear.FieldLUDecomposition; |
| import org.apache.commons.math4.legacy.linear.FieldMatrix; |
| import org.apache.commons.math4.legacy.core.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: |
| * <div style="white-space: pre"><code> |
| * 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 |
| * </code></div> |
| * |
| * <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: |
| * <div style="white-space: pre"><code> |
| * 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> |
| * </code></div> |
| * (we omit the k index in the notation for clarity). |
| * |
| * <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: |
| * <div style="white-space: pre"><code> |
| * r<sub>n</sub> = [ s<sub>2</sub>(n), s<sub>3</sub>(n) ... s<sub>k</sub>(n) ]<sup>T</sup> |
| * </code></div> |
| * (here again we omit the k index in the notation for clarity) |
| * |
| * <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. |
| * <div style="white-space: pre"><code> |
| * 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) |
| * </code></div> |
| * 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: |
| * <div style="white-space: pre"><code> |
| * q<sub>n</sub> = s<sub>1</sub>(n) u + P r<sub>n</sub> |
| * </code></div> |
| * 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>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>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>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 final 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<>(); |
| |
| /** 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<>(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 <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<>(); |
| CACHE.put(nSteps, map); |
| } |
| @SuppressWarnings("unchecked") |
| AdamsNordsieckFieldTransformer<T> t = (AdamsNordsieckFieldTransformer<T>) map.get(field); |
| if (t == null) { |
| t = new AdamsNordsieckFieldTransformer<>(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> |
| * @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<>(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<>(new Array2DRowFieldMatrix<>(a, false)); |
| final FieldMatrix<T> x = decomposition.getSolver().solve(new Array2DRowFieldMatrix<>(b, false)); |
| |
| // extract just the Nordsieck vector [s2 ... sk] |
| final Array2DRowFieldMatrix<T> 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). |
| * <p>The complete update of high order derivatives has a form similar to: |
| * <div style="white-space: pre"><code> |
| * 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> |
| * </code></div> |
| * this method computes the P<sup>-1</sup> A P r<sub>n</sub> part. |
| * @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(RealFieldElement[], RealFieldElement[], 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: |
| * <div style="white-space: pre"><code> |
| * 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> |
| * </code></div> |
| * this method computes the (s<sub>1</sub>(n) - s<sub>1</sub>(n+1)) P<sup>-1</sup> u part. |
| * <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]))); |
| } |
| } |
| } |
| |
| } |
