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; 19 20 import java.io.Serializable; 21 22 import org.apache.commons.math.ConvergenceException; 23 import org.apache.commons.math.FunctionEvaluationException; 24 import org.apache.commons.math.analysis.BrentSolver; 25 import org.apache.commons.math.analysis.UnivariateRealFunction; 26 import org.apache.commons.math.analysis.UnivariateRealSolver; 27 28 /** This class handles the state for one {@link SwitchingFunction 29 * switching function} during integration steps. 30 * 31 * <p>Each time the integrator proposes a step, the switching function 32 * should be checked. This class handles the state of one function 33 * during one integration step, with references to the state at the 34 * end of the preceding step. This information is used to determine if 35 * the function should trigger an event or not during the proposed 36 * step (and hence the step should be reduced to ensure the event 37 * occurs at a bound rather than inside the step).</p> 38 * 39 * @version $Revision: 620312 $ $Date: 2008-02-10 12:28:59 -0700 (Sun, 10 Feb 2008) $ 40 * @since 1.2 41 */ 42 class SwitchState implements Serializable { 43 44 /** Serializable version identifier. */ 45 private static final long serialVersionUID = -7307007422156119622L; 46 47 /** Switching function. */ 48 private SwitchingFunction function; 49 50 /** Maximal time interval between switching function checks. */ 51 private double maxCheckInterval; 52 53 /** Convergence threshold for event localisation. */ 54 private double convergence; 55 56 /** Upper limit in the iteration count for event localisation. */ 57 private int maxIterationCount; 58 59 /** Time at the beginning of the step. */ 60 private double t0; 61 62 /** Value of the switching function at the beginning of the step. */ 63 private double g0; 64 65 /** Simulated sign of g0 (we cheat when crossing events). */ 66 private boolean g0Positive; 67 68 /** Indicator of event expected during the step. */ 69 private boolean pendingEvent; 70 71 /** Occurrence time of the pending event. */ 72 private double pendingEventTime; 73 74 /** Occurrence time of the previous event. */ 75 private double previousEventTime; 76 77 /** Variation direction around pending event. 78 * (this is considered with respect to the integration direction) 79 */ 80 private boolean increasing; 81 82 /** Next action indicator. */ 83 private int nextAction; 84 85 /** Simple constructor. 86 * @param function switching function 87 * @param maxCheckInterval maximal time interval between switching 88 * function checks (this interval prevents missing sign changes in 89 * case the integration steps becomes very large) 90 * @param convergence convergence threshold in the event time search 91 * @param maxIterationCount upper limit of the iteration count in 92 * the event time search 93 */ 94 public SwitchState(SwitchingFunction function, double maxCheckInterval, 95 double convergence, int maxIterationCount) { 96 this.function = function; 97 this.maxCheckInterval = maxCheckInterval; 98 this.convergence = Math.abs(convergence); 99 this.maxIterationCount = maxIterationCount; 100 101 // some dummy values ... 102 t0 = Double.NaN; 103 g0 = Double.NaN; 104 g0Positive = true; 105 pendingEvent = false; 106 pendingEventTime = Double.NaN; 107 previousEventTime = Double.NaN; 108 increasing = true; 109 nextAction = SwitchingFunction.CONTINUE; 110 111 } 112 113 /** Reinitialize the beginning of the step. 114 * @param t0 value of the independent <i>time</i> variable at the 115 * beginning of the step 116 * @param y0 array containing the current value of the state vector 117 * at the beginning of the step 118 * @exception FunctionEvaluationException if the switching function 119 * value cannot be evaluated at the beginning of the step 120 */ 121 public void reinitializeBegin(double t0, double[] y0) 122 throws FunctionEvaluationException { 123 this.t0 = t0; 124 g0 = function.g(t0, y0); 125 g0Positive = (g0 >= 0); 126 } 127 128 /** Evaluate the impact of the proposed step on the switching function. 129 * @param interpolator step interpolator for the proposed step 130 * @return true if the switching function triggers an event before 131 * the end of the proposed step (this implies the step should be 132 * rejected) 133 * @exception DerivativeException if the interpolator fails to 134 * compute the function somewhere within the step 135 * @exception FunctionEvaluationException if the switching function 136 * cannot be evaluated 137 * @exception ConvergenceException if an event cannot be located 138 */ 139 public boolean evaluateStep(final StepInterpolator interpolator) 140 throws DerivativeException, FunctionEvaluationException, ConvergenceException { 141 142 try { 143 144 double t1 = interpolator.getCurrentTime(); 145 int n = Math.max(1, (int) Math.ceil(Math.abs(t1 - t0) / maxCheckInterval)); 146 double h = (t1 - t0) / n; 147 148 double ta = t0; 149 double ga = g0; 150 double tb = t0 + ((t1 > t0) ? convergence : -convergence); 151 for (int i = 0; i < n; ++i) { 152 153 // evaluate function value at the end of the substep 154 tb += h; 155 interpolator.setInterpolatedTime(tb); 156 double gb = function.g(tb, interpolator.getInterpolatedState()); 157 158 // check events occurrence 159 if (g0Positive ^ (gb >= 0)) { 160 // there is a sign change: an event is expected during this step 161 162 // variation direction, with respect to the integration direction 163 increasing = (gb >= ga); 164 165 UnivariateRealSolver solver = new BrentSolver(new UnivariateRealFunction() { 166 public double value(double t) throws FunctionEvaluationException { 167 try { 168 interpolator.setInterpolatedTime(t); 169 return function.g(t, interpolator.getInterpolatedState()); 170 } catch (DerivativeException e) { 171 throw new FunctionEvaluationException(t, e); 172 } 173 } 174 }); 175 solver.setAbsoluteAccuracy(convergence); 176 solver.setMaximalIterationCount(maxIterationCount); 177 double root = solver.solve(ta, tb); 178 if (Double.isNaN(previousEventTime) || (Math.abs(previousEventTime - root) > convergence)) { 179 pendingEventTime = root; 180 if (pendingEvent && (Math.abs(t1 - pendingEventTime) <= convergence)) { 181 // we were already waiting for this event which was 182 // found during a previous call for a step that was 183 // rejected, this step must now be accepted since it 184 // properly ends exactly at the event occurrence 185 return false; 186 } 187 // either we were not waiting for the event or it has 188 // moved in such a way the step cannot be accepted 189 pendingEvent = true; 190 return true; 191 } 192 193 } else { 194 // no sign change: there is no event for now 195 ta = tb; 196 ga = gb; 197 } 198 199 } 200 201 // no event during the whole step 202 pendingEvent = false; 203 pendingEventTime = Double.NaN; 204 return false; 205 206 } catch (FunctionEvaluationException e) { 207 Throwable cause = e.getCause(); 208 if ((cause != null) && (cause instanceof DerivativeException)) { 209 throw (DerivativeException) cause; 210 } 211 throw e; 212 } 213 214 } 215 216 /** Get the occurrence time of the event triggered in the current 217 * step. 218 * @return occurrence time of the event triggered in the current 219 * step. 220 */ 221 public double getEventTime() { 222 return pendingEventTime; 223 } 224 225 /** Acknowledge the fact the step has been accepted by the integrator. 226 * @param t value of the independent <i>time</i> variable at the 227 * end of the step 228 * @param y array containing the current value of the state vector 229 * at the end of the step 230 * @exception FunctionEvaluationException if the value of the switching 231 * function cannot be evaluated 232 */ 233 public void stepAccepted(double t, double[] y) 234 throws FunctionEvaluationException { 235 236 t0 = t; 237 g0 = function.g(t, y); 238 239 if (pendingEvent) { 240 // force the sign to its value "just after the event" 241 previousEventTime = t; 242 g0Positive = increasing; 243 nextAction = function.eventOccurred(t, y); 244 } else { 245 g0Positive = (g0 >= 0); 246 nextAction = SwitchingFunction.CONTINUE; 247 } 248 } 249 250 /** Check if the integration should be stopped at the end of the 251 * current step. 252 * @return true if the integration should be stopped 253 */ 254 public boolean stop() { 255 return nextAction == SwitchingFunction.STOP; 256 } 257 258 /** Let the switching function reset the state if it wants. 259 * @param t value of the independent <i>time</i> variable at the 260 * beginning of the next step 261 * @param y array were to put the desired state vector at the beginning 262 * of the next step 263 * @return true if the integrator should reset the derivatives too 264 */ 265 public boolean reset(double t, double[] y) { 266 267 if (! pendingEvent) { 268 return false; 269 } 270 271 if (nextAction == SwitchingFunction.RESET_STATE) { 272 function.resetState(t, y); 273 } 274 pendingEvent = false; 275 pendingEventTime = Double.NaN; 276 277 return (nextAction == SwitchingFunction.RESET_STATE) || 278 (nextAction == SwitchingFunction.RESET_DERIVATIVES); 279 280 } 281 282 }