三维微分方程组四阶龙格库塔求解

河北一帆

已于 2022-10-26 10:02:02 修改

阅读量729

点赞数 1

文章标签： c++ 算法开发语言

于 2022-10-26 09:36:56 首次发布

本文链接：https://blog.csdn.net/qq_41816368/article/details/127525958

版权

利用c++ boost库求解一阶微分方程计算船舶运动_河北一帆的博客-CSDN博客_c++求解微分方程

这篇博文中利用boost库求解线性三自由度船舶运动模型。但是boost库非常庞大，达到好几个G，不如自己写一个龙格库塔求解。网上搜索到的例子大多数是一维微分方程的，推广到三维其实很容易。三自由度船舶运动模型如下所示。

一维微分方程

$\frac{dx}{dt}=f(t,x)$

$x(t_0)=x_0$

四阶龙格库塔计算公式为

$K1=f(t_i,x_i)$

$K2=f(t_i+\frac{h}{2},x_i+\frac{hK1}{2})$

$K3=f(t_i+\frac{h}{2},x_i+\frac{hK2}{2})$

$K4=f(t_i+h,x_i+hK3)$

$x_{i+1}=x_{i}+h\frac{K1+2K2+2K3+K4}{6}$

$h$ 是求解步长。对于推广到三维，直接将三维向量 $\nu$ 替换 $x$

double h = 0.1;
double dt = 1;
Eigen::Vector3d nu_;
nu_ << this->nu[0], this->nu[1], this->nu[2];
for (int i = 0; i < 10; i++) {
    Eigen::Vector3d K1 = this->f(nu_, tau);
    Eigen::Vector3d K2 = this->f(nu_ + 0.5 * h * K1, tau);
    Eigen::Vector3d K3 = this->f(nu_ + 0.5*h*K2, tau);
    Eigen::Vector3d K4 = this->f(nu_ + h * K3, tau);
    nu_ = nu_ + h / 6 * (K1 + 2 * K2 + 2 * K3 + K4);
}

Eigen::Vector3d mss::Vessel::f(Eigen::Vector3d y, Eigen::Vector3d F) {
    return this->M.inverse() * (-this->D * y + F);
}