Добрый день)
Понадобился мне тут на днях класс, для решения систем обыкновенных диф. уравнений (ОДЕ) методом, желательно наиболее устойчивым)) Погуглив на эту тему, было решено остановиться на классическом методе Рунге-Кутта 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;
}
Вопрос: Как улучшить?
Заранее благодарю)
