View Javadoc

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