求解微分方程的数值解一般使用MATLAB等数值计算软件，其实C++也可以求解微分方程，需要用到odeint库，它是boost库的一部分。官方教程和示例比较晦涩，本文力求用较短的篇幅介绍它的基本用法,需要的朋友可以参考下面文章的具体内容

微分方程的标准形式为：


即：\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 代码实现

代码中有如下几个关键点：

  1. 要定义状态变量的类型state_type，定义为 std::vector<double> 即可
  2. 要用方程表示微分方程模型，和MATLAB中模型方程的写法非常类似
  3. 要写一个Observer以打印出计算结果，Observer函数也可以直接将数据写入文件中
  4. 要选择合适的求解器stepper，各种stepper的特点总结可以看 这里
  5. 要根据需要选择合适的集成函数，一般选择 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库求解微分方程内容请搜索编程学习网以前的文章希望大家以后多多支持编程学习网！