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