利用c++ boost库求解一阶微分方程计算船舶运动

船舶类的头文件 vessel.h

#pragma once
#include"Eigen/Dense"
#include"boost/numeric/odeint/stepper/runge_kutta_dopri5.hpp"

class Vessel {
public:
	static const size_t dim = 3;
	typedef Eigen::Vector3d state_type;
	struct VesselMotion {
		Eigen::Vector3d F;
		Eigen::Matrix3d D;
		Eigen::Matrix3d M_inv;
		VesselMotion(const Eigen::Matrix3d &M_inv, const Eigen::Matrix3d &D, Eigen::Vector3d &F) : M_inv(M_inv), D(D), F(F) {}
		void operator()(const state_type &x, state_type &dxdt, const double t) const {
			dxdt = M_inv * (-D * x + F);
		}
	};

	Vessel();
	void Move(Eigen::Vector3d &tauWind, Eigen::Vector3d &tauCurrent, Eigen::Vector3d &tauControl);

	Eigen::Vector3d eta;
	Eigen::Vector3d nu;
	Eigen::Vector3d nu_;
	Eigen::Matrix3d M;
	Eigen::Matrix3d D;
};

船舶类的成员函数 vessel.cpp

#include"vessel.h"

Vessel::Vessel() {
	M << 4, 0, 0,
		0, 4, 0,
		0, 0, 4;
	D << 0, 0, 0,
		0, 0, 0,
		0, 0, 0;
	eta << 0, 0, 0;
	nu << 0, 0, 0;
	nu_ << nu;
};

void Vessel::Move(Eigen::Vector3d &tauWind, Eigen::Vector3d &tauCurrent, Eigen::Vector3d &tauControl) {
	Eigen::Vector3d F = tauWind + tauCurrent + tauControl;
	const int steps = 10;
	const double dt = 0.1;
	state_type x = { this->nu[0], this->nu[1], this->nu[2] };
	boost::numeric::odeint::runge_kutta_dopri5<state_type> stepper;
	for (int step = 0; step < steps; ++step) {
		stepper.do_step(VesselMotion(this->M.inverse(), this->D, F), x, step*dt, dt);
	}
	this->nu_ << x[0], x[1], x[2];

	this->eta << (steps * dt / 2 * (this->nu + this->nu_));
	this->nu << this->nu_;
}

编写主程序求解 main.cpp

#include"vessel.h"

int main() {
	Vessel vessel;
	Eigen::Vector3d tauWind;
	Eigen::Vector3d tauCurrent;
	Eigen::Vector3d tauControl;
	tauWind << 1, 1, 1;
	tauCurrent << 1, 1, 1;
	tauControl << 4, 4, 4;
	vessel.Move(tauWind, tauCurrent, tauControl);
	
	system("pause");
	return 0;
}

其中，Eigen库和boost库在visual studio 2017 中的添加参见

Eigen库在VS2017下的配置与使用_Kerwines的博客-CSDN博客_vs配置eigen

VS配置boost库_今天也要debug的博客-CSDN博客_vs安装boost库