1dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond/* 2dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Licensed to the Apache Software Foundation (ASF) under one or more 3dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * contributor license agreements. See the NOTICE file distributed with 4dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * this work for additional information regarding copyright ownership. 5dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * The ASF licenses this file to You under the Apache License, Version 2.0 6dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * (the "License"); you may not use this file except in compliance with 7dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * the License. You may obtain a copy of the License at 8dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 9dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * http://www.apache.org/licenses/LICENSE-2.0 10dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 11dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Unless required by applicable law or agreed to in writing, software 12dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * distributed under the License is distributed on an "AS IS" BASIS, 13dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 14dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * See the License for the specific language governing permissions and 15dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * limitations under the License. 16dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 17dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 18dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondpackage org.apache.commons.math.ode.nonstiff; 19dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 20dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport java.io.IOException; 21dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport java.io.ObjectInput; 22dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport java.io.ObjectOutput; 23dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 24dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.ode.AbstractIntegrator; 25dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.ode.sampling.AbstractStepInterpolator; 26dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 27dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond/** This class represents an interpolator over the last step during an 28dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * ODE integration for Runge-Kutta and embedded Runge-Kutta integrators. 29dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 30dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @see RungeKuttaIntegrator 31dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @see EmbeddedRungeKuttaIntegrator 32dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 33dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @version $Revision: 811827 $ $Date: 2009-09-06 17:32:50 +0200 (dim. 06 sept. 2009) $ 34dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @since 1.2 35dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 36dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 37dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondabstract class RungeKuttaStepInterpolator 38dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond extends AbstractStepInterpolator { 39dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 40dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Slopes at the intermediate points */ 41dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond protected double[][] yDotK; 42dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 43dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Reference to the integrator. */ 44dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond protected AbstractIntegrator integrator; 45dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 46dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Simple constructor. 47dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * This constructor builds an instance that is not usable yet, the 48dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * {@link #reinitialize} method should be called before using the 49dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * instance in order to initialize the internal arrays. This 50dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * constructor is used only in order to delay the initialization in 51dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * some cases. The {@link RungeKuttaIntegrator} and {@link 52dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * EmbeddedRungeKuttaIntegrator} classes use the prototyping design 53dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * pattern to create the step interpolators by cloning an 54dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * uninitialized model and latter initializing the copy. 55dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 56dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond protected RungeKuttaStepInterpolator() { 57dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond super(); 58dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond yDotK = null; 59dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond integrator = null; 60dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 61dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 62dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Copy constructor. 63dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 64dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p>The copied interpolator should have been finalized before the 65dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * copy, otherwise the copy will not be able to perform correctly any 66dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * interpolation and will throw a {@link NullPointerException} 67dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * later. Since we don't want this constructor to throw the 68dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * exceptions finalization may involve and since we don't want this 69dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * method to modify the state of the copied interpolator, 70dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * finalization is <strong>not</strong> done automatically, it 71dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * remains under user control.</p> 72dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 73dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p>The copy is a deep copy: its arrays are separated from the 74dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * original arrays of the instance.</p> 75dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 76dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param interpolator interpolator to copy from. 77dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 78dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 79dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public RungeKuttaStepInterpolator(final RungeKuttaStepInterpolator interpolator) { 80dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 81dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond super(interpolator); 82dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 83dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (interpolator.currentState != null) { 84dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final int dimension = currentState.length; 85dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 86dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond yDotK = new double[interpolator.yDotK.length][]; 87dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int k = 0; k < interpolator.yDotK.length; ++k) { 88dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond yDotK[k] = new double[dimension]; 89dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond System.arraycopy(interpolator.yDotK[k], 0, 90dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond yDotK[k], 0, dimension); 91dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 92dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 93dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } else { 94dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond yDotK = null; 95dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 96dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 97dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // we cannot keep any reference to the equations in the copy 98dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // the interpolator should have been finalized before 99dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond integrator = null; 100dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 101dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 102dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 103dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Reinitialize the instance 104dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p>Some Runge-Kutta integrators need fewer functions evaluations 105dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * than their counterpart step interpolators. So the interpolator 106dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * should perform the last evaluations they need by themselves. The 107dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * {@link RungeKuttaIntegrator RungeKuttaIntegrator} and {@link 108dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * EmbeddedRungeKuttaIntegrator EmbeddedRungeKuttaIntegrator} 109dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * abstract classes call this method in order to let the step 110dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * interpolator perform the evaluations it needs. These evaluations 111dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * will be performed during the call to <code>doFinalize</code> if 112dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * any, i.e. only if the step handler either calls the {@link 113dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * AbstractStepInterpolator#finalizeStep finalizeStep} method or the 114dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * {@link AbstractStepInterpolator#getInterpolatedState 115dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * getInterpolatedState} method (for an interpolator which needs a 116dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * finalization) or if it clones the step interpolator.</p> 117dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param rkIntegrator integrator being used 118dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param y reference to the integrator array holding the state at 119dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * the end of the step 120dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param yDotArray reference to the integrator array holding all the 121dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * intermediate slopes 122dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param forward integration direction indicator 123dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 124dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public void reinitialize(final AbstractIntegrator rkIntegrator, 125dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[] y, final double[][] yDotArray, final boolean forward) { 126dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond reinitialize(y, forward); 127dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.yDotK = yDotArray; 128dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.integrator = rkIntegrator; 129dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 130dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 131dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** {@inheritDoc} */ 132dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond @Override 133dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public void writeExternal(final ObjectOutput out) 134dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throws IOException { 135dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 136dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // save the state of the base class 137dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond writeBaseExternal(out); 138dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 139dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // save the local attributes 140dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final int n = (currentState == null) ? -1 : currentState.length; 141dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final int kMax = (yDotK == null) ? -1 : yDotK.length; 142dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond out.writeInt(kMax); 143dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int k = 0; k < kMax; ++k) { 144dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int i = 0; i < n; ++i) { 145dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond out.writeDouble(yDotK[k][i]); 146dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 147dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 148dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 149dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // we do not save any reference to the equations 150dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 151dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 152dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 153dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** {@inheritDoc} */ 154dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond @Override 155dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public void readExternal(final ObjectInput in) 156dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throws IOException { 157dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 158dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // read the base class 159dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double t = readBaseExternal(in); 160dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 161dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // read the local attributes 162dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final int n = (currentState == null) ? -1 : currentState.length; 163dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final int kMax = in.readInt(); 164dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond yDotK = (kMax < 0) ? null : new double[kMax][]; 165dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int k = 0; k < kMax; ++k) { 166dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond yDotK[k] = (n < 0) ? null : new double[n]; 167dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int i = 0; i < n; ++i) { 168dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond yDotK[k][i] = in.readDouble(); 169dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 170dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 171dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 172dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond integrator = null; 173dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 174dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (currentState != null) { 175dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // we can now set the interpolated time and state 176dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond setInterpolatedTime(t); 177dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } else { 178dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond interpolatedTime = t; 179dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 180dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 181dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 182dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 183dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond} 184