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
18package org.apache.commons.math.ode.nonstiff;
19
20import org.apache.commons.math.ode.DerivativeException;
21import org.apache.commons.math.ode.AbstractIntegrator;
22import org.apache.commons.math.ode.sampling.StepInterpolator;
23
24/**
25 * This class represents an interpolator over the last step during an
26 * ODE integration for the 5(4) Dormand-Prince integrator.
27 *
28 * @see DormandPrince54Integrator
29 *
30 * @version $Revision: 1073158 $ $Date: 2011-02-21 22:46:52 +0100 (lun. 21 févr. 2011) $
31 * @since 1.2
32 */
33
34class DormandPrince54StepInterpolator
35  extends RungeKuttaStepInterpolator {
36
37    /** Last row of the Butcher-array internal weights, element 0. */
38    private static final double A70 =    35.0 /  384.0;
39
40    // element 1 is zero, so it is neither stored nor used
41
42    /** Last row of the Butcher-array internal weights, element 2. */
43    private static final double A72 =   500.0 / 1113.0;
44
45    /** Last row of the Butcher-array internal weights, element 3. */
46    private static final double A73 =   125.0 /  192.0;
47
48    /** Last row of the Butcher-array internal weights, element 4. */
49    private static final double A74 = -2187.0 / 6784.0;
50
51    /** Last row of the Butcher-array internal weights, element 5. */
52    private static final double A75 =    11.0 /   84.0;
53
54    /** Shampine (1986) Dense output, element 0. */
55    private static final double D0 =  -12715105075.0 /  11282082432.0;
56
57    // element 1 is zero, so it is neither stored nor used
58
59    /** Shampine (1986) Dense output, element 2. */
60    private static final double D2 =   87487479700.0 /  32700410799.0;
61
62    /** Shampine (1986) Dense output, element 3. */
63    private static final double D3 =  -10690763975.0 /   1880347072.0;
64
65    /** Shampine (1986) Dense output, element 4. */
66    private static final double D4 =  701980252875.0 / 199316789632.0;
67
68    /** Shampine (1986) Dense output, element 5. */
69    private static final double D5 =   -1453857185.0 /    822651844.0;
70
71    /** Shampine (1986) Dense output, element 6. */
72    private static final double D6 =      69997945.0 /     29380423.0;
73
74    /** Serializable version identifier */
75    private static final long serialVersionUID = 4104157279605906956L;
76
77    /** First vector for interpolation. */
78    private double[] v1;
79
80    /** Second vector for interpolation. */
81    private double[] v2;
82
83    /** Third vector for interpolation. */
84    private double[] v3;
85
86    /** Fourth vector for interpolation. */
87    private double[] v4;
88
89    /** Initialization indicator for the interpolation vectors. */
90    private boolean vectorsInitialized;
91
92  /** Simple constructor.
93   * This constructor builds an instance that is not usable yet, the
94   * {@link #reinitialize} method should be called before using the
95   * instance in order to initialize the internal arrays. This
96   * constructor is used only in order to delay the initialization in
97   * some cases. The {@link EmbeddedRungeKuttaIntegrator} uses the
98   * prototyping design pattern to create the step interpolators by
99   * cloning an uninitialized model and latter initializing the copy.
100   */
101  public DormandPrince54StepInterpolator() {
102    super();
103    v1 = null;
104    v2 = null;
105    v3 = null;
106    v4 = null;
107    vectorsInitialized = false;
108  }
109
110  /** Copy constructor.
111   * @param interpolator interpolator to copy from. The copy is a deep
112   * copy: its arrays are separated from the original arrays of the
113   * instance
114   */
115  public DormandPrince54StepInterpolator(final DormandPrince54StepInterpolator interpolator) {
116
117    super(interpolator);
118
119    if (interpolator.v1 == null) {
120
121      v1 = null;
122      v2 = null;
123      v3 = null;
124      v4 = null;
125      vectorsInitialized = false;
126
127    } else {
128
129      v1 = interpolator.v1.clone();
130      v2 = interpolator.v2.clone();
131      v3 = interpolator.v3.clone();
132      v4 = interpolator.v4.clone();
133      vectorsInitialized = interpolator.vectorsInitialized;
134
135    }
136
137  }
138
139  /** {@inheritDoc} */
140  @Override
141  protected StepInterpolator doCopy() {
142    return new DormandPrince54StepInterpolator(this);
143  }
144
145
146  /** {@inheritDoc} */
147  @Override
148  public void reinitialize(final AbstractIntegrator integrator,
149                           final double[] y, final double[][] yDotK, final boolean forward) {
150    super.reinitialize(integrator, y, yDotK, forward);
151    v1 = null;
152    v2 = null;
153    v3 = null;
154    v4 = null;
155    vectorsInitialized = false;
156  }
157
158  /** {@inheritDoc} */
159  @Override
160  public void storeTime(final double t) {
161    super.storeTime(t);
162    vectorsInitialized = false;
163  }
164
165  /** {@inheritDoc} */
166  @Override
167  protected void computeInterpolatedStateAndDerivatives(final double theta,
168                                          final double oneMinusThetaH)
169    throws DerivativeException {
170
171    if (! vectorsInitialized) {
172
173      if (v1 == null) {
174        v1 = new double[interpolatedState.length];
175        v2 = new double[interpolatedState.length];
176        v3 = new double[interpolatedState.length];
177        v4 = new double[interpolatedState.length];
178      }
179
180      // no step finalization is needed for this interpolator
181
182      // we need to compute the interpolation vectors for this time step
183      for (int i = 0; i < interpolatedState.length; ++i) {
184          final double yDot0 = yDotK[0][i];
185          final double yDot2 = yDotK[2][i];
186          final double yDot3 = yDotK[3][i];
187          final double yDot4 = yDotK[4][i];
188          final double yDot5 = yDotK[5][i];
189          final double yDot6 = yDotK[6][i];
190          v1[i] = A70 * yDot0 + A72 * yDot2 + A73 * yDot3 + A74 * yDot4 + A75 * yDot5;
191          v2[i] = yDot0 - v1[i];
192          v3[i] = v1[i] - v2[i] - yDot6;
193          v4[i] = D0 * yDot0 + D2 * yDot2 + D3 * yDot3 + D4 * yDot4 + D5 * yDot5 + D6 * yDot6;
194      }
195
196      vectorsInitialized = true;
197
198    }
199
200    // interpolate
201    final double eta = 1 - theta;
202    final double twoTheta = 2 * theta;
203    final double dot2 = 1 - twoTheta;
204    final double dot3 = theta * (2 - 3 * theta);
205    final double dot4 = twoTheta * (1 + theta * (twoTheta - 3));
206    for (int i = 0; i < interpolatedState.length; ++i) {
207      interpolatedState[i] =
208          currentState[i] - oneMinusThetaH * (v1[i] - theta * (v2[i] + theta * (v3[i] + eta * v4[i])));
209      interpolatedDerivatives[i] = v1[i] + dot2 * v2[i] + dot3 * v3[i] + dot4 * v4[i];
210    }
211
212  }
213
214}
215