switchstate.java

Apache的common math数学软件包

/*
 * Licensed to the Apache Software Foundation (ASF) under one or more
 * contributor license agreements.  See the NOTICE file distributed with
 * this work for additional information regarding copyright ownership.
 * The ASF licenses this file to You under the Apache License, Version 2.0
 * (the "License"); you may not use this file except in compliance with
 * the License.  You may obtain a copy of the License at
 *
 *      http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

package org.apache.commons.math.ode;

import java.io.Serializable;

import org.apache.commons.math.ConvergenceException;
import org.apache.commons.math.FunctionEvaluationException;
import org.apache.commons.math.analysis.BrentSolver;
import org.apache.commons.math.analysis.UnivariateRealFunction;
import org.apache.commons.math.analysis.UnivariateRealSolver;

/** This class handles the state for one {@link SwitchingFunction
 * switching function} during integration steps.
 *
 * <p>Each time the integrator proposes a step, the switching function
 * should be checked. This class handles the state of one function
 * during one integration step, with references to the state at the
 * end of the preceding step. This information is used to determine if
 * the function should trigger an event or not during the proposed
 * step (and hence the step should be reduced to ensure the event
 * occurs at a bound rather than inside the step).</p>
 *
 * @version $Revision: 620312 $ $Date: 2008-02-10 12:28:59 -0700 (Sun, 10 Feb 2008) $
 * @since 1.2
 */
class SwitchState implements Serializable {

  /** Serializable version identifier. */
  private static final long serialVersionUID = -7307007422156119622L;

  /** Switching function. */
  private SwitchingFunction function;

  /** Maximal time interval between switching function checks. */
  private double maxCheckInterval;

  /** Convergence threshold for event localisation. */
  private double convergence;

  /** Upper limit in the iteration count for event localisation. */
  private int maxIterationCount;

  /** Time at the beginning of the step. */
  private double t0;

  /** Value of the switching function at the beginning of the step. */
  private double g0;

  /** Simulated sign of g0 (we cheat when crossing events). */
  private boolean g0Positive;

  /** Indicator of event expected during the step. */
  private boolean pendingEvent;

  /** Occurrence time of the pending event. */
  private double pendingEventTime;

  /** Occurrence time of the previous event. */
  private double previousEventTime;

  /** Variation direction around pending event.
   *  (this is considered with respect to the integration direction)
   */
  private boolean increasing;

  /** Next action indicator. */
  private int nextAction;

  /** Simple constructor.
   * @param function switching function
   * @param maxCheckInterval maximal time interval between switching
   * function checks (this interval prevents missing sign changes in
   * case the integration steps becomes very large)
   * @param convergence convergence threshold in the event time search
   * @param maxIterationCount upper limit of the iteration count in
   * the event time search
   */
  public SwitchState(SwitchingFunction function, double maxCheckInterval,
                     double convergence, int maxIterationCount) {
    this.function          = function;
    this.maxCheckInterval  = maxCheckInterval;
    this.convergence       = Math.abs(convergence);
    this.maxIterationCount = maxIterationCount;

    // some dummy values ...
    t0                = Double.NaN;
    g0                = Double.NaN;
    g0Positive        = true;
    pendingEvent      = false;
    pendingEventTime  = Double.NaN;
    previousEventTime = Double.NaN;
    increasing        = true;
    nextAction        = SwitchingFunction.CONTINUE;

  }

  /** Reinitialize the beginning of the step.
   * @param t0 value of the independent <i>time</i> variable at the
   * beginning of the step
   * @param y0 array containing the current value of the state vector
   * at the beginning of the step
   * @exception FunctionEvaluationException if the switching function
   * value cannot be evaluated at the beginning of the step
   */
  public void reinitializeBegin(double t0, double[] y0)
    throws FunctionEvaluationException {
    this.t0 = t0;
    g0 = function.g(t0, y0);
    g0Positive = (g0 >= 0);
  }

  /** Evaluate the impact of the proposed step on the switching function.
   * @param interpolator step interpolator for the proposed step
   * @return true if the switching function triggers an event before
   * the end of the proposed step (this implies the step should be
   * rejected)
   * @exception DerivativeException if the interpolator fails to
   * compute the function somewhere within the step
   * @exception FunctionEvaluationException if the switching function
   * cannot be evaluated
   * @exception ConvergenceException if an event cannot be located
   */
  public boolean evaluateStep(final StepInterpolator interpolator)
    throws DerivativeException, FunctionEvaluationException, ConvergenceException {

    try {

      double t1 = interpolator.getCurrentTime();
      int    n  = Math.max(1, (int) Math.ceil(Math.abs(t1 - t0) / maxCheckInterval));
      double h  = (t1 - t0) / n;

      double ta = t0;
      double ga = g0;
      double tb = t0 + ((t1 > t0) ? convergence : -convergence);
      for (int i = 0; i < n; ++i) {

        // evaluate function value at the end of the substep
        tb += h;
        interpolator.setInterpolatedTime(tb);
        double gb = function.g(tb, interpolator.getInterpolatedState());

        // check events occurrence
        if (g0Positive ^ (gb >= 0)) {
          // there is a sign change: an event is expected during this step

          // variation direction, with respect to the integration direction
          increasing = (gb >= ga);

          UnivariateRealSolver solver = new BrentSolver(new UnivariateRealFunction() {
              public double value(double t) throws FunctionEvaluationException {
                  try {
                      interpolator.setInterpolatedTime(t);
                      return function.g(t, interpolator.getInterpolatedState());
                  } catch (DerivativeException e) {
                      throw new FunctionEvaluationException(t, e);
                  }
              }
          });
          solver.setAbsoluteAccuracy(convergence);
          solver.setMaximalIterationCount(maxIterationCount);
          double root = solver.solve(ta, tb);
          if (Double.isNaN(previousEventTime) || (Math.abs(previousEventTime - root) > convergence)) {
              pendingEventTime = root;
              if (pendingEvent && (Math.abs(t1 - pendingEventTime) <= convergence)) {
                  // we were already waiting for this event which was
                  // found during a previous call for a step that was
                  // rejected, this step must now be accepted since it
                  // properly ends exactly at the event occurrence
                  return false;
              }
              // either we were not waiting for the event or it has
              // moved in such a way the step cannot be accepted
              pendingEvent = true;
              return true;
          }

        } else {
          // no sign change: there is no event for now
          ta = tb;
          ga = gb;
        }

      }

      // no event during the whole step
      pendingEvent     = false;
      pendingEventTime = Double.NaN;
      return false;

    } catch (FunctionEvaluationException e) {
     Throwable cause = e.getCause();
      if ((cause != null) && (cause instanceof DerivativeException)) {
        throw (DerivativeException) cause;
      }
      throw e;
    }

  }

  /** Get the occurrence time of the event triggered in the current
   * step.
   * @return occurrence time of the event triggered in the current
   * step.
   */
  public double getEventTime() {
    return pendingEventTime;
  }

  /** Acknowledge the fact the step has been accepted by the integrator.
   * @param t value of the independent <i>time</i> variable at the
   * end of the step
   * @param y array containing the current value of the state vector
   * at the end of the step
   * @exception FunctionEvaluationException if the value of the switching
   * function cannot be evaluated
   */
  public void stepAccepted(double t, double[] y)
    throws FunctionEvaluationException {

    t0 = t;
    g0 = function.g(t, y);

    if (pendingEvent) {
      // force the sign to its value "just after the event"
      previousEventTime = t;
      g0Positive        = increasing;
      nextAction        = function.eventOccurred(t, y);
    } else {
      g0Positive = (g0 >= 0);
      nextAction = SwitchingFunction.CONTINUE;
    }
  }

  /** Check if the integration should be stopped at the end of the
   * current step.
   * @return true if the integration should be stopped
   */
  public boolean stop() {
    return nextAction == SwitchingFunction.STOP;
  }

  /** Let the switching function reset the state if it wants.
   * @param t value of the independent <i>time</i> variable at the
   * beginning of the next step
   * @param y array were to put the desired state vector at the beginning
   * of the next step
   * @return true if the integrator should reset the derivatives too
   */
  public boolean reset(double t, double[] y) {

    if (! pendingEvent) {
      return false;
    }

    if (nextAction == SwitchingFunction.RESET_STATE) {
      function.resetState(t, y);
    }
    pendingEvent      = false;
    pendingEventTime  = Double.NaN;

    return (nextAction == SwitchingFunction.RESET_STATE) ||
           (nextAction == SwitchingFunction.RESET_DERIVATIVES);

  }

}