Russian Qt Forum

Программирование => Алгоритмы => Тема начата: m_ax от Апрель 24, 2010, 14:22



Название: Классичкский метод Рунге-Кутта 4-го порядка (RK4)
Отправлено: m_ax от Апрель 24, 2010, 14:22
Добрый день)

Понадобился мне тут на днях класс, для решения систем обыкновенных диф. уравнений (ОДЕ) методом, желательно наиболее устойчивым)) Погуглив на эту тему, было решено остановиться на классическом методе Рунге-Кутта 4-го порядка.
В общем выкладываю реализацию, которая первая пришла в голову.. Но мне что  то она не оч нравиться((
Как возможно сделать наиболее правильно, чтоб с одной стороны этим классом было удобно и интуитивно понятно пользоваться, а с другой чтоб производительность не страдала, а напротив?)

abstractode.h
Код
C++ (Qt)
#ifndef ABSTRACTODE_H
#define ABSTRACTODE_H
 
#include <vector>
 
class AbstractODE
{
public:
   AbstractODE();
   virtual ~AbstractODE();
   std::vector<double> operator()(const std::vector<double> &initPoint,
                                  const double &t0, const double &t, const double &h=0.001) const;
   virtual int dimension() const = 0;
 
protected:
   virtual std::vector<double> function(const std::vector<double> &, const double &) const = 0;
 
private:
};
 
#endif // ABSTRACTODE_H
 
   

abstractode.cpp
Код
C++ (Qt)
#include "abstractode.h"
#include <cmath>
 
AbstractODE::AbstractODE()
{
}
//-----------------------------------------------------------------------------
 
AbstractODE::~AbstractODE()
{
}
//-----------------------------------------------------------------------------
 
std::vector<double> AbstractODE::operator ()(const std::vector<double> &initPoint,
                                            const double &t0, const double &t, const double &h) const
{
   std::vector<double> x = initPoint;
   if ((int)x.size() != dimension())
       x.resize(dimension());
 
   std::vector<double> k1(dimension());
   std::vector<double> k2(dimension());
   std::vector<double> k3(dimension());
   std::vector<double> k4(dimension());
 
   std::vector<double> cx(dimension());
 
   double _t = t0;
   double _h = fabs(h);
   if (t < t0)
       _h *= -1;
 
   bool isStop = false;
 
   do {
 
       if (fabs(t - _t) <= fabs(h)) {
           _h = t - _t;
           isStop = true;
       }
 
       for (int i = 0; i < dimension(); i++) {
           k1[i] = function(x, _t)[i];
       }
       for (int i = 0; i < dimension(); i++) {
           cx[i] = x[i] + 0.5 * _h * k1[i];
       }
       double tau = _t + 0.5 * _h;
       for (int i = 0; i < dimension(); i++) {
           k2[i] = function(cx, tau)[i];
       }
       for (int i = 0; i < dimension(); i++) {
           cx[i] = x[i] + 0.5 * _h * k2[i];
       }
       for (int i = 0; i < dimension(); i++) {
           k3[i] = function(cx, tau)[i];
       }
       for (int i = 0; i < dimension(); i++) {
           cx[i] = x[i] + _h * k3[i];
       }
       tau = _t + _h;
       for (int i = 0; i < dimension(); i++) {
           k4[i] = function(cx, tau)[i];
       }
 
       for (int i = 0; i < dimension(); i++) {
           x[i] = x[i] + _h * (k1[i] + 2.0 * (k2[i] + k3[i]) + k4[i])/6.0;
       }
 
       _t += _h;
 
   } while (!isStop);
 
   return x;
}
//-----------------------------------------------------------------------------
 

Используется енто так:
Вначале наследуемся от AbstractODE и переопределяем две виртуальные функции:
1) int dimension() const;
2) std::vector<double> function(const std::vector<double> &x, const double &) const;
пример:
Код
C++ (Qt)
#include <iostream>
#include <fstream>
#include "abstractode.h"
#include <cmath>
 
using namespace std;
 
class MyODE : public AbstractODE
{
public:
   int dimension() const { return 2; }
protected:
   std::vector<double> function(const std::vector<double> &x, const double &) const;
 
};
 
inline std::vector<double> MyODE::function(const std::vector<double> &x,
                                          const double &/*t*/) const
{
   static std::vector<double> res(dimension());
 
   res[0] = x[1];
   res[1] = -x[0];
 
   return res;
}
 
 
int main(int argc, char *argv[])
{
   MyODE myODE;
   std::vector<double> initPoint(myODE.dimension());
   initPoint[0] = 0.5;
   initPoint[1] = 0.0;
 
   const double t0 = 0.0;
   const double t = 10.0;
 
   std::vector<double> x = myODE(initPoint, t0, t, 1e-2);
 
   cout << " x[0] = " << x[0] << endl;
   cout << " x[1] = " << x[1] << endl;
 
   return 0;
}
 

Вопрос: Как улучшить?

Заранее благодарю)