numericalsolver.hpp

dysii是一款非常出色的滤波函数库

#ifndef INDII_ML_ODE_NUMERICALSOLVER_HPP
#define INDII_ML_ODE_NUMERICALSOLVER_HPP

#include "../aux/vector.hpp"

#include <gsl/gsl_odeiv.h>

namespace indii {
  namespace ml {
    namespace ode {

/**
 * Type of functions for calculating derivatives.
 */
typedef int df_t(double t, const double y[], double dydt[], void* params);

/**
 * Abstract numerical solver for a system of differential equations.
 *
 * @author Lawrence Murray <[e-mail]>
 * @version $Rev: 404 $
 * @date $Date: 2008-03-05 14:52:55 +0000 (Wed, 05 Mar 2008) $
 *
 * NumericalSolver should not be instantiated directly. Instead,
 * instantiate one of its derived classes.
 */
class NumericalSolver {
public:
  /**
   * Constructor.
   *
   * @param dimensions Dimensionality of the state space.
   *
   * The time is initialised to zero, but the state is uninitialised
   * and should be set with setVariable() or setState().
   */
  NumericalSolver(const unsigned int dimensions);

  /**
   * Constructor.
   *
   * @param y0 Initial state.
   *
   * The time is initialised to zero and the state to that given.
   */
  NumericalSolver(const indii::ml::aux::vector& y0);

  /**
   * Destructor.
   */
  virtual ~NumericalSolver();

  /**
   * Get the current time.
   */
  double getTime();

  /**
   * Set the current time.
   *
   * @param t Time.
   *
   * @see setVariable() or setState() to update the state for the new
   * time also.
   */
  void setTime(const double t);

  /**
   * Get the value of a state variable at the current time.
   *
   * @param index Index of the state variable to retrieve.
   *
   * @return The value of the state variable at the current time.
   */
  double getVariable(const unsigned int index);

  /**
   * Set the value of a state variable.
   *
   * @param index Index of the state variable to set.
   * @param value The value to which to set the state variable.
   */
  void setVariable(const unsigned int index, const double value);

  /**
   * Get the state at the current time.
   *
   * @return Current state of the system.
   */
  indii::ml::aux::vector getState();

  /**
   * Set the state at the current time.
   *
   * @param y New state of the system.
   */
  void setState(const indii::ml::aux::vector& y);

  /**
   * Advance the system one time step. The size of the step is chosen
   * to be optimal given the provided error bounds. State variables
   * are updated to this new time after completion of the step.
   *
   * @param upper Upper bound on time. This must be greater than the
   * current time. The current time is guaranteed not to exceed this
   * value after the step is complete. It may equal this time, and
   * indeed if step() is continuously called with the same upper bound
   * it will eventually do so.
   *
   * @return The new current time.
   */
  virtual double step(double upper);

  /**
   * Advance the system to a particular time.
   *
   * This convenience method internally calls step() as many times as
   * necessary to advance the system to the given time. This is useful
   * if fixed time steps are required rather than the variable steps
   * taken by step().
   *
   * @param to The time to which to advance the system. This must be greater
   * than the current time. The current time is guaranteed to match this value
   * at the end of the call.
   */
  void stepTo(double to);

  /**
   * Rewind the system one time step. This works in the same fashion
   * at step(), except that the step is backwards in time.
   *
   * @param lower Lower bound on time. This must be less than the
   * current time. The current time is guaranteed not to be below this
   * value after the step is complete. It may equal this time, and
   * indeed if stepBack() is continuously called with the same lower
   * bound it will eventually do so.
   *
   * @return The new current time.
   */
  virtual double stepBack(double lower);

  /**
   * Rewind the system to a specific time.
   *
   * This convenience method internally calls stepBack() as many times
   * as necessary to rewind the system to a given point in time. This
   * is useful if fixed time steps are required rather than the
   * variable steps taken by stepBack().
   *
   * @param to The time to which to rewind the system. This must be
   * less than the current time. The current time is guaranteed to
   * match this value at the end of the call.
   */
  void stepBackTo(double to);

  /**
   * Set the error bounds. Smaller error bounds produce more accurate
   * estimates, but reduce the size of time steps, so that more steps
   * are required in order to estimate functions over the same length
   * of time.
   *
   * @param maxAbsoluteError The maximum permitted absolute error of
   * estimated values compared to real values.
   * @param maxRelativeError The maximum permitted relative error of
   * estimated values compared to real values.
   */
  void setErrorBounds(double maxAbsoluteError = 1e-6,
		      double maxRelativeError = 1e-6);

  /**
   * Set the proposed size for the next time step. This is useful for
   * handling discontinuities, where step() may be called to advance
   * to the discontinuity, then setStepSize() used to propose a small
   * step size ensuring that the discontinuity is tiptoed across
   * rather than leapt.
   *
   * It should be called immediately before any call to step() or
   * stepTo(). Any other methods called in between may themselves adjust the
   * proposed step size.
   *
   * Note that this sets only the proposed step size -- it will be
   * optimised, within some constraints, relative to the permitted
   * error bounds.
   *
   * @param stepSize The proposed time step size for the next call to
   * step().
   */
  void setStepSize(double stepSize = 1.0e-7);

  /**
   * Convenience method to indicate a discontinuity at the current
   * time. Internally, this simply calls setStepSize() with the
   * default step size.
   */
  void setDiscontinuity();

protected:
  /**
   * Dimensionality of the system of differential equations.
   */
  const size_t dimensions;

  /**
   * \f$t\f$; the current time
   */
  double t;

  /**
   * Array of state variables of the system.
   */
  double* y;

  /**
   * Base time used for backward steps.
   */
  double base;

  /**
   * Proposed size for the next time step.
   */
  double stepSize;

  /**
   * GSL ordinary differential equations step structure.
   */
  gsl_odeiv_step* gslStep;

  /**
   * GSL ordinary differential equations control structure.
   */
  gsl_odeiv_control* gslControl;

  /**
   * GSL ordinary differential equations evolution structure.
   */
  gsl_odeiv_evolve* gslEvolve;
 
  /**
   * GSL system of ordinary differential equations for forward steps
   * in time.
   */
  gsl_odeiv_system gslForwardSystem;

  /**
   * GSL system of ordinary differential equations for backward steps
   * in time.
   */
  gsl_odeiv_system gslBackwardSystem;

  /**
   * Set function for calculating derivatives when moving forward in
   * time.
   *
   * @param f The function.
   *
   * Required by gsl_odeiv_system data type of GSL, which is in turn
   * passed to the gsl_odeiv_evolve_apply() function. Acts as a
   * wrapper to call the calculateDerivatives() function of a particular
   * DifferentialModel object.
   *
   * @see gsl_odeiv_system data type of GSL.
   * @see gsl_odeiv_evolve_apply() of GSL.
   */
  void setForwardFunction(df_t f);

  /**
   * Set function for calculating derivatives when moving forward in
   * time.
   *
   * @param f The function.
   *
   * @see setForwardFunction()
   */
  void setBackwardFunction(df_t f);

  /**
   * Set the stepping method.
   *
   * @param stepType The stepping method.
   */
  void setStepType(const gsl_odeiv_step_type* stepType);

  /**
   * Initialise the system. Should be called before any calls to
   * step() or stepTo().
   */
  void init();

  /**
   * Terminate the system. Should be called after work with the system
   * has been finished. Will be called by the destructor if it hasn't
   * been already.
   */
  void terminate();

  /**
   * Reset the model. This is used internally to reset GSL structures when the
   * next step will not be a continuation of the last, usually when the time
   * or state is changed with a call to setTime(), setVariable() or
   * setState().
   */
  virtual void reset();

};


    }
  }
}

#endif