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