Солнечные часы ODE_vector преобразование в броненосец для приближения ODE
Я пытаюсь использовать библиотеки решателей Sundials ODE, чтобы аппроксимировать динамику расширения к уравнениям конкуренции Лотки-Вольтерра с множеством видов. Например. в случае двух видов
dx1 / dt = r1 * x1 * (1 — (x1 + a12 * x2))
dx2 / dt = r2 * x2 * (1 — (x2 + a21 * x1))
или же…
dИкс/ дт = р * Икс * (1 — * Икс)
(где K = a11 = a22 = 1)
Из того, что я понимаю, механизм расчета ODE Sundials ожидает объект класса ODE_vector, элементы которого представляют различные переменные состояния системы, в этом случае x1 и x2, и у меня нет проблем с получением решателя для аппроксимации траекторий двух видов ‘популяции, если я запрограммирую каждый параметр (r1, r2, a12, a21) как уникальные объекты, как в следующем коде:
int community::dynamics(ODE_vector const & state, ODE_vector & time_derivative) {
double r1 = 0.5;
double r2 = 0.2;
double a12 = 0.2;
double a21 = 0.2;
time_derivative[0] = r1 * state[0] * (1 - (state[0] + a12 * state[1]));
time_derivative[1] = r2 * state[1] * (1 - (state[1] + a21 * state[0]));
return 0;
};
… где time_derivative а также state являются членами класса Sundials ODE,
Может кто-нибудь показать мне, как адаптировать код, чтобы привести систему ODE в матричную форму, то есть, чтобы r1, r2, a12 и a21 были в форме р а также ? Мне сказали конвертировать из ODE_vector в векторные / матричные объекты броненосца, используя «Расширенные конструкторы» (например, 1) но это одно, чтобы знать, и другое, чтобы успешно реализовать это !!
Спасибо!
Решение
Где ты думаешь о чем-то подобном? Я наивно это реализовал за 30 мин. Оптимизации много, но она никогда не будет такой же производительной, как урезанный код. Даже если вы использовали Armadillo, Blitz ++ и т. Д.
#include <algorithm>
#include <iostream>
#include <vector>
template<class T> class Matrix {
public:
Matrix() {}
Matrix(size_t m, size_t n) : m_(m), n_(n) { v_.resize(m*n); }
Matrix& operator*= (const Matrix& other) {
if (other.m_ != n_) {
throw std::length_error("dimesions don't match");
}
Matrix<T> nv(m_,other.n_);
for (size_t j = 0; j < other.n_; ++j) {
for (size_t i = 0; i < m_; ++i) {
for (size_t k = 0; k < n_; ++k)
nv(i,j) += operator()(i,k) * other(k,j);
}
}
*this = nv;
return *this;
}
Matrix operator* (const Matrix& other) const {
Matrix ret = *this;
return ret *= other;
}
Matrix& operator+= (const Matrix& other) {
if (other.m_ != m_ || other.n_ != n_) {
throw std::length_error("dimesions don't match");
}
std::transform(v_.begin(), v_.end(), other.v_.begin(), v_.begin(), std::plus<T>());
return *this;
}
Matrix operator+ (const Matrix& other) const {
Matrix ret = *this;
return ret += other;
}
Matrix operator- () const {
Matrix<T> res (m_,n_);
std::transform (v_.begin(), v_.end(), res.v_.begin(), std::negate<T>());
return res;
}
Matrix& operator-= (const Matrix& other) {
return *this += -other;
}
Matrix operator- (const Matrix& other) const {
Matrix ret = *this;
return ret -= other;
}
Matrix operator->* (const Matrix& other) const {
if (other.m_ != m_ || other.n_ != n_) {
throw std::length_error("dimesions don't match");
}
Matrix res = *this;
std::transform(res.v_.begin(), res.v_.end(), other.v_.begin(),
res.v_.begin(), std::multiplies<T>());
return res;
}Matrix operator=(std::vector<std::initializer_list<T>> const& other) {
n_ = other.size();
if (n_ == 0) {
throw std::bad_alloc();
}
m_ = other.front().size();
if (m_ == 0) {
throw std::bad_alloc();
}
for (auto const& i : other) {
if (i.size() != m_) {
throw std::bad_alloc();
}
}
v_.resize(m_*n_);
size_t j = 0;
for (const auto& i : other) {
std::copy(i.begin(), i.end(), v_.begin()+j);
j+=m_;
}
return *this;
}
T& operator() (size_t m, size_t n) {
return v_[n*m_+m];
}
const T& operator() (size_t m, size_t n) const {
return v_[n*m_+m];
}size_t m() const {return m_;}
size_t n() const {return n_;}
private:
size_t m_, n_;
std::vector<T> v_;
};
template<class T> inline std::ostream&
operator<< (std::ostream& os, const Matrix<T>& v) {
size_t i,j;
for (i = 0; i < v.m(); ++i) {
for (j = 0; j < v.n(); ++j) {
os << v(i,j);
os << " ";
}
if (i<v.m()-1)
os << std::endl;
}
return os;
}
int main() {
Matrix<double> A, r, x, one, dxdt;
A = {{1,.2},
{.2,1}};
r = {{0.5,0.2}};
x = {{0.1,-0.1}};
one = {{1,1}};
dxdt = (r ->* x ->* (one - A * r));
std::cout << dxdt << std::endl;
return 0;
}
Другие решения
Других решений пока нет …