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; 19 20 import java.io.ObjectInput; 21 import java.io.ObjectOutput; 22 import java.io.IOException; 23 24 /** This abstract class represents an interpolator over the last step 25 * during an ODE integration. 26 * 27 * <p>The various ODE integrators provide objects extending this class 28 * to the step handlers. The handlers can use these objects to 29 * retrieve the state vector at intermediate times between the 30 * previous and the current grid points (dense output).</p> 31 * 32 * @see FirstOrderIntegrator 33 * @see SecondOrderIntegrator 34 * @see StepHandler 35 * 36 * @version $Revision: 620312 $ $Date: 2008-02-10 12:28:59 -0700 (Sun, 10 Feb 2008) $ 37 * @since 1.2 38 * 39 */ 40 41 public abstract class AbstractStepInterpolator 42 implements StepInterpolator { 43 44 /** previous time */ 45 protected double previousTime; 46 47 /** current time */ 48 protected double currentTime; 49 50 /** current time step */ 51 protected double h; 52 53 /** current state */ 54 protected double[] currentState; 55 56 /** interpolated time */ 57 protected double interpolatedTime; 58 59 /** interpolated state */ 60 protected double[] interpolatedState; 61 62 /** indicate if the step has been finalized or not. */ 63 private boolean finalized; 64 65 /** integration direction. */ 66 private boolean forward; 67 68 /** Simple constructor. 69 * This constructor builds an instance that is not usable yet, the 70 * {@link #reinitialize} method should be called before using the 71 * instance in order to initialize the internal arrays. This 72 * constructor is used only in order to delay the initialization in 73 * some cases. As an example, the {@link 74 * EmbeddedRungeKuttaIntegrator} uses the prototyping design pattern 75 * to create the step interpolators by cloning an uninitialized 76 * model and latter initializing the copy. 77 */ 78 protected AbstractStepInterpolator() { 79 previousTime = Double.NaN; 80 currentTime = Double.NaN; 81 h = Double.NaN; 82 interpolatedTime = Double.NaN; 83 currentState = null; 84 interpolatedState = null; 85 finalized = false; 86 this.forward = true; 87 } 88 89 /** Simple constructor. 90 * @param y reference to the integrator array holding the state at 91 * the end of the step 92 * @param forward integration direction indicator 93 */ 94 protected AbstractStepInterpolator(double[] y, boolean forward) { 95 96 previousTime = Double.NaN; 97 currentTime = Double.NaN; 98 h = Double.NaN; 99 interpolatedTime = Double.NaN; 100 101 currentState = y; 102 interpolatedState = new double[y.length]; 103 104 finalized = false; 105 this.forward = forward; 106 107 } 108 109 /** Copy constructor. 110 111 * <p>The copied interpolator should have been finalized before the 112 * copy, otherwise the copy will not be able to perform correctly 113 * any derivative computation and will throw a {@link 114 * NullPointerException} later. Since we don't want this constructor 115 * to throw the exceptions finalization may involve and since we 116 * don't want this method to modify the state of the copied 117 * interpolator, finalization is <strong>not</strong> done 118 * automatically, it remains under user control.</p> 119 120 * <p>The copy is a deep copy: its arrays are separated from the 121 * original arrays of the instance.</p> 122 123 * @param interpolator interpolator to copy from. 124 125 */ 126 protected AbstractStepInterpolator(AbstractStepInterpolator interpolator) { 127 128 previousTime = interpolator.previousTime; 129 currentTime = interpolator.currentTime; 130 h = interpolator.h; 131 interpolatedTime = interpolator.interpolatedTime; 132 133 if (interpolator.currentState != null) { 134 currentState = (double[]) interpolator.currentState.clone(); 135 interpolatedState = (double[]) interpolator.interpolatedState.clone(); 136 } else { 137 currentState = null; 138 interpolatedState = null; 139 } 140 141 finalized = interpolator.finalized; 142 forward = interpolator.forward; 143 144 } 145 146 /** Reinitialize the instance 147 * @param y reference to the integrator array holding the state at 148 * the end of the step 149 * @param forward integration direction indicator 150 */ 151 protected void reinitialize(double[] y, boolean forward) { 152 153 previousTime = Double.NaN; 154 currentTime = Double.NaN; 155 h = Double.NaN; 156 interpolatedTime = Double.NaN; 157 158 currentState = y; 159 interpolatedState = new double[y.length]; 160 161 finalized = false; 162 this.forward = forward; 163 164 } 165 166 /** Copy the instance. 167 * <p>The copied instance is guaranteed to be independent from the 168 * original one. Both can be used with different settings for 169 * interpolated time without any side effect.</p> 170 * @return a deep copy of the instance, which can be used independently. 171 * @throws DerivativeException if this call induces an automatic 172 * step finalization that throws one 173 * @see #setInterpolatedTime(double) 174 */ 175 public StepInterpolator copy() throws DerivativeException { 176 177 // finalize the step before performing copy 178 finalizeStep(); 179 180 // create the new independent instance 181 return doCopy(); 182 183 } 184 185 /** Really copy the finalized instance. 186 * <p>This method is called by {@link #copy()} after the 187 * step has been finalized. It must perform a deep copy 188 * to have an new instance completely independent for the 189 * original instance. 190 * @return a copy of the finalized instance 191 */ 192 protected abstract StepInterpolator doCopy(); 193 194 /** Shift one step forward. 195 * Copy the current time into the previous time, hence preparing the 196 * interpolator for future calls to {@link #storeTime storeTime} 197 */ 198 public void shift() { 199 previousTime = currentTime; 200 } 201 202 /** Store the current step time. 203 * @param t current time 204 */ 205 public void storeTime(double t) { 206 207 currentTime = t; 208 h = currentTime - previousTime; 209 interpolatedTime = t; 210 System.arraycopy(currentState, 0, interpolatedState, 0, 211 currentState.length); 212 213 // the step is not finalized anymore 214 finalized = false; 215 216 } 217 218 /** 219 * Get the previous grid point time. 220 * @return previous grid point time 221 */ 222 public double getPreviousTime() { 223 return previousTime; 224 } 225 226 /** 227 * Get the current grid point time. 228 * @return current grid point time 229 */ 230 public double getCurrentTime() { 231 return currentTime; 232 } 233 234 /** 235 * Get the time of the interpolated point. 236 * If {@link #setInterpolatedTime} has not been called, it returns 237 * the current grid point time. 238 * @return interpolation point time 239 */ 240 public double getInterpolatedTime() { 241 return interpolatedTime; 242 } 243 244 /** 245 * Set the time of the interpolated point. 246 * <p>Setting the time outside of the current step is now allowed 247 * (it was not allowed up to version 5.4 of Mantissa), but should be 248 * used with care since the accuracy of the interpolator will 249 * probably be very poor far from this step. This allowance has been 250 * added to simplify implementation of search algorithms near the 251 * step endpoints.</p> 252 * @param time time of the interpolated point 253 * @throws DerivativeException if this call induces an automatic 254 * step finalization that throws one 255 */ 256 public void setInterpolatedTime(double time) 257 throws DerivativeException { 258 interpolatedTime = time; 259 double oneMinusThetaH = currentTime - interpolatedTime; 260 computeInterpolatedState((h - oneMinusThetaH) / h, oneMinusThetaH); 261 } 262 263 /** Check if the natural integration direction is forward. 264 * <p>This method provides the integration direction as specified by the 265 * integrator itself, it avoid some nasty problems in degenerated 266 * cases like null steps due to cancellation at step initialization, 267 * step control or switching function triggering.</p> 268 * @return true if the integration variable (time) increases during 269 * integration 270 */ 271 public boolean isForward() { 272 return forward; 273 } 274 275 /** Compute the state at the interpolated time. 276 * This is the main processing method that should be implemented by 277 * the derived classes to perform the interpolation. 278 * @param theta normalized interpolation abscissa within the step 279 * (theta is zero at the previous time step and one at the current time step) 280 * @param oneMinusThetaH time gap between the interpolated time and 281 * the current time 282 * @throws DerivativeException this exception is propagated to the caller if the 283 * underlying user function triggers one 284 */ 285 protected abstract void computeInterpolatedState(double theta, 286 double oneMinusThetaH) 287 throws DerivativeException; 288 289 /** 290 * Get the state vector of the interpolated point. 291 * @return state vector at time {@link #getInterpolatedTime} 292 */ 293 public double[] getInterpolatedState() { 294 return (double[]) interpolatedState.clone(); 295 } 296 297 298 /** 299 * Finalize the step. 300 301 * <p>Some embedded Runge-Kutta integrators need fewer functions 302 * evaluations than their counterpart step interpolators. These 303 * interpolators should perform the last evaluations they need by 304 * themselves only if they need them. This method triggers these 305 * extra evaluations. It can be called directly by the user step 306 * handler and it is called automatically if {@link 307 * #setInterpolatedTime} is called.</p> 308 309 * <p>Once this method has been called, <strong>no</strong> other 310 * evaluation will be performed on this step. If there is a need to 311 * have some side effects between the step handler and the 312 * differential equations (for example update some data in the 313 * equations once the step has been done), it is advised to call 314 * this method explicitly from the step handler before these side 315 * effects are set up. If the step handler induces no side effect, 316 * then this method can safely be ignored, it will be called 317 * transparently as needed.</p> 318 319 * <p><strong>Warning</strong>: since the step interpolator provided 320 * to the step handler as a parameter of the {@link 321 * StepHandler#handleStep handleStep} is valid only for the duration 322 * of the {@link StepHandler#handleStep handleStep} call, one cannot 323 * simply store a reference and reuse it later. One should first 324 * finalize the instance, then copy this finalized instance into a 325 * new object that can be kept.</p> 326 327 * <p>This method calls the protected <code>doFinalize</code> method 328 * if it has never been called during this step and set a flag 329 * indicating that it has been called once. It is the <code> 330 * doFinalize</code> method which should perform the evaluations. 331 * This wrapping prevents from calling <code>doFinalize</code> several 332 * times and hence evaluating the differential equations too often. 333 * Therefore, subclasses are not allowed not reimplement it, they 334 * should rather reimplement <code>doFinalize</code>.</p> 335 336 * @throws DerivativeException this exception is propagated to the 337 * caller if the underlying user function triggers one 338 339 */ 340 public final void finalizeStep() 341 throws DerivativeException { 342 if (! finalized) { 343 doFinalize(); 344 finalized = true; 345 } 346 } 347 348 /** 349 * Really finalize the step. 350 * The default implementation of this method does nothing. 351 * @throws DerivativeException this exception is propagated to the 352 * caller if the underlying user function triggers one 353 */ 354 protected void doFinalize() 355 throws DerivativeException { 356 } 357 358 /** Write the instance to an output channel. 359 * @param out output channel 360 * @exception IOException if the instance cannot be written 361 */ 362 public abstract void writeExternal(ObjectOutput out) 363 throws IOException; 364 365 /** Read the instance from an input channel. 366 * @param in input channel 367 * @exception IOException if the instance cannot be read 368 */ 369 public abstract void readExternal(ObjectInput in) 370 throws IOException; 371 372 /** Save the base state of the instance. 373 * This method performs step finalization if it has not been done 374 * before. 375 * @param out stream where to save the state 376 * @exception IOException in case of write error 377 */ 378 protected void writeBaseExternal(ObjectOutput out) 379 throws IOException { 380 381 out.writeInt(currentState.length); 382 out.writeDouble(previousTime); 383 out.writeDouble(currentTime); 384 out.writeDouble(h); 385 out.writeBoolean(forward); 386 387 for (int i = 0; i < currentState.length; ++i) { 388 out.writeDouble(currentState[i]); 389 } 390 391 out.writeDouble(interpolatedTime); 392 393 // we do not store the interpolated state, 394 // it will be recomputed as needed after reading 395 396 // finalize the step (and don't bother saving the now true flag) 397 try { 398 finalizeStep(); 399 } catch (DerivativeException e) { 400 throw new IOException(e.getMessage()); 401 } 402 403 } 404 405 /** Read the base state of the instance. 406 * This method does <strong>neither</strong> set the interpolated 407 * time nor state. It is up to the derived class to reset it 408 * properly calling the {@link #setInterpolatedTime} method later, 409 * once all rest of the object state has been set up properly. 410 * @param in stream where to read the state from 411 * @return interpolated time be set later by the caller 412 * @exception IOException in case of read error 413 */ 414 protected double readBaseExternal(ObjectInput in) 415 throws IOException { 416 417 int dimension = in.readInt(); 418 previousTime = in.readDouble(); 419 currentTime = in.readDouble(); 420 h = in.readDouble(); 421 forward = in.readBoolean(); 422 423 currentState = new double[dimension]; 424 for (int i = 0; i < currentState.length; ++i) { 425 currentState[i] = in.readDouble(); 426 } 427 428 // we do NOT handle the interpolated time and state here 429 interpolatedTime = Double.NaN; 430 interpolatedState = new double[dimension]; 431 432 finalized = true; 433 434 return in.readDouble(); 435 436 } 437 438 }