1 /*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17
18 package org.apache.commons.math4.legacy.ode.nonstiff;
19
20 import java.util.HashMap;
21 import java.util.Map;
22
23 import org.apache.commons.numbers.fraction.BigFraction;
24 import org.apache.commons.numbers.field.BigFractionField;
25 import org.apache.commons.math4.legacy.linear.Array2DRowRealMatrix;
26 import org.apache.commons.math4.legacy.linear.QRDecomposition;
27 import org.apache.commons.math4.legacy.linear.RealMatrix;
28 import org.apache.commons.math4.legacy.field.linalg.FieldDenseMatrix;
29 import org.apache.commons.math4.legacy.field.linalg.FieldDecompositionSolver;
30 import org.apache.commons.math4.legacy.field.linalg.FieldLUDecomposition;
31
32 /** Transformer to Nordsieck vectors for Adams integrators.
33 * <p>This class is used by {@link AdamsBashforthIntegrator Adams-Bashforth} and
34 * {@link AdamsMoultonIntegrator Adams-Moulton} integrators to convert between
35 * classical representation with several previous first derivatives and Nordsieck
36 * representation with higher order scaled derivatives.</p>
37 *
38 * <p>We define scaled derivatives s<sub>i</sub>(n) at step n as:
39 * <div style="white-space: pre"><code>
40 * s<sub>1</sub>(n) = h y'<sub>n</sub> for first derivative
41 * s<sub>2</sub>(n) = h<sup>2</sup>/2 y''<sub>n</sub> for second derivative
42 * s<sub>3</sub>(n) = h<sup>3</sup>/6 y'''<sub>n</sub> for third derivative
43 * ...
44 * s<sub>k</sub>(n) = h<sup>k</sup>/k! y<sup>(k)</sup><sub>n</sub> for k<sup>th</sup> derivative
45 * </code></div>
46 *
47 * <p>With the previous definition, the classical representation of multistep methods
48 * uses first derivatives only, i.e. it handles y<sub>n</sub>, s<sub>1</sub>(n) and
49 * q<sub>n</sub> where q<sub>n</sub> is defined as:
50 * <div style="white-space: pre"><code>
51 * 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>
52 * </code></div>
53 * (we omit the k index in the notation for clarity).
54 *
55 * <p>Another possible representation uses the Nordsieck vector with
56 * higher degrees scaled derivatives all taken at the same step, i.e it handles y<sub>n</sub>,
57 * s<sub>1</sub>(n) and r<sub>n</sub>) where r<sub>n</sub> is defined as:
58 * <div style="white-space: pre"><code>
59 * r<sub>n</sub> = [ s<sub>2</sub>(n), s<sub>3</sub>(n) ... s<sub>k</sub>(n) ]<sup>T</sup>
60 * </code></div>
61 * (here again we omit the k index in the notation for clarity)
62 *
63 * <p>Taylor series formulas show that for any index offset i, s<sub>1</sub>(n-i) can be
64 * computed from s<sub>1</sub>(n), s<sub>2</sub>(n) ... s<sub>k</sub>(n), the formula being exact
65 * for degree k polynomials.
66 * <div style="white-space: pre"><code>
67 * 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)
68 * </code></div>
69 * The previous formula can be used with several values for i to compute the transform between
70 * classical representation and Nordsieck vector at step end. The transform between r<sub>n</sub>
71 * and q<sub>n</sub> resulting from the Taylor series formulas above is:
72 * <div style="white-space: pre"><code>
73 * q<sub>n</sub> = s<sub>1</sub>(n) u + P r<sub>n</sub>
74 * </code></div>
75 * where u is the [ 1 1 ... 1 ]<sup>T</sup> vector and P is the (k-1)×(k-1) matrix built
76 * with the (j+1) (-i)<sup>j</sup> terms with i being the row number starting from 1 and j being
77 * the column number starting from 1:
78 * <pre>
79 * [ -2 3 -4 5 ... ]
80 * [ -4 12 -32 80 ... ]
81 * P = [ -6 27 -108 405 ... ]
82 * [ -8 48 -256 1280 ... ]
83 * [ ... ]
84 * </pre>
85 *
86 * <p>Changing -i into +i in the formula above can be used to compute a similar transform between
87 * classical representation and Nordsieck vector at step start. The resulting matrix is simply
88 * the absolute value of matrix P.</p>
89 *
90 * <p>For {@link AdamsBashforthIntegrator Adams-Bashforth} method, the Nordsieck vector
91 * at step n+1 is computed from the Nordsieck vector at step n as follows:
92 * <ul>
93 * <li>y<sub>n+1</sub> = y<sub>n</sub> + s<sub>1</sub>(n) + u<sup>T</sup> r<sub>n</sub></li>
94 * <li>s<sub>1</sub>(n+1) = h f(t<sub>n+1</sub>, y<sub>n+1</sub>)</li>
95 * <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>
96 * </ul>
97 * where A is a rows shifting matrix (the lower left part is an identity matrix):
98 * <pre>
99 * [ 0 0 ... 0 0 | 0 ]
100 * [ ---------------+---]
101 * [ 1 0 ... 0 0 | 0 ]
102 * A = [ 0 1 ... 0 0 | 0 ]
103 * [ ... | 0 ]
104 * [ 0 0 ... 1 0 | 0 ]
105 * [ 0 0 ... 0 1 | 0 ]
106 * </pre>
107 *
108 * <p>For {@link AdamsMoultonIntegrator Adams-Moulton} method, the predicted Nordsieck vector
109 * at step n+1 is computed from the Nordsieck vector at step n as follows:
110 * <ul>
111 * <li>Y<sub>n+1</sub> = y<sub>n</sub> + s<sub>1</sub>(n) + u<sup>T</sup> r<sub>n</sub></li>
112 * <li>S<sub>1</sub>(n+1) = h f(t<sub>n+1</sub>, Y<sub>n+1</sub>)</li>
113 * <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>
114 * </ul>
115 * From this predicted vector, the corrected vector is computed as follows:
116 * <ul>
117 * <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>
118 * <li>s<sub>1</sub>(n+1) = h f(t<sub>n+1</sub>, y<sub>n+1</sub>)</li>
119 * <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>
120 * </ul>
121 * where the upper case Y<sub>n+1</sub>, S<sub>1</sub>(n+1) and R<sub>n+1</sub> represent the
122 * predicted states whereas the lower case y<sub>n+1</sub>, s<sub>n+1</sub> and r<sub>n+1</sub>
123 * represent the corrected states.
124 *
125 * <p>We observe that both methods use similar update formulas. In both cases a P<sup>-1</sup>u
126 * vector and a P<sup>-1</sup> A P matrix are used that do not depend on the state,
127 * they only depend on k. This class handles these transformations.</p>
128 *
129 * @since 2.0
130 */
131 public final class AdamsNordsieckTransformer {
132
133 /** Cache for already computed coefficients. */
134 private static final Map<Integer, AdamsNordsieckTransformer> CACHE =
135 new HashMap<>();
136
137 /** Update matrix for the higher order derivatives h<sup>2</sup>/2 y'', h<sup>3</sup>/6 y''' ... */
138 private final Array2DRowRealMatrix update;
139
140 /** Update coefficients of the higher order derivatives wrt y'. */
141 private final double[] c1;
142
143 /** Simple constructor.
144 * @param n number of steps of the multistep method
145 * (excluding the one being computed)
146 */
147 private AdamsNordsieckTransformer(final int n) {
148 final int dim = n - 1;
149
150 // compute exact coefficients
151 final FieldDenseMatrix<BigFraction> bigP = buildP(dim);
152 final FieldDecompositionSolver<BigFraction> pSolver = FieldLUDecomposition.of(bigP).getSolver();
153
154 final FieldDenseMatrix<BigFraction> u = FieldDenseMatrix.create(BigFractionField.get(), dim, 1)
155 .fill(BigFraction.ONE);
156 final FieldDenseMatrix<BigFraction> bigC1 = pSolver.solve(u);
157
158 // update coefficients are computed by combining transform from
159 // Nordsieck to multistep, then shifting rows to represent step advance
160 // then applying inverse transform
161 final FieldDenseMatrix<BigFraction> shiftedP = bigP.copy();
162 for (int i = dim - 1; i > 0; --i) {
163 // shift rows
164 for (int j = 0; j < dim; j++) {
165 shiftedP.set(i, j, shiftedP.get(i - 1, j));
166 }
167 }
168 for (int j = 0; j < dim; j++) {
169 shiftedP.set(0, j, BigFraction.ZERO);
170 }
171
172 final FieldDenseMatrix<BigFraction> bigMSupdate = pSolver.solve(shiftedP);
173
174 // convert coefficients to double
175 final double[][] updateData = new double[dim][dim];
176 for (int i = 0; i < dim; i++) {
177 for (int j = 0; j < dim; j++) {
178 updateData[i][j] = bigMSupdate.get(i, j).doubleValue();
179 }
180 }
181
182 update = new Array2DRowRealMatrix(updateData, false);
183 c1 = new double[dim];
184 for (int i = 0; i < dim; ++i) {
185 c1[i] = bigC1.get(i, 0).doubleValue();
186 }
187 }
188
189 /** Get the Nordsieck transformer for a given number of steps.
190 * @param nSteps number of steps of the multistep method
191 * (excluding the one being computed)
192 * @return Nordsieck transformer for the specified number of steps
193 */
194 public static AdamsNordsieckTransformer getInstance(final int nSteps) {
195 synchronized(CACHE) {
196 AdamsNordsieckTransformer t = CACHE.get(nSteps);
197 if (t == null) {
198 t = new AdamsNordsieckTransformer(nSteps);
199 CACHE.put(nSteps, t);
200 }
201 return t;
202 }
203 }
204
205 /** Get the number of steps of the method
206 * (excluding the one being computed).
207 * @return number of steps of the method
208 * (excluding the one being computed)
209 * @deprecated as of 3.6, this method is not used anymore
210 */
211 @Deprecated
212 public int getNSteps() {
213 return c1.length;
214 }
215
216 /** Build the P matrix.
217 * <p>The P matrix general terms are shifted (j+1) (-i)<sup>j</sup> terms
218 * with i being the row number starting from 1 and j being the column
219 * number starting from 1:
220 * <pre>
221 * [ -2 3 -4 5 ... ]
222 * [ -4 12 -32 80 ... ]
223 * P = [ -6 27 -108 405 ... ]
224 * [ -8 48 -256 1280 ... ]
225 * [ ... ]
226 * </pre>
227 * @param rows number of rows of the matrix
228 * @return P matrix
229 */
230 private FieldDenseMatrix<BigFraction> buildP(final int rows) {
231 final FieldDenseMatrix<BigFraction> pData = FieldDenseMatrix.create(BigFractionField.get(),
232 rows, rows)
233 .fill(BigFraction.ZERO);
234
235 for (int i = 1; i <= rows; ++i) {
236 // build the P matrix elements from Taylor series formulas
237 final int factor = -i;
238 int aj = factor;
239 for (int j = 1; j <= rows; ++j) {
240 pData.set(i - 1, j - 1,
241 BigFraction.of(aj * (j + 1)));
242 aj *= factor;
243 }
244 }
245
246 return pData;
247 }
248
249 /** Initialize the high order scaled derivatives at step start.
250 * @param h step size to use for scaling
251 * @param t first steps times
252 * @param y first steps states
253 * @param yDot first steps derivatives
254 * @return Nordieck vector at start of first step (h<sup>2</sup>/2 y''<sub>n</sub>,
255 * h<sup>3</sup>/6 y'''<sub>n</sub> ... h<sup>k</sup>/k! y<sup>(k)</sup><sub>n</sub>)
256 */
257
258 public Array2DRowRealMatrix initializeHighOrderDerivatives(final double h, final double[] t,
259 final double[][] y,
260 final double[][] yDot) {
261
262 // using Taylor series with di = ti - t0, we get:
263 // y(ti) - y(t0) - di y'(t0) = di^2 / h^2 s2 + ... + di^k / h^k sk + O(h^k)
264 // y'(ti) - y'(t0) = 2 di / h^2 s2 + ... + k di^(k-1) / h^k sk + O(h^(k-1))
265 // we write these relations for i = 1 to i= 1+n/2 as a set of n + 2 linear
266 // equations depending on the Nordsieck vector [s2 ... sk rk], so s2 to sk correspond
267 // to the appropriately truncated Taylor expansion, and rk is the Taylor remainder.
268 // The goal is to have s2 to sk as accurate as possible considering the fact the sum is
269 // truncated and we don't want the error terms to be included in s2 ... sk, so we need
270 // to solve also for the remainder
271 final double[][] a = new double[c1.length + 1][c1.length + 1];
272 final double[][] b = new double[c1.length + 1][y[0].length];
273 final double[] y0 = y[0];
274 final double[] yDot0 = yDot[0];
275 for (int i = 1; i < y.length; ++i) {
276
277 final double di = t[i] - t[0];
278 final double ratio = di / h;
279 double dikM1Ohk = 1 / h;
280
281 // linear coefficients of equations
282 // y(ti) - y(t0) - di y'(t0) and y'(ti) - y'(t0)
283 final double[] aI = a[2 * i - 2];
284 final double[] aDotI = (2 * i - 1) < a.length ? a[2 * i - 1] : null;
285 for (int j = 0; j < aI.length; ++j) {
286 dikM1Ohk *= ratio;
287 aI[j] = di * dikM1Ohk;
288 if (aDotI != null) {
289 aDotI[j] = (j + 2) * dikM1Ohk;
290 }
291 }
292
293 // expected value of the previous equations
294 final double[] yI = y[i];
295 final double[] yDotI = yDot[i];
296 final double[] bI = b[2 * i - 2];
297 final double[] bDotI = (2 * i - 1) < b.length ? b[2 * i - 1] : null;
298 for (int j = 0; j < yI.length; ++j) {
299 bI[j] = yI[j] - y0[j] - di * yDot0[j];
300 if (bDotI != null) {
301 bDotI[j] = yDotI[j] - yDot0[j];
302 }
303 }
304 }
305
306 // solve the linear system to get the best estimate of the Nordsieck vector [s2 ... sk],
307 // with the additional terms s(k+1) and c grabbing the parts after the truncated Taylor expansion
308 final QRDecomposition decomposition = new QRDecomposition(new Array2DRowRealMatrix(a, false));
309 final RealMatrix x = decomposition.getSolver().solve(new Array2DRowRealMatrix(b, false));
310
311 // extract just the Nordsieck vector [s2 ... sk]
312 final Array2DRowRealMatrix truncatedX = new Array2DRowRealMatrix(x.getRowDimension() - 1, x.getColumnDimension());
313 for (int i = 0; i < truncatedX.getRowDimension(); ++i) {
314 for (int j = 0; j < truncatedX.getColumnDimension(); ++j) {
315 truncatedX.setEntry(i, j, x.getEntry(i, j));
316 }
317 }
318 return truncatedX;
319 }
320
321 /** Update the high order scaled derivatives for Adams integrators (phase 1).
322 * <p>The complete update of high order derivatives has a form similar to:
323 * <div style="white-space: pre"><code>
324 * 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>
325 * </code></div>
326 * this method computes the P<sup>-1</sup> A P r<sub>n</sub> part.
327 * @param highOrder high order scaled derivatives
328 * (h<sup>2</sup>/2 y'', ... h<sup>k</sup>/k! y(k))
329 * @return updated high order derivatives
330 * @see #updateHighOrderDerivativesPhase2(double[], double[], Array2DRowRealMatrix)
331 */
332 public Array2DRowRealMatrix updateHighOrderDerivativesPhase1(final Array2DRowRealMatrix highOrder) {
333 return update.multiply(highOrder);
334 }
335
336 /** Update the high order scaled derivatives Adams integrators (phase 2).
337 * <p>The complete update of high order derivatives has a form similar to:
338 * <div style="white-space: pre"><code>
339 * 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>
340 * </code></div>
341 * this method computes the (s<sub>1</sub>(n) - s<sub>1</sub>(n+1)) P<sup>-1</sup> u part.
342 * <p>Phase 1 of the update must already have been performed.</p>
343 * @param start first order scaled derivatives at step start
344 * @param end first order scaled derivatives at step end
345 * @param highOrder high order scaled derivatives, will be modified
346 * (h<sup>2</sup>/2 y'', ... h<sup>k</sup>/k! y(k))
347 * @see #updateHighOrderDerivativesPhase1(Array2DRowRealMatrix)
348 */
349 public void updateHighOrderDerivativesPhase2(final double[] start,
350 final double[] end,
351 final Array2DRowRealMatrix highOrder) {
352 final double[][] data = highOrder.getDataRef();
353 for (int i = 0; i < data.length; ++i) {
354 final double[] dataI = data[i];
355 final double c1I = c1[i];
356 for (int j = 0; j < dataI.length; ++j) {
357 dataI[j] += c1I * (start[j] - end[j]);
358 }
359 }
360 }
361 }