001 /*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements. See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License. You may obtain a copy of the License at
008 *
009 * http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017
018 package org.apache.commons.math3.ode.sampling;
019
020 import java.io.IOException;
021 import java.io.ObjectInput;
022 import java.io.ObjectOutput;
023 import java.util.Arrays;
024
025 import org.apache.commons.math3.exception.MaxCountExceededException;
026 import org.apache.commons.math3.linear.Array2DRowRealMatrix;
027 import org.apache.commons.math3.ode.EquationsMapper;
028 import org.apache.commons.math3.util.FastMath;
029
030 /**
031 * This class implements an interpolator for integrators using Nordsieck representation.
032 *
033 * <p>This interpolator computes dense output around the current point.
034 * The interpolation equation is based on Taylor series formulas.
035 *
036 * @see org.apache.commons.math3.ode.nonstiff.AdamsBashforthIntegrator
037 * @see org.apache.commons.math3.ode.nonstiff.AdamsMoultonIntegrator
038 * @version $Id: NordsieckStepInterpolator.java 1416643 2012-12-03 19:37:14Z tn $
039 * @since 2.0
040 */
041
042 public class NordsieckStepInterpolator extends AbstractStepInterpolator {
043
044 /** Serializable version identifier */
045 private static final long serialVersionUID = -7179861704951334960L;
046
047 /** State variation. */
048 protected double[] stateVariation;
049
050 /** Step size used in the first scaled derivative and Nordsieck vector. */
051 private double scalingH;
052
053 /** Reference time for all arrays.
054 * <p>Sometimes, the reference time is the same as previousTime,
055 * sometimes it is the same as currentTime, so we use a separate
056 * field to avoid any confusion.
057 * </p>
058 */
059 private double referenceTime;
060
061 /** First scaled derivative. */
062 private double[] scaled;
063
064 /** Nordsieck vector. */
065 private Array2DRowRealMatrix nordsieck;
066
067 /** Simple constructor.
068 * This constructor builds an instance that is not usable yet, the
069 * {@link AbstractStepInterpolator#reinitialize} method should be called
070 * before using the instance in order to initialize the internal arrays. This
071 * constructor is used only in order to delay the initialization in
072 * some cases.
073 */
074 public NordsieckStepInterpolator() {
075 }
076
077 /** Copy constructor.
078 * @param interpolator interpolator to copy from. The copy is a deep
079 * copy: its arrays are separated from the original arrays of the
080 * instance
081 */
082 public NordsieckStepInterpolator(final NordsieckStepInterpolator interpolator) {
083 super(interpolator);
084 scalingH = interpolator.scalingH;
085 referenceTime = interpolator.referenceTime;
086 if (interpolator.scaled != null) {
087 scaled = interpolator.scaled.clone();
088 }
089 if (interpolator.nordsieck != null) {
090 nordsieck = new Array2DRowRealMatrix(interpolator.nordsieck.getDataRef(), true);
091 }
092 if (interpolator.stateVariation != null) {
093 stateVariation = interpolator.stateVariation.clone();
094 }
095 }
096
097 /** {@inheritDoc} */
098 @Override
099 protected StepInterpolator doCopy() {
100 return new NordsieckStepInterpolator(this);
101 }
102
103 /** Reinitialize the instance.
104 * <p>Beware that all arrays <em>must</em> be references to integrator
105 * arrays, in order to ensure proper update without copy.</p>
106 * @param y reference to the integrator array holding the state at
107 * the end of the step
108 * @param forward integration direction indicator
109 * @param primaryMapper equations mapper for the primary equations set
110 * @param secondaryMappers equations mappers for the secondary equations sets
111 */
112 @Override
113 public void reinitialize(final double[] y, final boolean forward,
114 final EquationsMapper primaryMapper,
115 final EquationsMapper[] secondaryMappers) {
116 super.reinitialize(y, forward, primaryMapper, secondaryMappers);
117 stateVariation = new double[y.length];
118 }
119
120 /** Reinitialize the instance.
121 * <p>Beware that all arrays <em>must</em> be references to integrator
122 * arrays, in order to ensure proper update without copy.</p>
123 * @param time time at which all arrays are defined
124 * @param stepSize step size used in the scaled and nordsieck arrays
125 * @param scaledDerivative reference to the integrator array holding the first
126 * scaled derivative
127 * @param nordsieckVector reference to the integrator matrix holding the
128 * nordsieck vector
129 */
130 public void reinitialize(final double time, final double stepSize,
131 final double[] scaledDerivative,
132 final Array2DRowRealMatrix nordsieckVector) {
133 this.referenceTime = time;
134 this.scalingH = stepSize;
135 this.scaled = scaledDerivative;
136 this.nordsieck = nordsieckVector;
137
138 // make sure the state and derivatives will depend on the new arrays
139 setInterpolatedTime(getInterpolatedTime());
140
141 }
142
143 /** Rescale the instance.
144 * <p>Since the scaled and Nordiseck arrays are shared with the caller,
145 * this method has the side effect of rescaling this arrays in the caller too.</p>
146 * @param stepSize new step size to use in the scaled and nordsieck arrays
147 */
148 public void rescale(final double stepSize) {
149
150 final double ratio = stepSize / scalingH;
151 for (int i = 0; i < scaled.length; ++i) {
152 scaled[i] *= ratio;
153 }
154
155 final double[][] nData = nordsieck.getDataRef();
156 double power = ratio;
157 for (int i = 0; i < nData.length; ++i) {
158 power *= ratio;
159 final double[] nDataI = nData[i];
160 for (int j = 0; j < nDataI.length; ++j) {
161 nDataI[j] *= power;
162 }
163 }
164
165 scalingH = stepSize;
166
167 }
168
169 /**
170 * Get the state vector variation from current to interpolated state.
171 * <p>This method is aimed at computing y(t<sub>interpolation</sub>)
172 * -y(t<sub>current</sub>) accurately by avoiding the cancellation errors
173 * that would occur if the subtraction were performed explicitly.</p>
174 * <p>The returned vector is a reference to a reused array, so
175 * it should not be modified and it should be copied if it needs
176 * to be preserved across several calls.</p>
177 * @return state vector at time {@link #getInterpolatedTime}
178 * @see #getInterpolatedDerivatives()
179 * @exception MaxCountExceededException if the number of functions evaluations is exceeded
180 */
181 public double[] getInterpolatedStateVariation() throws MaxCountExceededException {
182 // compute and ignore interpolated state
183 // to make sure state variation is computed as a side effect
184 getInterpolatedState();
185 return stateVariation;
186 }
187
188 /** {@inheritDoc} */
189 @Override
190 protected void computeInterpolatedStateAndDerivatives(final double theta, final double oneMinusThetaH) {
191
192 final double x = interpolatedTime - referenceTime;
193 final double normalizedAbscissa = x / scalingH;
194
195 Arrays.fill(stateVariation, 0.0);
196 Arrays.fill(interpolatedDerivatives, 0.0);
197
198 // apply Taylor formula from high order to low order,
199 // for the sake of numerical accuracy
200 final double[][] nData = nordsieck.getDataRef();
201 for (int i = nData.length - 1; i >= 0; --i) {
202 final int order = i + 2;
203 final double[] nDataI = nData[i];
204 final double power = FastMath.pow(normalizedAbscissa, order);
205 for (int j = 0; j < nDataI.length; ++j) {
206 final double d = nDataI[j] * power;
207 stateVariation[j] += d;
208 interpolatedDerivatives[j] += order * d;
209 }
210 }
211
212 for (int j = 0; j < currentState.length; ++j) {
213 stateVariation[j] += scaled[j] * normalizedAbscissa;
214 interpolatedState[j] = currentState[j] + stateVariation[j];
215 interpolatedDerivatives[j] =
216 (interpolatedDerivatives[j] + scaled[j] * normalizedAbscissa) / x;
217 }
218
219 }
220
221 /** {@inheritDoc} */
222 @Override
223 public void writeExternal(final ObjectOutput out)
224 throws IOException {
225
226 // save the state of the base class
227 writeBaseExternal(out);
228
229 // save the local attributes
230 out.writeDouble(scalingH);
231 out.writeDouble(referenceTime);
232
233 final int n = (currentState == null) ? -1 : currentState.length;
234 if (scaled == null) {
235 out.writeBoolean(false);
236 } else {
237 out.writeBoolean(true);
238 for (int j = 0; j < n; ++j) {
239 out.writeDouble(scaled[j]);
240 }
241 }
242
243 if (nordsieck == null) {
244 out.writeBoolean(false);
245 } else {
246 out.writeBoolean(true);
247 out.writeObject(nordsieck);
248 }
249
250 // we don't save state variation, it will be recomputed
251
252 }
253
254 /** {@inheritDoc} */
255 @Override
256 public void readExternal(final ObjectInput in)
257 throws IOException, ClassNotFoundException {
258
259 // read the base class
260 final double t = readBaseExternal(in);
261
262 // read the local attributes
263 scalingH = in.readDouble();
264 referenceTime = in.readDouble();
265
266 final int n = (currentState == null) ? -1 : currentState.length;
267 final boolean hasScaled = in.readBoolean();
268 if (hasScaled) {
269 scaled = new double[n];
270 for (int j = 0; j < n; ++j) {
271 scaled[j] = in.readDouble();
272 }
273 } else {
274 scaled = null;
275 }
276
277 final boolean hasNordsieck = in.readBoolean();
278 if (hasNordsieck) {
279 nordsieck = (Array2DRowRealMatrix) in.readObject();
280 } else {
281 nordsieck = null;
282 }
283
284 if (hasScaled && hasNordsieck) {
285 // we can now set the interpolated time and state
286 stateVariation = new double[n];
287 setInterpolatedTime(t);
288 } else {
289 stateVariation = null;
290 }
291
292 }
293
294 }