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