← Back to branch summary

~ubuntu-branches/ubuntu/maverick/commons-math/maverick

« back to all changes in this revision

Viewing changes to src/main/java/org/apache/commons/math/ode/nonstiff/RungeKuttaIntegrator.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/RungeKuttaIntegrator.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
 
 
21
import org.apache.commons.math.ode.AbstractIntegrator;
 
22
import org.apache.commons.math.ode.DerivativeException;
 
23
import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
 
24
import org.apache.commons.math.ode.IntegratorException;
 
25
import org.apache.commons.math.ode.events.CombinedEventsManager;
 
26
import org.apache.commons.math.ode.sampling.AbstractStepInterpolator;
 
27
import org.apache.commons.math.ode.sampling.DummyStepInterpolator;
 
28
import org.apache.commons.math.ode.sampling.StepHandler;
 
29
 
 
30
/**
 
31
 * This class implements the common part of all fixed step Runge-Kutta
 
32
 * integrators for Ordinary Differential Equations.
 
33
 *
 
34
 * <p>These methods are explicit Runge-Kutta methods, their Butcher
 
35
 * arrays are as follows :
 
36
 * <pre>
 
37
 *    0  |
 
38
 *   c2  | a21
 
39
 *   c3  | a31  a32
 
40
 *   ... |        ...
 
41
 *   cs  | as1  as2  ...  ass-1
 
42
 *       |--------------------------
 
43
 *       |  b1   b2  ...   bs-1  bs
 
44
 * </pre>
 
45
 * </p>
 
46
 *
 
47
 * @see EulerIntegrator
 
48
 * @see ClassicalRungeKuttaIntegrator
 
49
 * @see GillIntegrator
 
50
 * @see MidpointIntegrator
 
51
 * @version $Revision: 785473 $ $Date: 2009-06-17 00:02:35 -0400 (Wed, 17 Jun 2009) $
 
52
 * @since 1.2
 
53
 */
 
54
 
 
55
public abstract class RungeKuttaIntegrator extends AbstractIntegrator {
 
56
 
 
57
  /** Simple constructor.
 
58
   * Build a Runge-Kutta integrator with the given
 
59
   * step. The default step handler does nothing.
 
60
   * @param name name of the method
 
61
   * @param c time steps from Butcher array (without the first zero)
 
62
   * @param a internal weights from Butcher array (without the first empty row)
 
63
   * @param b propagation weights for the high order method from Butcher array
 
64
   * @param prototype prototype of the step interpolator to use
 
65
   * @param step integration step
 
66
   */
 
67
  protected RungeKuttaIntegrator(final String name,
 
68
                                 final double[] c, final double[][] a, final double[] b,
 
69
                                 final RungeKuttaStepInterpolator prototype,
 
70
                                 final double step) {
 
71
    super(name);
 
72
    this.c          = c;
 
73
    this.a          = a;
 
74
    this.b          = b;
 
75
    this.prototype  = prototype;
 
76
    this.step       = Math.abs(step);
 
77
  }
 
78
 
 
79
  /** {@inheritDoc} */
 
80
  public double integrate(final FirstOrderDifferentialEquations equations,
 
81
                          final double t0, final double[] y0,
 
82
                          final double t, final double[] y)
 
83
  throws DerivativeException, IntegratorException {
 
84
 
 
85
    sanityChecks(equations, t0, y0, t, y);
 
86
    setEquations(equations);
 
87
    resetEvaluations();
 
88
    final boolean forward = (t > t0);
 
89
 
 
90
    // create some internal working arrays
 
91
    final int stages = c.length + 1;
 
92
    if (y != y0) {
 
93
      System.arraycopy(y0, 0, y, 0, y0.length);
 
94
    }
 
95
    final double[][] yDotK = new double[stages][];
 
96
    for (int i = 0; i < stages; ++i) {
 
97
      yDotK [i] = new double[y0.length];
 
98
    }
 
99
    final double[] yTmp = new double[y0.length];
 
100
 
 
101
    // set up an interpolator sharing the integrator arrays
 
102
    AbstractStepInterpolator interpolator;
 
103
    if (requiresDenseOutput() || (! eventsHandlersManager.isEmpty())) {
 
104
      final RungeKuttaStepInterpolator rki = (RungeKuttaStepInterpolator) prototype.copy();
 
105
      rki.reinitialize(this, yTmp, yDotK, forward);
 
106
      interpolator = rki;
 
107
    } else {
 
108
      interpolator = new DummyStepInterpolator(yTmp, forward);
 
109
    }
 
110
    interpolator.storeTime(t0);
 
111
 
 
112
    // set up integration control objects
 
113
    stepStart = t0;
 
114
    stepSize  = forward ? step : -step;
 
115
    for (StepHandler handler : stepHandlers) {
 
116
        handler.reset();
 
117
    }
 
118
    CombinedEventsManager manager = addEndTimeChecker(t0, t, eventsHandlersManager);
 
119
    boolean lastStep = false;
 
120
 
 
121
    // main integration loop
 
122
    while (!lastStep) {
 
123
 
 
124
      interpolator.shift();
 
125
 
 
126
      for (boolean loop = true; loop;) {
 
127
 
 
128
        // first stage
 
129
        computeDerivatives(stepStart, y, yDotK[0]);
 
130
 
 
131
        // next stages
 
132
        for (int k = 1; k < stages; ++k) {
 
133
 
 
134
          for (int j = 0; j < y0.length; ++j) {
 
135
            double sum = a[k-1][0] * yDotK[0][j];
 
136
            for (int l = 1; l < k; ++l) {
 
137
              sum += a[k-1][l] * yDotK[l][j];
 
138
            }
 
139
            yTmp[j] = y[j] + stepSize * sum;
 
140
          }
 
141
 
 
142
          computeDerivatives(stepStart + c[k-1] * stepSize, yTmp, yDotK[k]);
 
143
 
 
144
        }
 
145
 
 
146
        // estimate the state at the end of the step
 
147
        for (int j = 0; j < y0.length; ++j) {
 
148
          double sum    = b[0] * yDotK[0][j];
 
149
          for (int l = 1; l < stages; ++l) {
 
150
            sum    += b[l] * yDotK[l][j];
 
151
          }
 
152
          yTmp[j] = y[j] + stepSize * sum;
 
153
        }
 
154
 
 
155
        // discrete events handling
 
156
        interpolator.storeTime(stepStart + stepSize);
 
157
        if (manager.evaluateStep(interpolator)) {
 
158
            final double dt = manager.getEventTime() - stepStart;
 
159
            if (Math.abs(dt) <= Math.ulp(stepStart)) {
 
160
                // rejecting the step would lead to a too small next step, we accept it
 
161
                loop = false;
 
162
            } else {
 
163
                // reject the step to match exactly the next switch time
 
164
                stepSize = dt;
 
165
            }
 
166
        } else {
 
167
          loop = false;
 
168
        }
 
169
 
 
170
      }
 
171
 
 
172
      // the step has been accepted
 
173
      final double nextStep = stepStart + stepSize;
 
174
      System.arraycopy(yTmp, 0, y, 0, y0.length);
 
175
      manager.stepAccepted(nextStep, y);
 
176
      lastStep = manager.stop();
 
177
 
 
178
      // provide the step data to the step handler
 
179
      interpolator.storeTime(nextStep);
 
180
      for (StepHandler handler : stepHandlers) {
 
181
          handler.handleStep(interpolator, lastStep);
 
182
      }
 
183
      stepStart = nextStep;
 
184
 
 
185
      if (manager.reset(stepStart, y) && ! lastStep) {
 
186
        // some events handler has triggered changes that
 
187
        // invalidate the derivatives, we need to recompute them
 
188
        computeDerivatives(stepStart, y, yDotK[0]);
 
189
      }
 
190
 
 
191
      // make sure step size is set to default before next step
 
192
      stepSize = forward ? step : -step;
 
193
 
 
194
    }
 
195
 
 
196
    final double stopTime = stepStart;
 
197
    stepStart = Double.NaN;
 
198
    stepSize  = Double.NaN;
 
199
    return stopTime;
 
200
 
 
201
  }
 
202
 
 
203
  /** Time steps from Butcher array (without the first zero). */
 
204
  private double[] c;
 
205
 
 
206
  /** Internal weights from Butcher array (without the first empty row). */
 
207
  private double[][] a;
 
208
 
 
209
  /** External weights for the high order method from Butcher array. */
 
210
  private double[] b;
 
211
 
 
212
  /** Prototype of the step interpolator. */
 
213
  private RungeKuttaStepInterpolator prototype;
 
214
                                         
 
215
  /** Integration step. */
 
216
  private double step;
 
217
 
 
218
}