001/* Explicit variable step size Runge-Kutta 4(5) ODE solver. 002 003 Copyright (c) 2004-2014 The Regents of the University of California. 004 All rights reserved. 005 Permission is hereby granted, without written agreement and without 006 license or royalty fees, to use, copy, modify, and distribute this 007 software and its documentation for any purpose, provided that the above 008 copyright notice and the following two paragraphs appear in all copies 009 of this software. 010 011 IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA BE LIABLE TO ANY PARTY 012 FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES 013 ARISING OUT OF THE USE OF THIS SOFTWARE AND ITS DOCUMENTATION, EVEN IF 014 THE UNIVERSITY OF CALIFORNIA HAS BEEN ADVISED OF THE POSSIBILITY OF 015 SUCH DAMAGE. 016 017 THE UNIVERSITY OF CALIFORNIA SPECIFICALLY DISCLAIMS ANY WARRANTIES, 018 INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF 019 MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE SOFTWARE 020 PROVIDED HEREUNDER IS ON AN "AS IS" BASIS, AND THE UNIVERSITY OF 021 CALIFORNIA HAS NO OBLIGATION TO PROVIDE MAINTENANCE, SUPPORT, UPDATES, 022 ENHANCEMENTS, OR MODIFICATIONS. 023 024 PT_COPYRIGHT_VERSION_2 025 COPYRIGHTENDKEY 026 027 */ 028package ptolemy.domains.continuous.kernel.solver; 029 030import ptolemy.domains.continuous.kernel.ContinuousIntegrator; 031import ptolemy.domains.continuous.kernel.ContinuousODESolver; 032import ptolemy.kernel.util.IllegalActionException; 033import ptolemy.kernel.util.InvalidStateException; 034 035/////////////////////////////////////////////////////////////////// 036//// ExplicitRK45Solver 037 038/** 039 This class implements a fourth-order Runge-Kutta ODE solving method. 040 The algorithm was introduced in "A Variable Order Runge-Kutta 041 Method for Initial Value Problems with Rapidly Varying Right-Hand Sides" 042 by J. R. Cash and Alan H. Karp, ACM Transactions on Mathematical Software, 043 vol 16, pp. 201-222, 1990. For completeness, a brief explanation of the 044 algorithm is explained below. 045 <p> 046 For an ODE of the form: 047 <pre> 048 dx(t)/dt = f(x(t), t), x(0) = x0 049 </pre> 050 it does the following: 051 <pre> 052 K0 = f(x(n), tn); 053 K1 = f(x(n) + 0.2*K0*h, tn + 0.2*h); 054 K2 = f(x(n) + (3.0/40*K0 + 9.0/40*K1)*h, tn + 0.3*h); 055 K3 = f(x(n) + (0.3*K0 - 0.9*K1 + 1.2*K2)*h, tn + 0.6*h); 056 K4 = f(x(n) + (-11/54*K0 + 5.0/2*K1 -70/27*K2 + 35/27*K3)*h, tn + 1.0*h); 057 K5 = f(x(n) + (1631/55296*K0 + 175/512*K1 + 575/13824*K2 + 3544275/110592*K3 + 253/4096*K4)*h, tn + 7/8*h); 058 x(n+1) = x(n)+(37/378*K0 + 250/621*K2 + 125.0/594*K3 + 512.0/1771*K5)*h; 059 </pre>, 060 and error control: 061 <pre> 062 LTE = [(37.0/378 - 2825.0/27648)*K0 + (250.0/621 - 18575.0/48384)*K2 + 063 (125.0/594 - 13525.0/55296)*K3 + (0.0 - 277.0/14336)*K4 + 064 (512.0/1771 - 0.25)*K5]*h. 065 </pre> 066 <P> 067 If the LTE is less than the error tolerance, then this step size h is 068 considered successful, and the next integration step size h' is predicted as: 069 <pre> 070 h' = h * Math.pow((ErrorTolerance/LTE), 1.0/5.0) 071 </pre> 072 This is a fourth order method, but uses a fifth order procedure to estimate 073 the local truncation error. 074 <p> 075 It takes 6 steps for this solver to resolve a state with an integration 076 step size. 077 078 @author Haiyang Zheng, Edward A. Lee 079 @version $Id$ 080 @since Ptolemy II 6.0 081 @Pt.ProposedRating Green (hyzheng) 082 @Pt.AcceptedRating Green (hyzheng) 083 */ 084public class ExplicitRK45Solver extends ContinuousODESolver { 085 086 /////////////////////////////////////////////////////////////////// 087 //// public methods //// 088 089 /** Return the number of time increments plus one (to store the 090 * truncation error). 091 * @return The number of time increments plus one. 092 */ 093 @Override 094 public final int getIntegratorAuxVariableCount() { 095 // Allow one for the truncation error 096 return _TIME_INCREMENTS.length + 1; 097 } 098 099 /** Fire the given integrator. This method performs the ODE solving 100 * algorithm described in the class comment. 101 * @param integrator The integrator of that calls this method. 102 * @exception IllegalActionException If there is no director, or can not 103 * read input, or can not send output. 104 */ 105 @Override 106 public void integratorIntegrate(ContinuousIntegrator integrator) 107 throws IllegalActionException { 108 double xn = integrator.getState(); 109 double outputValue; 110 double h = _director.getCurrentStepSize(); 111 double[] k = integrator.getAuxVariables(); 112 integrator.setAuxVariables(_roundCount, integrator.getDerivative()); 113 114 switch (_roundCount) { 115 case 0: 116 outputValue = xn + h * k[0] * _B[0][0]; 117 break; 118 119 case 1: 120 outputValue = xn + h * (k[0] * _B[1][0] + k[1] * _B[1][1]); 121 break; 122 123 case 2: 124 outputValue = xn 125 + h * (k[0] * _B[2][0] + k[1] * _B[2][1] + k[2] * _B[2][2]); 126 break; 127 128 case 3: 129 outputValue = xn + h * (k[0] * _B[3][0] + k[1] * _B[3][1] 130 + k[2] * _B[3][2] + k[3] * _B[3][3]); 131 break; 132 133 case 4: 134 outputValue = xn + h * (k[0] * _B[4][0] + k[1] * _B[4][1] 135 + k[2] * _B[4][2] + k[3] * _B[4][3] + k[4] * _B[4][4]); 136 break; 137 138 case 5: 139 outputValue = xn + h * (k[0] * _B[5][0] + k[1] * _B[5][1] 140 + k[2] * _B[5][2] + k[3] * _B[5][3] + k[4] * _B[5][4] 141 + k[5] * _B[5][5]); 142 break; 143 144 case 6: 145 outputValue = integrator.getTentativeState(); 146 return; 147 148 default: 149 throw new InvalidStateException("Execution sequence out of range."); 150 } 151 152 integrator.setTentativeState(outputValue); 153 } 154 155 /** Return true if the integration is accurate for the given 156 * integrator. This estimates the local truncation error for that 157 * integrator and compare it with the error tolerance. 158 * 159 * @param integrator The integrator of that calls this method. 160 * @return True if the integration is successful. 161 */ 162 @Override 163 public boolean integratorIsAccurate(ContinuousIntegrator integrator) { 164 double tolerance = _director.getErrorTolerance(); 165 double h = _director.getCurrentStepSize(); 166 double[] k = integrator.getAuxVariables(); 167 double error = h * Math.abs(k[0] * _E[0] + k[1] * _E[1] + k[2] * _E[2] 168 + k[3] * _E[3] + k[4] * _E[4] + k[5] * _E[5]); 169 170 integrator.setAuxVariables(_ERROR_INDEX, error); 171 172 if (_isDebugging()) { 173 _debug("Integrator: " + integrator.getName() 174 + " local truncation error = " + error); 175 } 176 if (error < tolerance) { 177 if (_isDebugging()) { 178 _debug("Integrator: " + integrator.getName() 179 + " report a success."); 180 } 181 return true; 182 } else { 183 if (_isDebugging()) { 184 _debug("Integrator: " + integrator.getName() 185 + " reports a failure."); 186 } 187 return false; 188 } 189 } 190 191 /** Predict the next step size for the integrators executed under this 192 * solver. This uses the algorithm in the class comments 193 * to predict the next step size based on the current estimation 194 * of the local truncation error. 195 * 196 * @param integrator The integrator that calls this method. 197 * @return The next step size suggested by the given integrator. 198 */ 199 @Override 200 public double integratorSuggestedStepSize(ContinuousIntegrator integrator) { 201 double error = integrator.getAuxVariables()[_ERROR_INDEX]; 202 double h = _director.getCurrentStepSize(); 203 double tolerance = _director.getErrorTolerance(); 204 double newh = 5.0 * h; 205 206 if (error > tolerance) { 207 newh = h * Math.pow(tolerance / error, 1.0 / _ORDER); 208 } 209 210 if (_isDebugging()) { 211 _debug("integrator: " + integrator.getName() 212 + " suggests next step size = " + newh); 213 } 214 215 return newh; 216 } 217 218 /////////////////////////////////////////////////////////////////// 219 //// protected methods //// 220 221 /** Return the current round. 222 * @return The current round. 223 */ 224 @Override 225 protected int _getRound() { 226 return _roundCount; 227 } 228 229 /** Get the current round factor. If the 230 * step is finished, then return 1.0. 231 */ 232 @Override 233 protected final double _getRoundTimeIncrement() { 234 return _TIME_INCREMENTS[_roundCount]; 235 } 236 237 /** Return true if the current integration step is finished. 238 * This method will return true if _incrementRound() has been 239 * called 6 or more times since _reset(). 240 * @see #_reset() 241 */ 242 @Override 243 protected final boolean _isStepFinished() { 244 return _roundCount >= _TIME_INCREMENTS.length; 245 } 246 247 /** Reset the solver, indicating to it that we are starting an 248 * integration step. This method resets the round counter. 249 */ 250 @Override 251 protected final void _reset() { 252 _roundCount = 0; 253 } 254 255 /** Set the round for the next integration step. 256 * @param round The round for the next integration step. 257 */ 258 @Override 259 protected void _setRound(int round) { 260 _roundCount = round; 261 } 262 263 /////////////////////////////////////////////////////////////////// 264 //// private variables //// 265 266 /** The ratio of time increments within one integration step. */ 267 protected static final double[] _TIME_INCREMENTS = { 0.2, 0.3, 0.6, 1.0, 268 0.875, 1.0, 1.0 }; 269 270 /** B coefficients */ 271 private static final double[][] _B = { { 0.2 }, { 3.0 / 40, 9.0 / 40 }, 272 { 0.3, -0.9, 1.2 }, { -11.0 / 54, 5.0 / 2, -70.0 / 27, 35.0 / 27 }, 273 { 1631.0 / 55296, 175.0 / 512, 575.0 / 13824, 44275.0 / 110592, 274 253.0 / 4096 }, 275 { 37.0 / 378, 0.0, 250.0 / 621, 125.0 / 594, 0.0, 512.0 / 1771 } }; 276 277 /** E coefficients */ 278 private static final double[] _E = { 37.0 / 378 - 2825.0 / 27648, 0.0, 279 250.0 / 621 - 18575.0 / 48384, 125.0 / 594 - 13525.0 / 55296, 280 0.0 - 277.0 / 14336, 512.0 / 1771 - 0.25 }; 281 282 /** The index of the error stored in the auxiliary variables. */ 283 private static final int _ERROR_INDEX = _TIME_INCREMENTS.length; 284 285 /** The order of the algorithm. */ 286 private static final int _ORDER = 5; 287 288 /** The round counter. */ 289 private int _roundCount = 0; 290}