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