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.math3.ode.nonstiff; 019 020import java.util.Arrays; 021import java.util.HashMap; 022import java.util.Map; 023 024import org.apache.commons.math3.Field; 025import org.apache.commons.math3.RealFieldElement; 026import org.apache.commons.math3.linear.Array2DRowFieldMatrix; 027import org.apache.commons.math3.linear.ArrayFieldVector; 028import org.apache.commons.math3.linear.FieldDecompositionSolver; 029import org.apache.commons.math3.linear.FieldLUDecomposition; 030import org.apache.commons.math3.linear.FieldMatrix; 031import org.apache.commons.math3.util.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 * <pre> 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 * </pre></p> 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 * <pre> 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 * </pre> 054 * (we omit the k index in the notation for clarity).</p> 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 * <pre> 060 * r<sub>n</sub> = [ s<sub>2</sub>(n), s<sub>3</sub>(n) ... s<sub>k</sub>(n) ]<sup>T</sup> 061 * </pre> 062 * (here again we omit the k index in the notation for clarity) 063 * </p> 064 * 065 * <p>Taylor series formulas show that for any index offset i, s<sub>1</sub>(n-i) can be 066 * computed from s<sub>1</sub>(n), s<sub>2</sub>(n) ... s<sub>k</sub>(n), the formula being exact 067 * for degree k polynomials. 068 * <pre> 069 * 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) 070 * </pre> 071 * The previous formula can be used with several values for i to compute the transform between 072 * classical representation and Nordsieck vector at step end. The transform between r<sub>n</sub> 073 * and q<sub>n</sub> resulting from the Taylor series formulas above is: 074 * <pre> 075 * q<sub>n</sub> = s<sub>1</sub>(n) u + P r<sub>n</sub> 076 * </pre> 077 * where u is the [ 1 1 ... 1 ]<sup>T</sup> vector and P is the (k-1)×(k-1) matrix built 078 * with the (j+1) (-i)<sup>j</sup> terms with i being the row number starting from 1 and j being 079 * the column number starting from 1: 080 * <pre> 081 * [ -2 3 -4 5 ... ] 082 * [ -4 12 -32 80 ... ] 083 * P = [ -6 27 -108 405 ... ] 084 * [ -8 48 -256 1280 ... ] 085 * [ ... ] 086 * </pre></p> 087 * 088 * <p>Changing -i into +i in the formula above can be used to compute a similar transform between 089 * classical representation and Nordsieck vector at step start. The resulting matrix is simply 090 * the absolute value of matrix P.</p> 091 * 092 * <p>For {@link AdamsBashforthIntegrator Adams-Bashforth} method, the Nordsieck vector 093 * at step n+1 is computed from the Nordsieck vector at step n as follows: 094 * <ul> 095 * <li>y<sub>n+1</sub> = y<sub>n</sub> + s<sub>1</sub>(n) + u<sup>T</sup> r<sub>n</sub></li> 096 * <li>s<sub>1</sub>(n+1) = h f(t<sub>n+1</sub>, y<sub>n+1</sub>)</li> 097 * <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> 098 * </ul> 099 * where A is a rows shifting matrix (the lower left part is an identity matrix): 100 * <pre> 101 * [ 0 0 ... 0 0 | 0 ] 102 * [ ---------------+---] 103 * [ 1 0 ... 0 0 | 0 ] 104 * A = [ 0 1 ... 0 0 | 0 ] 105 * [ ... | 0 ] 106 * [ 0 0 ... 1 0 | 0 ] 107 * [ 0 0 ... 0 1 | 0 ] 108 * </pre></p> 109 * 110 * <p>For {@link AdamsMoultonIntegrator Adams-Moulton} method, the predicted Nordsieck vector 111 * at step n+1 is computed from the Nordsieck vector at step n as follows: 112 * <ul> 113 * <li>Y<sub>n+1</sub> = y<sub>n</sub> + s<sub>1</sub>(n) + u<sup>T</sup> r<sub>n</sub></li> 114 * <li>S<sub>1</sub>(n+1) = h f(t<sub>n+1</sub>, Y<sub>n+1</sub>)</li> 115 * <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> 116 * </ul> 117 * From this predicted vector, the corrected vector is computed as follows: 118 * <ul> 119 * <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> 120 * <li>s<sub>1</sub>(n+1) = h f(t<sub>n+1</sub>, y<sub>n+1</sub>)</li> 121 * <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> 122 * </ul> 123 * where the upper case Y<sub>n+1</sub>, S<sub>1</sub>(n+1) and R<sub>n+1</sub> represent the 124 * predicted states whereas the lower case y<sub>n+1</sub>, s<sub>n+1</sub> and r<sub>n+1</sub> 125 * represent the corrected states.</p> 126 * 127 * <p>We observe that both methods use similar update formulas. In both cases a P<sup>-1</sup>u 128 * vector and a P<sup>-1</sup> A P matrix are used that do not depend on the state, 129 * they only depend on k. This class handles these transformations.</p> 130 * 131 * @param <T> the type of the field elements 132 * @since 3.6 133 */ 134public class AdamsNordsieckFieldTransformer<T extends RealFieldElement<T>> { 135 136 /** Cache for already computed coefficients. */ 137 private static final Map<Integer, 138 Map<Field<? extends RealFieldElement<?>>, 139 AdamsNordsieckFieldTransformer<? extends RealFieldElement<?>>>> CACHE = 140 new HashMap<Integer, Map<Field<? extends RealFieldElement<?>>, 141 AdamsNordsieckFieldTransformer<? extends RealFieldElement<?>>>>(); 142 143 /** Field to which the time and state vector elements belong. */ 144 private final Field<T> field; 145 146 /** Update matrix for the higher order derivatives h<sup>2</sup>/2 y'', h<sup>3</sup>/6 y''' ... */ 147 private final Array2DRowFieldMatrix<T> update; 148 149 /** Update coefficients of the higher order derivatives wrt y'. */ 150 private final T[] c1; 151 152 /** Simple constructor. 153 * @param field field to which the time and state vector elements belong 154 * @param n number of steps of the multistep method 155 * (excluding the one being computed) 156 */ 157 private AdamsNordsieckFieldTransformer(final Field<T> field, final int n) { 158 159 this.field = field; 160 final int rows = n - 1; 161 162 // compute coefficients 163 FieldMatrix<T> bigP = buildP(rows); 164 FieldDecompositionSolver<T> pSolver = 165 new FieldLUDecomposition<T>(bigP).getSolver(); 166 167 T[] u = MathArrays.buildArray(field, rows); 168 Arrays.fill(u, field.getOne()); 169 c1 = pSolver.solve(new ArrayFieldVector<T>(u, false)).toArray(); 170 171 // update coefficients are computed by combining transform from 172 // Nordsieck to multistep, then shifting rows to represent step advance 173 // then applying inverse transform 174 T[][] shiftedP = bigP.getData(); 175 for (int i = shiftedP.length - 1; i > 0; --i) { 176 // shift rows 177 shiftedP[i] = shiftedP[i - 1]; 178 } 179 shiftedP[0] = MathArrays.buildArray(field, rows); 180 Arrays.fill(shiftedP[0], field.getZero()); 181 update = new Array2DRowFieldMatrix<T>(pSolver.solve(new Array2DRowFieldMatrix<T>(shiftedP, false)).getData()); 182 183 } 184 185 /** Get the Nordsieck transformer for a given field and number of steps. 186 * @param field field to which the time and state vector elements belong 187 * @param nSteps number of steps of the multistep method 188 * (excluding the one being computed) 189 * @return Nordsieck transformer for the specified field and number of steps 190 * @param <T> the type of the field elements 191 */ 192 @SuppressWarnings("unchecked") 193 public static <T extends RealFieldElement<T>> AdamsNordsieckFieldTransformer<T> 194 getInstance(final Field<T> field, final int nSteps) { 195 synchronized(CACHE) { 196 Map<Field<? extends RealFieldElement<?>>, 197 AdamsNordsieckFieldTransformer<? extends RealFieldElement<?>>> map = CACHE.get(nSteps); 198 if (map == null) { 199 map = new HashMap<Field<? extends RealFieldElement<?>>, 200 AdamsNordsieckFieldTransformer<? extends RealFieldElement<?>>>(); 201 CACHE.put(nSteps, map); 202 } 203 @SuppressWarnings("rawtypes") // use rawtype to avoid compilation problems with java 1.5 204 AdamsNordsieckFieldTransformer t = map.get(field); 205 if (t == null) { 206 t = new AdamsNordsieckFieldTransformer<T>(field, nSteps); 207 map.put(field, (AdamsNordsieckFieldTransformer<T>) t); 208 } 209 return (AdamsNordsieckFieldTransformer<T>) t; 210 211 } 212 } 213 214 /** Build the P matrix. 215 * <p>The P matrix general terms are shifted (j+1) (-i)<sup>j</sup> terms 216 * with i being the row number starting from 1 and j being the column 217 * number starting from 1: 218 * <pre> 219 * [ -2 3 -4 5 ... ] 220 * [ -4 12 -32 80 ... ] 221 * P = [ -6 27 -108 405 ... ] 222 * [ -8 48 -256 1280 ... ] 223 * [ ... ] 224 * </pre></p> 225 * @param rows number of rows of the matrix 226 * @return P matrix 227 */ 228 private FieldMatrix<T> buildP(final int rows) { 229 230 final T[][] pData = MathArrays.buildArray(field, rows, rows); 231 232 for (int i = 1; i <= pData.length; ++i) { 233 // build the P matrix elements from Taylor series formulas 234 final T[] pI = pData[i - 1]; 235 final int factor = -i; 236 T aj = field.getZero().add(factor); 237 for (int j = 1; j <= pI.length; ++j) { 238 pI[j - 1] = aj.multiply(j + 1); 239 aj = aj.multiply(factor); 240 } 241 } 242 243 return new Array2DRowFieldMatrix<T>(pData, false); 244 245 } 246 247 /** Initialize the high order scaled derivatives at step start. 248 * @param h step size to use for scaling 249 * @param t first steps times 250 * @param y first steps states 251 * @param yDot first steps derivatives 252 * @return Nordieck vector at start of first step (h<sup>2</sup>/2 y''<sub>n</sub>, 253 * h<sup>3</sup>/6 y'''<sub>n</sub> ... h<sup>k</sup>/k! y<sup>(k)</sup><sub>n</sub>) 254 */ 255 256 public Array2DRowFieldMatrix<T> initializeHighOrderDerivatives(final T h, final T[] t, 257 final T[][] y, 258 final T[][] yDot) { 259 260 // using Taylor series with di = ti - t0, we get: 261 // y(ti) - y(t0) - di y'(t0) = di^2 / h^2 s2 + ... + di^k / h^k sk + O(h^k) 262 // y'(ti) - y'(t0) = 2 di / h^2 s2 + ... + k di^(k-1) / h^k sk + O(h^(k-1)) 263 // we write these relations for i = 1 to i= 1+n/2 as a set of n + 2 linear 264 // equations depending on the Nordsieck vector [s2 ... sk rk], so s2 to sk correspond 265 // to the appropriately truncated Taylor expansion, and rk is the Taylor remainder. 266 // The goal is to have s2 to sk as accurate as possible considering the fact the sum is 267 // truncated and we don't want the error terms to be included in s2 ... sk, so we need 268 // to solve also for the remainder 269 final T[][] a = MathArrays.buildArray(field, c1.length + 1, c1.length + 1); 270 final T[][] b = MathArrays.buildArray(field, c1.length + 1, y[0].length); 271 final T[] y0 = y[0]; 272 final T[] yDot0 = yDot[0]; 273 for (int i = 1; i < y.length; ++i) { 274 275 final T di = t[i].subtract(t[0]); 276 final T ratio = di.divide(h); 277 T dikM1Ohk = h.reciprocal(); 278 279 // linear coefficients of equations 280 // y(ti) - y(t0) - di y'(t0) and y'(ti) - y'(t0) 281 final T[] aI = a[2 * i - 2]; 282 final T[] aDotI = (2 * i - 1) < a.length ? a[2 * i - 1] : null; 283 for (int j = 0; j < aI.length; ++j) { 284 dikM1Ohk = dikM1Ohk.multiply(ratio); 285 aI[j] = di.multiply(dikM1Ohk); 286 if (aDotI != null) { 287 aDotI[j] = dikM1Ohk.multiply(j + 2); 288 } 289 } 290 291 // expected value of the previous equations 292 final T[] yI = y[i]; 293 final T[] yDotI = yDot[i]; 294 final T[] bI = b[2 * i - 2]; 295 final T[] bDotI = (2 * i - 1) < b.length ? b[2 * i - 1] : null; 296 for (int j = 0; j < yI.length; ++j) { 297 bI[j] = yI[j].subtract(y0[j]).subtract(di.multiply(yDot0[j])); 298 if (bDotI != null) { 299 bDotI[j] = yDotI[j].subtract(yDot0[j]); 300 } 301 } 302 303 } 304 305 // solve the linear system to get the best estimate of the Nordsieck vector [s2 ... sk], 306 // with the additional terms s(k+1) and c grabbing the parts after the truncated Taylor expansion 307 final FieldLUDecomposition<T> decomposition = new FieldLUDecomposition<T>(new Array2DRowFieldMatrix<T>(a, false)); 308 final FieldMatrix<T> x = decomposition.getSolver().solve(new Array2DRowFieldMatrix<T>(b, false)); 309 310 // extract just the Nordsieck vector [s2 ... sk] 311 final Array2DRowFieldMatrix<T> truncatedX = 312 new Array2DRowFieldMatrix<T>(field, 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 322 /** Update the high order scaled derivatives for Adams integrators (phase 1). 323 * <p>The complete update of high order derivatives has a form similar to: 324 * <pre> 325 * 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> 326 * </pre> 327 * this method computes the P<sup>-1</sup> A P r<sub>n</sub> part.</p> 328 * @param highOrder high order scaled derivatives 329 * (h<sup>2</sup>/2 y'', ... h<sup>k</sup>/k! y(k)) 330 * @return updated high order derivatives 331 * @see #updateHighOrderDerivativesPhase2(RealFieldElement[], RealFieldElement[], Array2DRowFieldMatrix) 332 */ 333 public Array2DRowFieldMatrix<T> updateHighOrderDerivativesPhase1(final Array2DRowFieldMatrix<T> highOrder) { 334 return update.multiply(highOrder); 335 } 336 337 /** Update the high order scaled derivatives Adams integrators (phase 2). 338 * <p>The complete update of high order derivatives has a form similar to: 339 * <pre> 340 * 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> 341 * </pre> 342 * this method computes the (s<sub>1</sub>(n) - s<sub>1</sub>(n+1)) P<sup>-1</sup> u part.</p> 343 * <p>Phase 1 of the update must already have been performed.</p> 344 * @param start first order scaled derivatives at step start 345 * @param end first order scaled derivatives at step end 346 * @param highOrder high order scaled derivatives, will be modified 347 * (h<sup>2</sup>/2 y'', ... h<sup>k</sup>/k! y(k)) 348 * @see #updateHighOrderDerivativesPhase1(Array2DRowFieldMatrix) 349 */ 350 public void updateHighOrderDerivativesPhase2(final T[] start, 351 final T[] end, 352 final Array2DRowFieldMatrix<T> highOrder) { 353 final T[][] data = highOrder.getDataRef(); 354 for (int i = 0; i < data.length; ++i) { 355 final T[] dataI = data[i]; 356 final T c1I = c1[i]; 357 for (int j = 0; j < dataI.length; ++j) { 358 dataI[j] = dataI[j].add(c1I.multiply(start[j].subtract(end[j]))); 359 } 360 } 361 } 362 363}