/* * 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.fraction.BigFraction; import org.apache.commons.math3.linear.Array2DRowFieldMatrix; import org.apache.commons.math3.linear.Array2DRowRealMatrix; 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.linear.MatrixUtils; import org.apache.commons.math3.linear.QRDecomposition; import org.apache.commons.math3.linear.RealMatrix; /** 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: *
* [ 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: *
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.
* * @since 2.0 */ public class AdamsNordsieckTransformer { /** Cache for already computed coefficients. */ private static final MapThe 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
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(double[], double[], Array2DRowRealMatrix) */ public Array2DRowRealMatrix updateHighOrderDerivativesPhase1(final Array2DRowRealMatrix 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(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]); } } } }