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