用C++的odeint库求解微分方程
微分方程的标准形式为：
即：\dot{\boldsymbol{x}} = \boldsymbol{f}(\boldsymbol{x}, t),\, \boldsymbol{x}(0) = \boldsymbol{x_0}
这是一阶微分方程组， \boldsymbol{x} 和 \boldsymbol{f}(\boldsymbol{x}, t) 均为向量。如果要求解高阶微分方程，需要先转换成一阶微分方程组后再用odeint求解。
1、集成方程
API中最重要的是集成函数（integrate functions），一共有5种，它们的调用接口很类似。 integrate_const 的函数调用方式为：
integrate_const(stepper, system, x0, t0, t1, dt, observer)
其中：
stepper是求解器，也就是所使用的数值算法（例如Runge-Kutta或Euler法）system是待求解的微分方程x0是初始条件t0和 t1 分别是初始时间和结束时间dt是时间间隔，它重要与否取决于求解器的类型observer是每N个时间间隔调用一次的函数，可用来打印实时的解，该参数是可选的，如果没有此参数，集成函数会从 t0 计算到 t1 ，不产生任何输出就返回
给定初始状态 x0 ，集成函数从初始时间 t0 到结束时间 t1 不断地调用给定的 stepper ，计算微分方程在不同时刻的解，用户还可以提供 observer 以分析某个时刻的状态值。具体选择哪个集成函数取决于你想要什么类型的结果，也就是调用 observer 的频次。
integrate_const 每过相等的时间间隔 dt 会调用一次 observer ，语法为：
integrate_const(stepper, system, x0, t0, t1, dt, observer)
integrate_n_steps 和前面的类似，但它不需要知道结束时间，它只需要知道要计算的步数，语法为：
integrate_n_steps(stepper, system, x0, t0, dt, n, observer)
integrate_times 计算在用户给定时间点的值，语法为：
integrate_times(stepper, system, x0, times_start, times_end, dt, observer) integrate_times(stepper, system, x0, time_range, dt, observer)
integrate_adaptive 用于需要在每个时间间隔调用 observer 的场合，语法为：
integrate_adaptive(stepper, system, x0, t0, t1, dt, observer)
integrate 是最方便的集成函数, 不需要指定 stepper ，简单快捷，语法为：
integrate(system, x0, t0, t1, dt, observer)
求解器stepper的选择（比如自适应方式会根据误差修改时间间隔）会改变计算的具体实现方式, 但是observer的调用（也就是你的输出结果）依然遵循上述规则。
2、求解单摆模型
2.1 微分方程标准化
现在求单摆系统微分方程的解，以得出单摆角度随时间变化的规律。其微分方程
即：\ddot{\theta}(t) = -\mu \dot{\theta}(t) - \frac{g}{L} \sin \theta(t)
即：\begin{cases} \dot{\theta}(t) & = \omega(t) \\ \dot{\omega}(t) & = -\mu \omega(t) - g/L \sin \theta(t) \end{cases}
令状态变量
即：\boldsymbol{x} = \begin{bmatrix} x_1(t)\\ x_2(t) \end{bmatrix} = \begin{bmatrix} \theta(t)\\ \omega(t) \end{bmatrix}
微分方程组变为
即：\frac{\mathrm{d}\boldsymbol{x}}{\mathrm{d}t}= \frac{\mathrm{d}}{\mathrm{d}t} \begin{bmatrix} x_1(t)\\ x_2(t) \end{bmatrix} = \begin{bmatrix} x_2(t)\\ -\mu x_2(t) - g/L \sin x_1(t) \end{bmatrix}
2.2 代码实现
代码中有如下几个关键点：
- 要定义状态变量的类型
state_type，定义为std::vector<double>即可 - 要用方程表示微分方程模型，和
MATLAB中模型方程的写法非常类似 - 要写一个
Observer以打印出计算结果，Observer函数也可以直接将数据写入文件中 - 要选择合适的求解器
stepper，各种stepper的特点总结可以看 这里 - 要根据需要选择合适的集成函数，一般选择
integrate_const即可满足要求
下面的代码可作为标准模板使用：
#include <iostream>
#include <cmath>
#include <boost/numeric/odeint.hpp>
using namespace std;
using namespace boost::numeric::odeint;
const double g = 9.81; // 重力加速度
const double L = 1.00; // 摆线长度
const double mu = 0.80; // 阻力系数
// 定义状态变量的类型
typedef std::vector<double> state_type;
// 要求解的微分方程
void pendulum(const state_type &x, state_type &dxdt, double t)
{
dxdt[0] = x[1];
dxdt[1] = -mu*x[1] - g/L*sin(x[0]);
}
// Observer打印状态值
void write_pendulum(const state_type &x, const double t)
{
cout << t << '\t' << x[0] << '\t' << x[1] << endl;
}
int main(int argc, char **argv)
{
// 初始条件，二维向量
state_type x = {0.10 , 0.00};
// 求解方法为runge_kutta4
integrate_const(runge_kutta4<state_type>(), pendulum, x , 0.0 , 5.0 , 0.01 , write_pendulum);
}
编译该程序依赖boost库，要在 CMakeLists.txt 中添加相应的内容。编译成功后运行，会得到如下的结果：
0 0.1 0 0.01 0.0999512 -0.009753 0.02 0.0998052 -0.0194188 0.03 0.0995631 -0.0289887 0.04 0.0992258 -0.0384542 0.05 0.0987944 -0.0478069 0.06 0.0982701 -0.0570385 0.07 0.0976541 -0.0661412 0.08 0.0969477 -0.075107 0.09 0.0961524 -0.0839283 0.1 0.0952696 -0.0925977 0.11 0.094301 -0.101108 ---- many lines ommitted ----
可以将输出数据重定向到文本文件 data.txt 中，然后使用Python等脚本语言提取数据并画图显示。下面是实现该功能的参考代码：
import numpy as np
import matplotlib.pyplot as plt
lines = tuple(open("data.txt", 'r')) # 读取文件行到tuple中
rows = len(lines)
time = np.zeros(rows)
theta = np.zeros(rows)
omega = np.zeros(rows)
for r in range(rows):
[str1, str2, str3] = lines[r].split()
time[r] = float(str1)
theta[r] = float(str2)
omega[r] = float(str3)
plt.plot(time, theta, time, omega) # 角度和角速度变化
# plt.plot(theta, omega) # 相图
plt.show()
到此这篇关于用C++的odeint库求解微分方程的文章就介绍到这了,更多相关C++的odeint库求解微分方程内容请搜索脚本之家以前的文章或继续浏览下面的相关文章希望大家以后多多支持脚本之家！
相关文章
这篇文章主要为大家详细介绍了OpenCV实现霍夫变换直线检测，文中示例代码介绍的非常详细，具有一定的参考价值，感兴趣的小伙伴们可以参考一下2021-06-06
这篇文章主要介绍了Windows程序内部运行机制实例详解,对于学习Windows程序设计来说是非常重要的一课,需要的朋友可以参考下2014-08-08
这篇文章主要介绍了C语言中变长数组的使用，文中通过示例代码介绍的非常详细，对大家的学习或者工作具有一定的参考学习价值，需要的朋友们下面随着小编来一起学习学习吧2023-02-02
这篇文章主要为大家介绍了C C++ 题解LeetCode1417重新格式化字符串，有需要的朋友可以借鉴参考下，希望能够有所帮助，祝大家多多进步，早日升职加薪2022-10-10
这篇文章主要为大家详细介绍了C++实现简单的职工工资管理系统，文中示例代码介绍的非常详细，具有一定的参考价值，感兴趣的小伙伴们可以参考一下2022-03-03
这篇文章主要为大家详细介绍了用C语言实现通讯录，文中示例代码介绍的非常详细，具有一定的参考价值，感兴趣的小伙伴们可以参考一下2022-06-06
这篇文章主要介绍了C++虚函数的工作原理的的相关资料，文中讲解非常细致，代码帮助大家更好的理解和学习，感兴趣的朋友可以了解下2020-06-06
本文主要介绍了当给你一个字符串s，其中包含字母顺序打乱的用英文单词表示的若干数字（0-9）时，如何通过Leetcode按升序返回原始的数字。感兴趣的童鞋可以来看看2021-11-11
这篇文章主要为大家详细介绍了C++基于CMD命令行实现扫雷小游戏，文中示例代码介绍的非常详细，具有一定的参考价值，感兴趣的小伙伴们可以参考一下2022-05-05
C++ WideCharToMultiByte()函数案例详解
这篇文章主要介绍了C++ WideCharToMultiByte()函数案例详解,本篇文章通过简要的案例,讲解了该项技术的了解与使用,以下就是详细内容,需要的朋友可以参考下2021-08-08
最新评论