Raymond | dee0849 | 2015-04-02 10:43:13 -0700 | [diff] [blame] | 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.math.ode.nonstiff; |
| 19 | |
| 20 | import java.util.Arrays; |
| 21 | import java.util.HashMap; |
| 22 | import java.util.Map; |
| 23 | |
| 24 | import org.apache.commons.math.fraction.BigFraction; |
| 25 | import org.apache.commons.math.linear.Array2DRowFieldMatrix; |
| 26 | import org.apache.commons.math.linear.Array2DRowRealMatrix; |
| 27 | import org.apache.commons.math.linear.DefaultFieldMatrixChangingVisitor; |
| 28 | import org.apache.commons.math.linear.FieldDecompositionSolver; |
| 29 | import org.apache.commons.math.linear.FieldLUDecompositionImpl; |
| 30 | import org.apache.commons.math.linear.FieldMatrix; |
| 31 | import org.apache.commons.math.linear.MatrixUtils; |
| 32 | |
| 33 | /** Transformer to Nordsieck vectors for Adams integrators. |
| 34 | * <p>This class i 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 | * <pre> |
| 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(k)<sub>n</sub> for k<sup>th</sup> derivative |
| 46 | * </pre></p> |
| 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 | * <pre> |
| 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 | * </pre> |
| 54 | * (we omit the k index in the notation for clarity).</p> |
| 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 | * <pre> |
| 60 | * r<sub>n</sub> = [ s<sub>2</sub>(n), s<sub>3</sub>(n) ... s<sub>k</sub>(n) ]<sup>T</sup> |
| 61 | * </pre> |
| 62 | * (here again we omit the k index in the notation for clarity) |
| 63 | * </p> |
| 64 | * |
| 65 | * <p>Taylor series formulas show that for any index offset i, s<sub>1</sub>(n-i) can be |
| 66 | * computed from s<sub>1</sub>(n), s<sub>2</sub>(n) ... s<sub>k</sub>(n), the formula being exact |
| 67 | * for degree k polynomials. |
| 68 | * <pre> |
| 69 | * 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) |
| 70 | * </pre> |
| 71 | * The previous formula can be used with several values for i to compute the transform between |
| 72 | * classical representation and Nordsieck vector at step end. The transform between r<sub>n</sub> |
| 73 | * and q<sub>n</sub> resulting from the Taylor series formulas above is: |
| 74 | * <pre> |
| 75 | * q<sub>n</sub> = s<sub>1</sub>(n) u + P r<sub>n</sub> |
| 76 | * </pre> |
| 77 | * where u is the [ 1 1 ... 1 ]<sup>T</sup> vector and P is the (k-1)×(k-1) matrix built |
| 78 | * with the j (-i)<sup>j-1</sup> terms: |
| 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></p> |
| 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></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 | } |