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: 791723 $ $Date: 2009-07-07 03:07:23 -0400 (Tue, 07 Jul 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 /** Step size used in the first scaled derivative and Nordsieck vector. */ 046 private double scalingH; 047 048 /** Reference time for all arrays. 049 * <p>Sometimes, the reference time is the same as previousTime, 050 * sometimes it is the same as currentTime, so we use a separate 051 * field to avoid any confusion. 052 * </p> 053 */ 054 private double referenceTime; 055 056 /** First scaled derivative. */ 057 private double[] scaled; 058 059 /** Nordsieck vector. */ 060 private Array2DRowRealMatrix nordsieck; 061 062 /** State variation. */ 063 protected double[] stateVariation; 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 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 }