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