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 org.apache.commons.math.ode.DerivativeException; 21dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.ode.AbstractIntegrator; 22dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.ode.sampling.StepInterpolator; 23dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 24dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond/** 25dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * This class represents an interpolator over the last step during an 26dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * ODE integration for the 5(4) Dormand-Prince integrator. 27dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 28dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @see DormandPrince54Integrator 29dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 30dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @version $Revision: 1073158 $ $Date: 2011-02-21 22:46:52 +0100 (lun. 21 févr. 2011) $ 31dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @since 1.2 32dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 33dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 34dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondclass DormandPrince54StepInterpolator 35dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond extends RungeKuttaStepInterpolator { 36dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 37dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Last row of the Butcher-array internal weights, element 0. */ 38dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private static final double A70 = 35.0 / 384.0; 39dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 40dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // element 1 is zero, so it is neither stored nor used 41dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 42dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Last row of the Butcher-array internal weights, element 2. */ 43dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private static final double A72 = 500.0 / 1113.0; 44dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 45dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Last row of the Butcher-array internal weights, element 3. */ 46dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private static final double A73 = 125.0 / 192.0; 47dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 48dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Last row of the Butcher-array internal weights, element 4. */ 49dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private static final double A74 = -2187.0 / 6784.0; 50dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 51dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Last row of the Butcher-array internal weights, element 5. */ 52dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private static final double A75 = 11.0 / 84.0; 53dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 54dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Shampine (1986) Dense output, element 0. */ 55dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private static final double D0 = -12715105075.0 / 11282082432.0; 56dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 57dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // element 1 is zero, so it is neither stored nor used 58dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 59dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Shampine (1986) Dense output, element 2. */ 60dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private static final double D2 = 87487479700.0 / 32700410799.0; 61dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 62dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Shampine (1986) Dense output, element 3. */ 63dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private static final double D3 = -10690763975.0 / 1880347072.0; 64dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 65dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Shampine (1986) Dense output, element 4. */ 66dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private static final double D4 = 701980252875.0 / 199316789632.0; 67dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 68dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Shampine (1986) Dense output, element 5. */ 69dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private static final double D5 = -1453857185.0 / 822651844.0; 70dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 71dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Shampine (1986) Dense output, element 6. */ 72dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private static final double D6 = 69997945.0 / 29380423.0; 73dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 74dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Serializable version identifier */ 75dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private static final long serialVersionUID = 4104157279605906956L; 76dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 77dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** First vector for interpolation. */ 78dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private double[] v1; 79dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 80dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Second vector for interpolation. */ 81dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private double[] v2; 82dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 83dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Third vector for interpolation. */ 84dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private double[] v3; 85dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 86dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Fourth vector for interpolation. */ 87dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private double[] v4; 88dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 89dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Initialization indicator for the interpolation vectors. */ 90dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private boolean vectorsInitialized; 91dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 92dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Simple constructor. 93dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * This constructor builds an instance that is not usable yet, the 94dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * {@link #reinitialize} method should be called before using the 95dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * instance in order to initialize the internal arrays. This 96dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * constructor is used only in order to delay the initialization in 97dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * some cases. The {@link EmbeddedRungeKuttaIntegrator} uses the 98dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * prototyping design pattern to create the step interpolators by 99dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * cloning an uninitialized model and latter initializing the copy. 100dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 101dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public DormandPrince54StepInterpolator() { 102dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond super(); 103dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond v1 = null; 104dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond v2 = null; 105dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond v3 = null; 106dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond v4 = null; 107dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond vectorsInitialized = false; 108dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 109dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 110dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Copy constructor. 111dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param interpolator interpolator to copy from. The copy is a deep 112dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * copy: its arrays are separated from the original arrays of the 113dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * instance 114dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 115dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public DormandPrince54StepInterpolator(final DormandPrince54StepInterpolator interpolator) { 116dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 117dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond super(interpolator); 118dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 119dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (interpolator.v1 == null) { 120dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 121dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond v1 = null; 122dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond v2 = null; 123dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond v3 = null; 124dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond v4 = null; 125dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond vectorsInitialized = false; 126dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 127dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } else { 128dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 129dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond v1 = interpolator.v1.clone(); 130dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond v2 = interpolator.v2.clone(); 131dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond v3 = interpolator.v3.clone(); 132dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond v4 = interpolator.v4.clone(); 133dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond vectorsInitialized = interpolator.vectorsInitialized; 134dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 135dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 136dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 137dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 138dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 139dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** {@inheritDoc} */ 140dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond @Override 141dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond protected StepInterpolator doCopy() { 142dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return new DormandPrince54StepInterpolator(this); 143dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 144dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 145dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 146dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** {@inheritDoc} */ 147dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond @Override 148dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public void reinitialize(final AbstractIntegrator integrator, 149dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[] y, final double[][] yDotK, final boolean forward) { 150dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond super.reinitialize(integrator, y, yDotK, forward); 151dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond v1 = null; 152dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond v2 = null; 153dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond v3 = null; 154dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond v4 = null; 155dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond vectorsInitialized = false; 156dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 157dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 158dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** {@inheritDoc} */ 159dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond @Override 160dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public void storeTime(final double t) { 161dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond super.storeTime(t); 162dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond vectorsInitialized = false; 163dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 164dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 165dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** {@inheritDoc} */ 166dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond @Override 167dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond protected void computeInterpolatedStateAndDerivatives(final double theta, 168dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double oneMinusThetaH) 169dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throws DerivativeException { 170dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 171dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (! vectorsInitialized) { 172dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 173dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (v1 == null) { 174dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond v1 = new double[interpolatedState.length]; 175dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond v2 = new double[interpolatedState.length]; 176dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond v3 = new double[interpolatedState.length]; 177dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond v4 = new double[interpolatedState.length]; 178dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 179dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 180dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // no step finalization is needed for this interpolator 181dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 182dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // we need to compute the interpolation vectors for this time step 183dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int i = 0; i < interpolatedState.length; ++i) { 184dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double yDot0 = yDotK[0][i]; 185dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double yDot2 = yDotK[2][i]; 186dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double yDot3 = yDotK[3][i]; 187dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double yDot4 = yDotK[4][i]; 188dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double yDot5 = yDotK[5][i]; 189dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double yDot6 = yDotK[6][i]; 190dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond v1[i] = A70 * yDot0 + A72 * yDot2 + A73 * yDot3 + A74 * yDot4 + A75 * yDot5; 191dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond v2[i] = yDot0 - v1[i]; 192dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond v3[i] = v1[i] - v2[i] - yDot6; 193dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond v4[i] = D0 * yDot0 + D2 * yDot2 + D3 * yDot3 + D4 * yDot4 + D5 * yDot5 + D6 * yDot6; 194dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 195dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 196dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond vectorsInitialized = true; 197dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 198dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 199dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 200dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // interpolate 201dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double eta = 1 - theta; 202dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double twoTheta = 2 * theta; 203dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double dot2 = 1 - twoTheta; 204dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double dot3 = theta * (2 - 3 * theta); 205dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double dot4 = twoTheta * (1 + theta * (twoTheta - 3)); 206dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int i = 0; i < interpolatedState.length; ++i) { 207dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond interpolatedState[i] = 208dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond currentState[i] - oneMinusThetaH * (v1[i] - theta * (v2[i] + theta * (v3[i] + eta * v4[i]))); 209dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond interpolatedDerivatives[i] = v1[i] + dot2 * v2[i] + dot3 * v3[i] + dot4 * v4[i]; 210dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 211dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 212dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 213dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 214dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond} 215