日晷ODE_vector变换为模逼近

Question

我试图使用Sundials ODE求解器库来近似于扩展到多物种Lotka竞争方程的动力学。例如在两个物种的情况下

dx1/dt = r1 * x1 * (1 - (x1 + a12 * x2))

dx2 2/dt= r2 * x2 * (1 - (x2 + a21 * x1))

或者..。

dx/dt = r * x * (1 - A * x)

(其中K= a11 = a22 = 1)

据我所知，Sundials ODE求解器期望一个类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的成员。

有人能告诉我如何修改代码，使之成为矩阵形式，即r1、r2、a12和a21是以r和A的形式存在的吗？我被告知使用“高级构造函数”(例如，ODE_vector )将从1转换为armadillo向量/矩阵对象，但要成功地实现这一点又是另一回事！！

谢谢!

Answer 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;
}