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 018 package org.apache.commons.math.ode.nonstiff; 019 020 import java.util.Arrays; 021 import java.util.HashMap; 022 import java.util.Map; 023 024 import org.apache.commons.math.fraction.BigFraction; 025 import org.apache.commons.math.linear.Array2DRowFieldMatrix; 026 import org.apache.commons.math.linear.Array2DRowRealMatrix; 027 import org.apache.commons.math.linear.DefaultFieldMatrixChangingVisitor; 028 import org.apache.commons.math.linear.FieldDecompositionSolver; 029 import org.apache.commons.math.linear.FieldLUDecompositionImpl; 030 import org.apache.commons.math.linear.FieldMatrix; 031 import org.apache.commons.math.linear.MatrixUtils; 032 033 /** Transformer to Nordsieck vectors for Adams integrators. 034 * <p>This class i 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(k)<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</sub> j (-i)<sup>j-1</sup> s<sub>j</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 (-i)<sup>j-1</sup> terms: 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></p> 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></p> 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.</p> 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 * @version $Revision: 810196 $ $Date: 2009-09-01 21:47:46 +0200 (mar. 01 sept. 2009) $ 131 * @since 2.0 132 */ 133 public class AdamsNordsieckTransformer { 134 135 /** Cache for already computed coefficients. */ 136 private static final Map<Integer, AdamsNordsieckTransformer> CACHE = 137 new HashMap<Integer, AdamsNordsieckTransformer>(); 138 139 /** Initialization matrix for the higher order derivatives wrt y'', y''' ... */ 140 private final Array2DRowRealMatrix initialization; 141 142 /** Update matrix for the higher order derivatives h<sup>2</sup>/2y'', h<sup>3</sup>/6 y''' ... */ 143 private final Array2DRowRealMatrix update; 144 145 /** Update coefficients of the higher order derivatives wrt y'. */ 146 private final double[] c1; 147 148 /** Simple constructor. 149 * @param nSteps number of steps of the multistep method 150 * (excluding the one being computed) 151 */ 152 private AdamsNordsieckTransformer(final int nSteps) { 153 154 // compute exact coefficients 155 FieldMatrix<BigFraction> bigP = buildP(nSteps); 156 FieldDecompositionSolver<BigFraction> pSolver = 157 new FieldLUDecompositionImpl<BigFraction>(bigP).getSolver(); 158 159 BigFraction[] u = new BigFraction[nSteps]; 160 Arrays.fill(u, BigFraction.ONE); 161 BigFraction[] bigC1 = pSolver.solve(u); 162 163 // update coefficients are computed by combining transform from 164 // Nordsieck to multistep, then shifting rows to represent step advance 165 // then applying inverse transform 166 BigFraction[][] shiftedP = bigP.getData(); 167 for (int i = shiftedP.length - 1; i > 0; --i) { 168 // shift rows 169 shiftedP[i] = shiftedP[i - 1]; 170 } 171 shiftedP[0] = new BigFraction[nSteps]; 172 Arrays.fill(shiftedP[0], BigFraction.ZERO); 173 FieldMatrix<BigFraction> bigMSupdate = 174 pSolver.solve(new Array2DRowFieldMatrix<BigFraction>(shiftedP, false)); 175 176 // initialization coefficients, computed from a R matrix = abs(P) 177 bigP.walkInOptimizedOrder(new DefaultFieldMatrixChangingVisitor<BigFraction>(BigFraction.ZERO) { 178 /** {@inheritDoc} */ 179 @Override 180 public BigFraction visit(int row, int column, BigFraction value) { 181 return ((column & 0x1) == 0x1) ? value : value.negate(); 182 } 183 }); 184 FieldMatrix<BigFraction> bigRInverse = 185 new FieldLUDecompositionImpl<BigFraction>(bigP).getSolver().getInverse(); 186 187 // convert coefficients to double 188 initialization = MatrixUtils.bigFractionMatrixToRealMatrix(bigRInverse); 189 update = MatrixUtils.bigFractionMatrixToRealMatrix(bigMSupdate); 190 c1 = new double[nSteps]; 191 for (int i = 0; i < nSteps; ++i) { 192 c1[i] = bigC1[i].doubleValue(); 193 } 194 195 } 196 197 /** Get the Nordsieck transformer for a given number of steps. 198 * @param nSteps number of steps of the multistep method 199 * (excluding the one being computed) 200 * @return Nordsieck transformer for the specified number of steps 201 */ 202 public static AdamsNordsieckTransformer getInstance(final int nSteps) { 203 synchronized(CACHE) { 204 AdamsNordsieckTransformer t = CACHE.get(nSteps); 205 if (t == null) { 206 t = new AdamsNordsieckTransformer(nSteps); 207 CACHE.put(nSteps, t); 208 } 209 return t; 210 } 211 } 212 213 /** Get the number of steps of the method 214 * (excluding the one being computed). 215 * @return number of steps of the method 216 * (excluding the one being computed) 217 */ 218 public int getNSteps() { 219 return c1.length; 220 } 221 222 /** Build the P matrix. 223 * <p>The P matrix general terms are shifted j (-i)<sup>j-1</sup> terms: 224 * <pre> 225 * [ -2 3 -4 5 ... ] 226 * [ -4 12 -32 80 ... ] 227 * P = [ -6 27 -108 405 ... ] 228 * [ -8 48 -256 1280 ... ] 229 * [ ... ] 230 * </pre></p> 231 * @param nSteps number of steps of the multistep method 232 * (excluding the one being computed) 233 * @return P matrix 234 */ 235 private FieldMatrix<BigFraction> buildP(final int nSteps) { 236 237 final BigFraction[][] pData = new BigFraction[nSteps][nSteps]; 238 239 for (int i = 0; i < pData.length; ++i) { 240 // build the P matrix elements from Taylor series formulas 241 final BigFraction[] pI = pData[i]; 242 final int factor = -(i + 1); 243 int aj = factor; 244 for (int j = 0; j < pI.length; ++j) { 245 pI[j] = new BigFraction(aj * (j + 2)); 246 aj *= factor; 247 } 248 } 249 250 return new Array2DRowFieldMatrix<BigFraction>(pData, false); 251 252 } 253 254 /** Initialize the high order scaled derivatives at step start. 255 * @param first first scaled derivative at step start 256 * @param multistep scaled derivatives after step start (hy'1, ..., hy'k-1) 257 * will be modified 258 * @return high order derivatives at step start 259 */ 260 public Array2DRowRealMatrix initializeHighOrderDerivatives(final double[] first, 261 final double[][] multistep) { 262 for (int i = 0; i < multistep.length; ++i) { 263 final double[] msI = multistep[i]; 264 for (int j = 0; j < first.length; ++j) { 265 msI[j] -= first[j]; 266 } 267 } 268 return initialization.multiply(new Array2DRowRealMatrix(multistep, false)); 269 } 270 271 /** Update the high order scaled derivatives for Adams integrators (phase 1). 272 * <p>The complete update of high order derivatives has a form similar to: 273 * <pre> 274 * 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> 275 * </pre> 276 * this method computes the P<sup>-1</sup> A P r<sub>n</sub> part.</p> 277 * @param highOrder high order scaled derivatives 278 * (h<sup>2</sup>/2 y'', ... h<sup>k</sup>/k! y(k)) 279 * @return updated high order derivatives 280 * @see #updateHighOrderDerivativesPhase2(double[], double[], Array2DRowRealMatrix) 281 */ 282 public Array2DRowRealMatrix updateHighOrderDerivativesPhase1(final Array2DRowRealMatrix highOrder) { 283 return update.multiply(highOrder); 284 } 285 286 /** Update the high order scaled derivatives Adams integrators (phase 2). 287 * <p>The complete update of high order derivatives has a form similar to: 288 * <pre> 289 * 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> 290 * </pre> 291 * this method computes the (s<sub>1</sub>(n) - s<sub>1</sub>(n+1)) P<sup>-1</sup> u part.</p> 292 * <p>Phase 1 of the update must already have been performed.</p> 293 * @param start first order scaled derivatives at step start 294 * @param end first order scaled derivatives at step end 295 * @param highOrder high order scaled derivatives, will be modified 296 * (h<sup>2</sup>/2 y'', ... h<sup>k</sup>/k! y(k)) 297 * @see #updateHighOrderDerivativesPhase1(Array2DRowRealMatrix) 298 */ 299 public void updateHighOrderDerivativesPhase2(final double[] start, 300 final double[] end, 301 final Array2DRowRealMatrix highOrder) { 302 final double[][] data = highOrder.getDataRef(); 303 for (int i = 0; i < data.length; ++i) { 304 final double[] dataI = data[i]; 305 final double c1I = c1[i]; 306 for (int j = 0; j < dataI.length; ++j) { 307 dataI[j] += c1I * (start[j] - end[j]); 308 } 309 } 310 } 311 312 }