001/* Explicit variable step size Runge-Kutta 2(3) ODE solver.
002
003 Copyright (c) 1998-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//// ExplicitRK23Solver
037
038/**
039 This class implements the Explicit Runge-Kutta 2(3) ODE solving method.
040 For an ODE of the form:
041 <pre>
042 dx/dt = f(x, t), x(0) = x0
043 </pre>
044 it does the following:
045 <pre>
046 K0 = f(x(n), tn);
047 K1 = f(x(n)+0.5*h*K0, tn+0.5*h);
048 K2 = f(x(n)+0.75*h*K1, tn+0.75*h);
049 x(n+1) = x(n)+(2/9)*h*K0+(1/3)*h*K0+(4/9)*h*K2;
050 K3 = f(x(n+1), tn+h);
051 </pre>,
052 and error control:
053 <pre>
054 LTE = h*[(-5.0/72.0)*K0 + (1.0/12.0)*K1 + (1.0/9.0)*K2 + (-1.0/8.0)*K3]
055 </pre>
056 <P>
057 If the LTE is less than the error tolerance, then this step is considered
058 successful, and the next integration step is predicted as:
059 <pre>
060 h' = 0.8*Math.pow((ErrorTolerance/LTE), 1.0/3.0)
061 </pre>
062 This is a second order method, but uses a third order procedure to estimate
063 the local truncation error.
064
065 @author  Jie Liu, Haiyang Zheng, Edward A. Lee
066 @version $Id$
067 @since Ptolemy II 6.0
068 @Pt.ProposedRating Green (hyzheng)
069 @Pt.AcceptedRating Green (hyzheng)
070 */
071public class ExplicitRK23Solver extends ContinuousODESolver {
072
073    ///////////////////////////////////////////////////////////////////
074    ////                         public methods                    ////
075
076    /** Return the number of time increments plus one (to store the
077     *  truncation error).
078     *  @return The number of time increments plus one.
079     */
080    @Override
081    public final int getIntegratorAuxVariableCount() {
082        // Allow one for the truncation error
083        return _TIME_INCREMENTS.length + 1;
084    }
085
086    /** Fire the given integrator. This method performs the ODE solving
087     *  algorithm described in the class comment.
088     *  @param integrator The integrator of that calls this method.
089     *  @exception IllegalActionException If there is no director, or can not
090     *  read input, or can not send output.
091     */
092    @Override
093    public void integratorIntegrate(ContinuousIntegrator integrator)
094            throws IllegalActionException {
095        double outputValue;
096        double xn = integrator.getState();
097        double h = _director.getCurrentStepSize();
098        double[] k = integrator.getAuxVariables();
099        integrator.setAuxVariables(_roundCount, integrator.getDerivative());
100
101        switch (_roundCount) {
102        case 0:
103            outputValue = xn + h * k[0] * _B[0][0];
104            break;
105
106        case 1:
107            outputValue = xn + h * (k[0] * _B[1][0] + k[1] * _B[1][1]);
108            break;
109
110        case 2:
111            outputValue = xn
112                    + h * (k[0] * _B[2][0] + k[1] * _B[2][1] + k[2] * _B[2][2]);
113            break;
114
115        case 3:
116            outputValue = integrator.getTentativeState();
117            return;
118
119        default:
120            throw new InvalidStateException("Execution sequence out of range.");
121        }
122
123        integrator.setTentativeState(outputValue);
124    }
125
126    /** Return true if the integration is accurate for the given
127     *  integrator. This estimates the local truncation error for that
128     *  integrator and compare it with the error tolerance.
129     *
130     *  @param integrator The integrator of that calls this method.
131     *  @return True if the integration is successful.
132     */
133    @Override
134    public boolean integratorIsAccurate(ContinuousIntegrator integrator) {
135        double tolerance = _director.getErrorTolerance();
136        double h = _director.getCurrentStepSize();
137        double[] k = integrator.getAuxVariables();
138        double error = h * Math
139                .abs(k[0] * _E[0] + k[1] * _E[1] + k[2] * _E[2] + k[3] * _E[3]);
140
141        integrator.setAuxVariables(_ERROR_INDEX, error);
142        if (_isDebugging()) {
143            _debug("Integrator: " + integrator.getName()
144                    + " local truncation error = " + error);
145        }
146
147        if (error < tolerance) {
148            if (_isDebugging()) {
149                _debug("Integrator: " + integrator.getName()
150                        + " report a success.");
151            }
152            return true;
153        } else {
154            if (_isDebugging()) {
155                _debug("Integrator: " + integrator.getName()
156                        + " reports a failure.");
157            }
158            return false;
159        }
160    }
161
162    /** Provide the predictedStepSize() method for the integrators
163     *  under this solver. It uses the algorithm in the class comments
164     *  to predict the next step size based on the current estimation
165     *  of the local truncation error.
166     *
167     *  @param integrator The integrator of that calls this method.
168     *  @return The next step size suggested by the given integrator.
169     */
170    @Override
171    public double integratorSuggestedStepSize(ContinuousIntegrator integrator) {
172        double error = integrator.getAuxVariables()[_ERROR_INDEX];
173        double h = _director.getCurrentStepSize();
174        double tolerance = _director.getErrorTolerance();
175        double newh = 5.0 * h;
176
177        if (error > tolerance) {
178            newh = 0.8 * Math.pow(tolerance / error, 1.0 / _ORDER);
179            if (newh > h) {
180                newh = 0.5 * h;
181            }
182        }
183
184        if (_isDebugging()) {
185            _debug("integrator: " + integrator.getName()
186                    + " suggests next step size = " + newh);
187        }
188        return newh;
189    }
190
191    ///////////////////////////////////////////////////////////////////
192    ////                         protected methods                 ////
193
194    /** Return the current round.
195     *  @return The current round.
196     */
197    @Override
198    protected int _getRound() {
199        return _roundCount;
200    }
201
202    /** Get the current round factor.
203     */
204    @Override
205    protected final double _getRoundTimeIncrement() {
206        return _TIME_INCREMENTS[_roundCount];
207    }
208
209    /** Return true if the current integration step is finished.
210     *  This method will return true if _incrementRound() has been
211     *  called 4 or more times since _reset().
212     *  @see #_reset()
213     */
214    @Override
215    protected final boolean _isStepFinished() {
216        return _roundCount >= _TIME_INCREMENTS.length;
217    }
218
219    /** Reset the solver, indicating to it that we are starting an
220     *  integration step. This method resets the round counter.
221     */
222    @Override
223    protected final void _reset() {
224        _roundCount = 0;
225    }
226
227    /** Set the round for the next integration step.
228     *  @param round The round for the next integration step.
229     */
230    @Override
231    protected void _setRound(int round) {
232        _roundCount = round;
233    }
234
235    ///////////////////////////////////////////////////////////////////
236    ////                         protected variables               ////
237
238    /** The ratio of time increments within one integration step. */
239    protected static final double[] _TIME_INCREMENTS = { 0.5, 0.75, 1.0, 1.0 };
240
241    ///////////////////////////////////////////////////////////////////
242    ////                         private variables                 ////
243
244    /** B coefficients. */
245    private static final double[][] _B = { { 0.5 }, { 0, 0.75 },
246            { 2.0 / 9.0, 1.0 / 3.0, 4.0 / 9.0 } };
247
248    /** E coefficients. */
249    private static final double[] _E = { -5.0 / 72.0, 1.0 / 12.0, 1.0 / 9.0,
250            -1.0 / 8.0 };
251
252    /** The index of the error stored in the auxiliary variables. */
253    private static final int _ERROR_INDEX = _TIME_INCREMENTS.length;
254
255    /** The order of the algorithm. */
256    private static final int _ORDER = 3;
257
258    /** The round counter. */
259    private int _roundCount = 0;
260}