summaryrefslogtreecommitdiffstats
path: root/kscreensaver/kdesavers/rkodesolver.h
diff options
context:
space:
mode:
Diffstat (limited to 'kscreensaver/kdesavers/rkodesolver.h')
-rw-r--r--kscreensaver/kdesavers/rkodesolver.h187
1 files changed, 187 insertions, 0 deletions
diff --git a/kscreensaver/kdesavers/rkodesolver.h b/kscreensaver/kdesavers/rkodesolver.h
new file mode 100644
index 00000000..d1f9d556
--- /dev/null
+++ b/kscreensaver/kdesavers/rkodesolver.h
@@ -0,0 +1,187 @@
+//============================================================================
+//
+// Ordinary differential equation solver using the Runge-Kutta method.
+// $Id$
+// Copyright (C) 2004 Georg Drenkhahn
+//
+// This file is free software; you can redistribute it and/or modify it under
+// the terms of the GNU General Public License version 2 as published by the
+// Free Software Foundation.
+//
+//============================================================================
+
+#ifndef RKODESOLVER_H
+#define RKODESOLVER_H
+
+// STL headers
+#include <valarray>
+
+/** @brief Solver class to integrate a first-order ordinary differential
+ * equation (ODE) by means of a 6. order Runge-Kutta method.
+ *
+ * The ODE system must be given as the derivative
+ * dy/dx = f(x,y)
+ * with x in R and y in R^n.
+ *
+ * Within this class the function f() is a pure virtual function, which must be
+ * reimplemented in a derived class.
+ *
+ * No other special data type for vectors or matrices are needed besides the STL
+ * class std::valarray. */
+template<typename T>
+class RkOdeSolver
+{
+ public:
+ /** @brief Constructor
+ * @param x Initial integration parameter
+ * @param y Initial function values of function to integrate
+ * @param dx Initial guess for step size. Will be automatically adjusted to
+ * guarantee required precision.
+ * @param eps Relative precision
+ *
+ * Initialises the solver with start conditions. */
+ RkOdeSolver(const T& x=0.0,
+ const std::valarray<T>& y=std::valarray<T>(0),
+ const T& dx=0,
+ const T& eps=1e-6);
+ /** @brief Destructor */
+ virtual ~RkOdeSolver(void);
+
+ /** @brief Integrates the ordinary differential equation from the current x
+ * value to x+@a dx.
+ * @param dx x-interval size to integrate over starting from x. dx may be
+ * negative.
+ *
+ * The integration is performed by calling rkStepCheck() repeatedly until the
+ * desired x value is reached. */
+ void integrate(const T& dx);
+
+ // Accessors
+
+ // get/set x value
+ /** @brief Get current x value.
+ * @return Reference of x value. */
+ const T& X(void) const;
+ /** @brief Set current x value.
+ * @param a The value to be set. */
+ void X(const T& a);
+
+ // get/set y value
+ /** @brief Get current y value.
+ * @return Reference of y vector. */
+ const std::valarray<T>& Y(void) const;
+ /** @brief Set current y values.
+ * @param a The vector to be set. */
+ void Y(const std::valarray<T>& a);
+
+ /** @brief Get current dy/dx value.
+ * @return Reference of dy/dx vector. */
+ const std::valarray<T>& dYdX(void) const;
+
+ // get/set dx value
+ /** @brief Get current estimated step size dX.
+ * @return Reference of dX value. */
+ const T& dX(void) const;
+ /** @brief Set estimated step size dX.
+ * @param a The value to be set. */
+ void dX(const T& a);
+
+ // get/set eps value
+ /** @brief Get current presision.
+ * @return Reference of precision value. */
+ const T& Eps(void) const;
+ /** @brief Set estimated presision.
+ * @param a The value to be set. */
+ void Eps(const T& a);
+
+ protected:
+ // purely virtual function which is integrated
+ /** @brief ODE function
+ * @param x Integration value
+ * @param y Function value
+ * @return Derivation
+ *
+ * This purely virtual function returns the value of dy/dx for the given
+ * parameter values of x and y. */
+ virtual std::valarray<T>
+ f(const T& x, const std::valarray<T>& y) const = 0;
+
+ private:
+ /** @brief Perform one integration step with a tolerable relative error given
+ * by ::mErr.
+ * @param dx Maximal step size, may be positive or negative depending on
+ * integration direction.
+ * @return Flag indicating if made absolute integration step was equal |@a dx
+ * | (true) less than |@a dx | (false).
+ *
+ * A new estimate for the maximum next step size is saved to ::mStep. The
+ * new values for x, y and f are saved in ::mX, ::mY and ::mDydx. */
+ bool rkStepCheck(const T& dx);
+ /** @brief Perform one Runge-Kutta integration step forward with step size
+ * ::mStep
+ * @param dx Step size relative to current x value.
+ * @param yerr Reference to vector in which the estimated error made in y is
+ * returned.
+ * @return The y value after the step at x+@a dx.
+ *
+ * Stored current x,y values are not adjusted. */
+ std::valarray<T> rkStep(const T& dx, std::valarray<T>& yerr) const;
+
+ /** current x value */
+ T mX;
+ /** current y value */
+ std::valarray<T> mY;
+ /** current value of dy/dx */
+ std::valarray<T> mDydx;
+
+ /** allowed relative error */
+ T mEps;
+ /** estimated step size for next Runge-Kutta step */
+ T mStep;
+};
+
+// inline accessors
+
+template<typename T>
+inline const T&
+RkOdeSolver<T>::X(void) const
+{
+ return mX;
+}
+
+template<typename T>
+inline void
+RkOdeSolver<T>::X(const T &a)
+{
+ mX = a;
+}
+
+template<typename T>
+inline const std::valarray<T>&
+RkOdeSolver<T>::Y(void) const
+{
+ return mY;
+}
+
+template<typename T>
+inline const std::valarray<T>&
+RkOdeSolver<T>::dYdX(void) const
+{
+ return mDydx;
+}
+
+template<typename T>
+inline const T&
+RkOdeSolver<T>::dX(void) const
+{
+ return mStep;
+}
+
+template<typename T>
+inline const T&
+RkOdeSolver<T>::Eps(void) const
+{
+ return mEps;
+}
+
+#endif