001 /*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements. See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License. You may obtain a copy of the License at
008 *
009 * http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017
018 package org.apache.commons.math3.ode.nonstiff;
019
020 import java.util.Arrays;
021 import java.util.HashMap;
022 import java.util.Map;
023
024 import org.apache.commons.math3.fraction.BigFraction;
025 import org.apache.commons.math3.linear.Array2DRowFieldMatrix;
026 import org.apache.commons.math3.linear.Array2DRowRealMatrix;
027 import org.apache.commons.math3.linear.ArrayFieldVector;
028 import org.apache.commons.math3.linear.FieldDecompositionSolver;
029 import org.apache.commons.math3.linear.FieldLUDecomposition;
030 import org.apache.commons.math3.linear.FieldMatrix;
031 import org.apache.commons.math3.linear.MatrixUtils;
032 import org.apache.commons.math3.linear.QRDecomposition;
033 import org.apache.commons.math3.linear.RealMatrix;
034
035 /** Transformer to Nordsieck vectors for Adams integrators.
036 * <p>This class is used by {@link AdamsBashforthIntegrator Adams-Bashforth} and
037 * {@link AdamsMoultonIntegrator Adams-Moulton} integrators to convert between
038 * classical representation with several previous first derivatives and Nordsieck
039 * representation with higher order scaled derivatives.</p>
040 *
041 * <p>We define scaled derivatives s<sub>i</sub>(n) at step n as:
042 * <pre>
043 * s<sub>1</sub>(n) = h y'<sub>n</sub> for first derivative
044 * s<sub>2</sub>(n) = h<sup>2</sup>/2 y''<sub>n</sub> for second derivative
045 * s<sub>3</sub>(n) = h<sup>3</sup>/6 y'''<sub>n</sub> for third derivative
046 * ...
047 * s<sub>k</sub>(n) = h<sup>k</sup>/k! y<sup>(k)</sup><sub>n</sub> for k<sup>th</sup> derivative
048 * </pre></p>
049 *
050 * <p>With the previous definition, the classical representation of multistep methods
051 * uses first derivatives only, i.e. it handles y<sub>n</sub>, s<sub>1</sub>(n) and
052 * q<sub>n</sub> where q<sub>n</sub> is defined as:
053 * <pre>
054 * 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>
055 * </pre>
056 * (we omit the k index in the notation for clarity).</p>
057 *
058 * <p>Another possible representation uses the Nordsieck vector with
059 * higher degrees scaled derivatives all taken at the same step, i.e it handles y<sub>n</sub>,
060 * s<sub>1</sub>(n) and r<sub>n</sub>) where r<sub>n</sub> is defined as:
061 * <pre>
062 * r<sub>n</sub> = [ s<sub>2</sub>(n), s<sub>3</sub>(n) ... s<sub>k</sub>(n) ]<sup>T</sup>
063 * </pre>
064 * (here again we omit the k index in the notation for clarity)
065 * </p>
066 *
067 * <p>Taylor series formulas show that for any index offset i, s<sub>1</sub>(n-i) can be
068 * computed from s<sub>1</sub>(n), s<sub>2</sub>(n) ... s<sub>k</sub>(n), the formula being exact
069 * for degree k polynomials.
070 * <pre>
071 * s<sub>1</sub>(n-i) = s<sub>1</sub>(n) + ∑<sub>j>1</sub> j (-i)<sup>j-1</sup> s<sub>j</sub>(n)
072 * </pre>
073 * The previous formula can be used with several values for i to compute the transform between
074 * classical representation and Nordsieck vector at step end. The transform between r<sub>n</sub>
075 * and q<sub>n</sub> resulting from the Taylor series formulas above is:
076 * <pre>
077 * q<sub>n</sub> = s<sub>1</sub>(n) u + P r<sub>n</sub>
078 * </pre>
079 * where u is the [ 1 1 ... 1 ]<sup>T</sup> vector and P is the (k-1)×(k-1) matrix built
080 * with the j (-i)<sup>j-1</sup> terms:
081 * <pre>
082 * [ -2 3 -4 5 ... ]
083 * [ -4 12 -32 80 ... ]
084 * P = [ -6 27 -108 405 ... ]
085 * [ -8 48 -256 1280 ... ]
086 * [ ... ]
087 * </pre></p>
088 *
089 * <p>Changing -i into +i in the formula above can be used to compute a similar transform between
090 * classical representation and Nordsieck vector at step start. The resulting matrix is simply
091 * the absolute value of matrix P.</p>
092 *
093 * <p>For {@link AdamsBashforthIntegrator Adams-Bashforth} method, the Nordsieck vector
094 * at step n+1 is computed from the Nordsieck vector at step n as follows:
095 * <ul>
096 * <li>y<sub>n+1</sub> = y<sub>n</sub> + s<sub>1</sub>(n) + u<sup>T</sup> r<sub>n</sub></li>
097 * <li>s<sub>1</sub>(n+1) = h f(t<sub>n+1</sub>, y<sub>n+1</sub>)</li>
098 * <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>
099 * </ul>
100 * where A is a rows shifting matrix (the lower left part is an identity matrix):
101 * <pre>
102 * [ 0 0 ... 0 0 | 0 ]
103 * [ ---------------+---]
104 * [ 1 0 ... 0 0 | 0 ]
105 * A = [ 0 1 ... 0 0 | 0 ]
106 * [ ... | 0 ]
107 * [ 0 0 ... 1 0 | 0 ]
108 * [ 0 0 ... 0 1 | 0 ]
109 * </pre></p>
110 *
111 * <p>For {@link AdamsMoultonIntegrator Adams-Moulton} method, the predicted Nordsieck vector
112 * at step n+1 is computed from the Nordsieck vector at step n as follows:
113 * <ul>
114 * <li>Y<sub>n+1</sub> = y<sub>n</sub> + s<sub>1</sub>(n) + u<sup>T</sup> r<sub>n</sub></li>
115 * <li>S<sub>1</sub>(n+1) = h f(t<sub>n+1</sub>, Y<sub>n+1</sub>)</li>
116 * <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>
117 * </ul>
118 * From this predicted vector, the corrected vector is computed as follows:
119 * <ul>
120 * <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>
121 * <li>s<sub>1</sub>(n+1) = h f(t<sub>n+1</sub>, y<sub>n+1</sub>)</li>
122 * <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>
123 * </ul>
124 * where the upper case Y<sub>n+1</sub>, S<sub>1</sub>(n+1) and R<sub>n+1</sub> represent the
125 * predicted states whereas the lower case y<sub>n+1</sub>, s<sub>n+1</sub> and r<sub>n+1</sub>
126 * represent the corrected states.</p>
127 *
128 * <p>We observe that both methods use similar update formulas. In both cases a P<sup>-1</sup>u
129 * vector and a P<sup>-1</sup> A P matrix are used that do not depend on the state,
130 * they only depend on k. This class handles these transformations.</p>
131 *
132 * @version $Id: AdamsNordsieckTransformer.java 1416643 2012-12-03 19:37:14Z tn $
133 * @since 2.0
134 */
135 public class AdamsNordsieckTransformer {
136
137 /** Cache for already computed coefficients. */
138 private static final Map<Integer, AdamsNordsieckTransformer> CACHE =
139 new HashMap<Integer, AdamsNordsieckTransformer>();
140
141 /** Update matrix for the higher order derivatives h<sup>2</sup>/2y'', h<sup>3</sup>/6 y''' ... */
142 private final Array2DRowRealMatrix update;
143
144 /** Update coefficients of the higher order derivatives wrt y'. */
145 private final double[] c1;
146
147 /** Simple constructor.
148 * @param nSteps number of steps of the multistep method
149 * (excluding the one being computed)
150 */
151 private AdamsNordsieckTransformer(final int nSteps) {
152
153 // compute exact coefficients
154 FieldMatrix<BigFraction> bigP = buildP(nSteps);
155 FieldDecompositionSolver<BigFraction> pSolver =
156 new FieldLUDecomposition<BigFraction>(bigP).getSolver();
157
158 BigFraction[] u = new BigFraction[nSteps];
159 Arrays.fill(u, BigFraction.ONE);
160 BigFraction[] bigC1 = pSolver
161 .solve(new ArrayFieldVector<BigFraction>(u, false)).toArray();
162
163 // update coefficients are computed by combining transform from
164 // Nordsieck to multistep, then shifting rows to represent step advance
165 // then applying inverse transform
166 BigFraction[][] shiftedP = bigP.getData();
167 for (int i = shiftedP.length - 1; i > 0; --i) {
168 // shift rows
169 shiftedP[i] = shiftedP[i - 1];
170 }
171 shiftedP[0] = new BigFraction[nSteps];
172 Arrays.fill(shiftedP[0], BigFraction.ZERO);
173 FieldMatrix<BigFraction> bigMSupdate =
174 pSolver.solve(new Array2DRowFieldMatrix<BigFraction>(shiftedP, false));
175
176 // convert coefficients to double
177 update = MatrixUtils.bigFractionMatrixToRealMatrix(bigMSupdate);
178 c1 = new double[nSteps];
179 for (int i = 0; i < nSteps; ++i) {
180 c1[i] = bigC1[i].doubleValue();
181 }
182
183 }
184
185 /** Get the Nordsieck transformer for a given number of steps.
186 * @param nSteps number of steps of the multistep method
187 * (excluding the one being computed)
188 * @return Nordsieck transformer for the specified number of steps
189 */
190 public static AdamsNordsieckTransformer getInstance(final int nSteps) {
191 synchronized(CACHE) {
192 AdamsNordsieckTransformer t = CACHE.get(nSteps);
193 if (t == null) {
194 t = new AdamsNordsieckTransformer(nSteps);
195 CACHE.put(nSteps, t);
196 }
197 return t;
198 }
199 }
200
201 /** Get the number of steps of the method
202 * (excluding the one being computed).
203 * @return number of steps of the method
204 * (excluding the one being computed)
205 */
206 public int getNSteps() {
207 return c1.length;
208 }
209
210 /** Build the P matrix.
211 * <p>The P matrix general terms are shifted j (-i)<sup>j-1</sup> terms:
212 * <pre>
213 * [ -2 3 -4 5 ... ]
214 * [ -4 12 -32 80 ... ]
215 * P = [ -6 27 -108 405 ... ]
216 * [ -8 48 -256 1280 ... ]
217 * [ ... ]
218 * </pre></p>
219 * @param nSteps number of steps of the multistep method
220 * (excluding the one being computed)
221 * @return P matrix
222 */
223 private FieldMatrix<BigFraction> buildP(final int nSteps) {
224
225 final BigFraction[][] pData = new BigFraction[nSteps][nSteps];
226
227 for (int i = 0; i < pData.length; ++i) {
228 // build the P matrix elements from Taylor series formulas
229 final BigFraction[] pI = pData[i];
230 final int factor = -(i + 1);
231 int aj = factor;
232 for (int j = 0; j < pI.length; ++j) {
233 pI[j] = new BigFraction(aj * (j + 2));
234 aj *= factor;
235 }
236 }
237
238 return new Array2DRowFieldMatrix<BigFraction>(pData, false);
239
240 }
241
242 /** Initialize the high order scaled derivatives at step start.
243 * @param h step size to use for scaling
244 * @param t first steps times
245 * @param y first steps states
246 * @param yDot first steps derivatives
247 * @return Nordieck vector at first step (h<sup>2</sup>/2 y''<sub>n</sub>,
248 * h<sup>3</sup>/6 y'''<sub>n</sub> ... h<sup>k</sup>/k! y<sup>(k)</sup><sub>n</sub>)
249 */
250 public Array2DRowRealMatrix initializeHighOrderDerivatives(final double h, final double[] t,
251 final double[][] y,
252 final double[][] yDot) {
253
254 // using Taylor series with di = ti - t0, we get:
255 // y(ti) - y(t0) - di y'(t0) = di^2 / h^2 s2 + ... + di^k / h^k sk + O(h^(k+1))
256 // y'(ti) - y'(t0) = 2 di / h^2 s2 + ... + k di^(k-1) / h^k sk + O(h^k)
257 // we write these relations for i = 1 to i= n-1 as a set of 2(n-1) linear
258 // equations depending on the Nordsieck vector [s2 ... sk]
259 final double[][] a = new double[2 * (y.length - 1)][c1.length];
260 final double[][] b = new double[2 * (y.length - 1)][y[0].length];
261 final double[] y0 = y[0];
262 final double[] yDot0 = yDot[0];
263 for (int i = 1; i < y.length; ++i) {
264
265 final double di = t[i] - t[0];
266 final double ratio = di / h;
267 double dikM1Ohk = 1 / h;
268
269 // linear coefficients of equations
270 // y(ti) - y(t0) - di y'(t0) and y'(ti) - y'(t0)
271 final double[] aI = a[2 * i - 2];
272 final double[] aDotI = a[2 * i - 1];
273 for (int j = 0; j < aI.length; ++j) {
274 dikM1Ohk *= ratio;
275 aI[j] = di * dikM1Ohk;
276 aDotI[j] = (j + 2) * dikM1Ohk;
277 }
278
279 // expected value of the previous equations
280 final double[] yI = y[i];
281 final double[] yDotI = yDot[i];
282 final double[] bI = b[2 * i - 2];
283 final double[] bDotI = b[2 * i - 1];
284 for (int j = 0; j < yI.length; ++j) {
285 bI[j] = yI[j] - y0[j] - di * yDot0[j];
286 bDotI[j] = yDotI[j] - yDot0[j];
287 }
288
289 }
290
291 // solve the rectangular system in the least square sense
292 // to get the best estimate of the Nordsieck vector [s2 ... sk]
293 QRDecomposition decomposition;
294 decomposition = new QRDecomposition(new Array2DRowRealMatrix(a, false));
295 RealMatrix x = decomposition.getSolver().solve(new Array2DRowRealMatrix(b, false));
296 return new Array2DRowRealMatrix(x.getData(), false);
297 }
298
299 /** Update the high order scaled derivatives for Adams integrators (phase 1).
300 * <p>The complete update of high order derivatives has a form similar to:
301 * <pre>
302 * 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>
303 * </pre>
304 * this method computes the P<sup>-1</sup> A P r<sub>n</sub> part.</p>
305 * @param highOrder high order scaled derivatives
306 * (h<sup>2</sup>/2 y'', ... h<sup>k</sup>/k! y(k))
307 * @return updated high order derivatives
308 * @see #updateHighOrderDerivativesPhase2(double[], double[], Array2DRowRealMatrix)
309 */
310 public Array2DRowRealMatrix updateHighOrderDerivativesPhase1(final Array2DRowRealMatrix highOrder) {
311 return update.multiply(highOrder);
312 }
313
314 /** Update the high order scaled derivatives Adams integrators (phase 2).
315 * <p>The complete update of high order derivatives has a form similar to:
316 * <pre>
317 * 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>
318 * </pre>
319 * this method computes the (s<sub>1</sub>(n) - s<sub>1</sub>(n+1)) P<sup>-1</sup> u part.</p>
320 * <p>Phase 1 of the update must already have been performed.</p>
321 * @param start first order scaled derivatives at step start
322 * @param end first order scaled derivatives at step end
323 * @param highOrder high order scaled derivatives, will be modified
324 * (h<sup>2</sup>/2 y'', ... h<sup>k</sup>/k! y(k))
325 * @see #updateHighOrderDerivativesPhase1(Array2DRowRealMatrix)
326 */
327 public void updateHighOrderDerivativesPhase2(final double[] start,
328 final double[] end,
329 final Array2DRowRealMatrix highOrder) {
330 final double[][] data = highOrder.getDataRef();
331 for (int i = 0; i < data.length; ++i) {
332 final double[] dataI = data[i];
333 final double c1I = c1[i];
334 for (int j = 0; j < dataI.length; ++j) {
335 dataI[j] += c1I * (start[j] - end[j]);
336 }
337 }
338 }
339
340 }