1

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

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

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

或者...

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

（其中 K = a11 = a22 = 1）

据我了解，日晷 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是日晷班的成员ODE。

谁能告诉我如何修改代码以将 ODE 系统转换为矩阵形式，即 r1、r2、a12 和 a21 的形式为r和A？我被告知要使用“高级构造函数”（例如1）从 ODE_vector 转换为犰狳矢量/矩阵对象，但知道这一点是一回事，而要成功实现它又是另一回事！

谢谢！

4

1 回答 1

0

你在哪里想到这样的事情？我在 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;
}
于 2018-07-06T08:57:27.817 回答