← Back to branch summary

~ubuntu-branches/ubuntu/oneiric/commons-math/oneiric

« back to all changes in this revision

Viewing changes to src/main/java/org/apache/commons/math/ode/nonstiff/AdamsNordsieckTransformer.java

  • Committer: Bazaar Package Importer
  • Author(s): Damien Raude-Morvan
  • Date: 2009-08-22 01:13:25 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20090822011325-hi4peq1ua5weguwn
Tags: 2.0-1
* New upstream release.
* Set Maintainer field to Debian Java Team
* Add myself as Uploaders
* Switch to Quilt patch system:
  - Refresh all patchs
  - Remove B-D on dpatch, Add B-D on quilt
  - Include patchsys-quilt.mk in debian/rules
* Bump Standards-Version to 3.8.3:
  - Add a README.source to describe patch system
* Maven POMs:
  - Add a Build-Depends-Indep dependency on maven-repo-helper
  - Use mh_installpom and mh_installjar to install the POM and the jar to the
    Maven repository
* Use default-jdk/jre:
  - Depends on java5-runtime-headless
  - Build-Depends on default-jdk
  - Use /usr/lib/jvm/default-java as JAVA_HOME
* Move api documentation to /usr/share/doc/libcommons-math-java/api
* Build-Depends on junit4 instead of junit

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/main/java/org/apache/commons/math/ode/nonstiff/AdamsNordsieckTransformer.java
 
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) + &sum;<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)&times;(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 ... &plusmn;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: 790374 $ $Date: 2009-07-01 16:57:20 -0400 (Wed, 01 Jul 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
}