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