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.sampling;
19
20 import java.io.IOException;
21 import java.io.ObjectInput;
22 import java.io.ObjectOutput;
23 import java.util.Arrays;
24
25 import org.apache.commons.math.linear.Array2DRowRealMatrix;
26 import org.apache.commons.math.ode.DerivativeException;
27
28 /**
29 * This class implements an interpolator for integrators using Nordsieck representation.
30 *
31 * <p>This interpolator computes dense output around the current point.
32 * The interpolation equation is based on Taylor series formulas.
33 *
34 * @see org.apache.commons.math.ode.nonstiff.AdamsBashforthIntegrator
35 * @see org.apache.commons.math.ode.nonstiff.AdamsMoultonIntegrator
36 * @version $Revision: 791723 $ $Date: 2009-07-07 03:07:23 -0400 (Tue, 07 Jul 2009) $
37 * @since 2.0
38 */
39
40 public class NordsieckStepInterpolator extends AbstractStepInterpolator {
41
42 /** Serializable version identifier */
43 private static final long serialVersionUID = -7179861704951334960L;
44
45 /** Step size used in the first scaled derivative and Nordsieck vector. */
46 private double scalingH;
47
48 /** Reference time for all arrays.
49 * <p>Sometimes, the reference time is the same as previousTime,
50 * sometimes it is the same as currentTime, so we use a separate
51 * field to avoid any confusion.
52 * </p>
53 */
54 private double referenceTime;
55
56 /** First scaled derivative. */
57 private double[] scaled;
58
59 /** Nordsieck vector. */
60 private Array2DRowRealMatrix nordsieck;
61
62 /** State variation. */
63 protected double[] stateVariation;
64
65 /** Simple constructor.
66 * This constructor builds an instance that is not usable yet, the
67 * {@link AbstractStepInterpolator#reinitialize} method should be called
68 * before using the instance in order to initialize the internal arrays. This
69 * constructor is used only in order to delay the initialization in
70 * some cases.
71 */
72 public NordsieckStepInterpolator() {
73 }
74
75 /** Copy constructor.
76 * @param interpolator interpolator to copy from. The copy is a deep
77 * copy: its arrays are separated from the original arrays of the
78 * instance
79 */
80 public NordsieckStepInterpolator(final NordsieckStepInterpolator interpolator) {
81 super(interpolator);
82 scalingH = interpolator.scalingH;
83 referenceTime = interpolator.referenceTime;
84 if (interpolator.scaled != null) {
85 scaled = interpolator.scaled.clone();
86 }
87 if (interpolator.nordsieck != null) {
88 nordsieck = new Array2DRowRealMatrix(interpolator.nordsieck.getDataRef(), true);
89 }
90 if (interpolator.stateVariation != null) {
91 stateVariation = interpolator.stateVariation.clone();
92 }
93 }
94
95 /** {@inheritDoc} */
96 @Override
97 protected StepInterpolator doCopy() {
98 return new NordsieckStepInterpolator(this);
99 }
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 referenceTime time at which all arrays are defined
118 * @param scalingH step size used in the scaled and nordsieck arrays
119 * @param scaled reference to the integrator array holding the first
120 * scaled derivative
121 * @param nordsieck reference to the integrator matrix holding the
122 * nordsieck vector
123 */
124 public void reinitialize(final double referenceTime, final double scalingH,
125 final double[] scaled, final Array2DRowRealMatrix nordsieck) {
126 this.referenceTime = referenceTime;
127 this.scalingH = scalingH;
128 this.scaled = scaled;
129 this.nordsieck = nordsieck;
130
131 // make sure the state and derivatives will depend on the new arrays
132 setInterpolatedTime(getInterpolatedTime());
133
134 }
135
136 /** Rescale the instance.
137 * <p>Since the scaled and Nordiseck arrays are shared with the caller,
138 * this method has the side effect of rescaling this arrays in the caller too.</p>
139 * @param scalingH new step size to use in the scaled and nordsieck arrays
140 */
141 public void rescale(final double scalingH) {
142
143 final double ratio = scalingH / this.scalingH;
144 for (int i = 0; i < scaled.length; ++i) {
145 scaled[i] *= ratio;
146 }
147
148 final double[][] nData = nordsieck.getDataRef();
149 double power = ratio;
150 for (int i = 0; i < nData.length; ++i) {
151 power *= ratio;
152 final double[] nDataI = nData[i];
153 for (int j = 0; j < nDataI.length; ++j) {
154 nDataI[j] *= power;
155 }
156 }
157
158 this.scalingH = scalingH;
159
160 }
161
162 /**
163 * Get the state vector variation from current to interpolated state.
164 * <p>This method is aimed at computing y(t<sub>interpolation</sub>)
165 * -y(t<sub>current</sub>) accurately by avoiding the cancellation errors
166 * that would occur if the subtraction were performed explicitly.</p>
167 * <p>The returned vector is a reference to a reused array, so
168 * it should not be modified and it should be copied if it needs
169 * to be preserved across several calls.</p>
170 * @return state vector at time {@link #getInterpolatedTime}
171 * @see #getInterpolatedDerivatives()
172 * @throws DerivativeException if this call induces an automatic
173 * step finalization that throws one
174 */
175 public double[] getInterpolatedStateVariation()
176 throws DerivativeException {
177 // compute and ignore interpolated state
178 // to make sure state variation is computed as a side effect
179 getInterpolatedState();
180 return stateVariation;
181 }
182
183 /** {@inheritDoc} */
184 @Override
185 protected void computeInterpolatedStateAndDerivatives(final double theta, final double oneMinusThetaH) {
186
187 final double x = interpolatedTime - referenceTime;
188 final double normalizedAbscissa = x / scalingH;
189
190 Arrays.fill(stateVariation, 0.0);
191 Arrays.fill(interpolatedDerivatives, 0.0);
192
193 // apply Taylor formula from high order to low order,
194 // for the sake of numerical accuracy
195 final double[][] nData = nordsieck.getDataRef();
196 for (int i = nData.length - 1; i >= 0; --i) {
197 final int order = i + 2;
198 final double[] nDataI = nData[i];
199 final double power = Math.pow(normalizedAbscissa, order);
200 for (int j = 0; j < nDataI.length; ++j) {
201 final double d = nDataI[j] * power;
202 stateVariation[j] += d;
203 interpolatedDerivatives[j] += order * d;
204 }
205 }
206
207 for (int j = 0; j < currentState.length; ++j) {
208 stateVariation[j] += scaled[j] * normalizedAbscissa;
209 interpolatedState[j] = currentState[j] + stateVariation[j];
210 interpolatedDerivatives[j] =
211 (interpolatedDerivatives[j] + scaled[j] * normalizedAbscissa) / x;
212 }
213
214 }
215
216 /** {@inheritDoc} */
217 @Override
218 public void writeExternal(final ObjectOutput out)
219 throws IOException {
220
221 // save the state of the base class
222 writeBaseExternal(out);
223
224 // save the local attributes
225 out.writeDouble(scalingH);
226 out.writeDouble(referenceTime);
227
228 final int n = (currentState == null) ? -1 : currentState.length;
229 if (scaled == null) {
230 out.writeBoolean(false);
231 } else {
232 out.writeBoolean(true);
233 for (int j = 0; j < n; ++j) {
234 out.writeDouble(scaled[j]);
235 }
236 }
237
238 if (nordsieck == null) {
239 out.writeBoolean(false);
240 } else {
241 out.writeBoolean(true);
242 out.writeObject(nordsieck);
243 }
244
245 // we don't save state variation, it will be recomputed
246
247 }
248
249 /** {@inheritDoc} */
250 @Override
251 public void readExternal(final ObjectInput in)
252 throws IOException, ClassNotFoundException {
253
254 // read the base class
255 final double t = readBaseExternal(in);
256
257 // read the local attributes
258 scalingH = in.readDouble();
259 referenceTime = in.readDouble();
260
261 final int n = (currentState == null) ? -1 : currentState.length;
262 final boolean hasScaled = in.readBoolean();
263 if (hasScaled) {
264 scaled = new double[n];
265 for (int j = 0; j < n; ++j) {
266 scaled[j] = in.readDouble();
267 }
268 } else {
269 scaled = null;
270 }
271
272 final boolean hasNordsieck = in.readBoolean();
273 if (hasNordsieck) {
274 nordsieck = (Array2DRowRealMatrix) in.readObject();
275 } else {
276 nordsieck = null;
277 }
278
279 if (hasScaled && hasNordsieck) {
280 // we can now set the interpolated time and state
281 stateVariation = new double[n];
282 setInterpolatedTime(t);
283 } else {
284 stateVariation = null;
285 }
286
287 }
288
289 }