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 018package org.apache.commons.math4.legacy.ode.nonstiff; 019 020import java.util.Arrays; 021import java.util.HashMap; 022import java.util.Map; 023 024import org.apache.commons.math4.legacy.core.Field; 025import org.apache.commons.math4.legacy.core.RealFieldElement; 026import org.apache.commons.math4.legacy.linear.Array2DRowFieldMatrix; 027import org.apache.commons.math4.legacy.linear.ArrayFieldVector; 028import org.apache.commons.math4.legacy.linear.FieldDecompositionSolver; 029import org.apache.commons.math4.legacy.linear.FieldLUDecomposition; 030import org.apache.commons.math4.legacy.linear.FieldMatrix; 031import org.apache.commons.math4.legacy.core.MathArrays; 032 033/** Transformer to Nordsieck vectors for Adams integrators. 034 * <p>This class is used by {@link AdamsBashforthIntegrator Adams-Bashforth} and 035 * {@link AdamsMoultonIntegrator Adams-Moulton} integrators to convert between 036 * classical representation with several previous first derivatives and Nordsieck 037 * representation with higher order scaled derivatives.</p> 038 * 039 * <p>We define scaled derivatives s<sub>i</sub>(n) at step n as: 040 * <div style="white-space: pre"><code> 041 * s<sub>1</sub>(n) = h y'<sub>n</sub> for first derivative 042 * s<sub>2</sub>(n) = h<sup>2</sup>/2 y''<sub>n</sub> for second derivative 043 * s<sub>3</sub>(n) = h<sup>3</sup>/6 y'''<sub>n</sub> for third derivative 044 * ... 045 * s<sub>k</sub>(n) = h<sup>k</sup>/k! y<sup>(k)</sup><sub>n</sub> for k<sup>th</sup> derivative 046 * </code></div> 047 * 048 * <p>With the previous definition, the classical representation of multistep methods 049 * uses first derivatives only, i.e. it handles y<sub>n</sub>, s<sub>1</sub>(n) and 050 * q<sub>n</sub> where q<sub>n</sub> is defined as: 051 * <div style="white-space: pre"><code> 052 * 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> 053 * </code></div> 054 * (we omit the k index in the notation for clarity). 055 * 056 * <p>Another possible representation uses the Nordsieck vector with 057 * higher degrees scaled derivatives all taken at the same step, i.e it handles y<sub>n</sub>, 058 * s<sub>1</sub>(n) and r<sub>n</sub>) where r<sub>n</sub> is defined as: 059 * <div style="white-space: pre"><code> 060 * r<sub>n</sub> = [ s<sub>2</sub>(n), s<sub>3</sub>(n) ... s<sub>k</sub>(n) ]<sup>T</sup> 061 * </code></div> 062 * (here again we omit the k index in the notation for clarity) 063 * 064 * <p>Taylor series formulas show that for any index offset i, s<sub>1</sub>(n-i) can be 065 * computed from s<sub>1</sub>(n), s<sub>2</sub>(n) ... s<sub>k</sub>(n), the formula being exact 066 * for degree k polynomials. 067 * <div style="white-space: pre"><code> 068 * 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) 069 * </code></div> 070 * The previous formula can be used with several values for i to compute the transform between 071 * classical representation and Nordsieck vector at step end. The transform between r<sub>n</sub> 072 * and q<sub>n</sub> resulting from the Taylor series formulas above is: 073 * <div style="white-space: pre"><code> 074 * q<sub>n</sub> = s<sub>1</sub>(n) u + P r<sub>n</sub> 075 * </code></div> 076 * where u is the [ 1 1 ... 1 ]<sup>T</sup> vector and P is the (k-1)×(k-1) matrix built 077 * with the (j+1) (-i)<sup>j</sup> terms with i being the row number starting from 1 and j being 078 * the column number starting from 1: 079 * <pre> 080 * [ -2 3 -4 5 ... ] 081 * [ -4 12 -32 80 ... ] 082 * P = [ -6 27 -108 405 ... ] 083 * [ -8 48 -256 1280 ... ] 084 * [ ... ] 085 * </pre> 086 * 087 * <p>Changing -i into +i in the formula above can be used to compute a similar transform between 088 * classical representation and Nordsieck vector at step start. The resulting matrix is simply 089 * the absolute value of matrix P.</p> 090 * 091 * <p>For {@link AdamsBashforthIntegrator Adams-Bashforth} method, the Nordsieck vector 092 * at step n+1 is computed from the Nordsieck vector at step n as follows: 093 * <ul> 094 * <li>y<sub>n+1</sub> = y<sub>n</sub> + s<sub>1</sub>(n) + u<sup>T</sup> r<sub>n</sub></li> 095 * <li>s<sub>1</sub>(n+1) = h f(t<sub>n+1</sub>, y<sub>n+1</sub>)</li> 096 * <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> 097 * </ul> 098 * where A is a rows shifting matrix (the lower left part is an identity matrix): 099 * <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 */ 133public 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}