用C++的odeint库求解微分方程

发表于 2021-09-17 分类于程序设计 Waline：本文字数： 1.4k 阅读时长 ≈ 5 分钟

介绍

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

所求解的微分方程的标准形式为 \[ \dot{\mathbf{x}} = \mathbf{f}(\mathbf{x}, t),\, \mathbf{x}(0) = \mathbf{x_0} \]

这是一阶微分方程组，\(\mathbf{x}\)和\(\mathbf{f}(\mathbf{x}, t)\)均为向量。高阶微分方程要先转换成一阶微分方程组，然后再用odeint求解。

integrate函数

API中最重要的是integrate函数，一共有5种，它们的调用接口很类似。integrate_const的函数调用方式为：

1	integrate_const(stepper, system, x0, t0, t1, dt, observer)

其中：

stepper是求解器，也就是所使用的数值算法（例如Runge-Kutta法、Euler法）
system是待求解的微分方程
x0是初始条件
t0和t1分别是初始时间和结束时间
dt是时间间隔，它重要与否取决于求解器的类型
observer是每N个时间间隔调用一次的函数，可用来打印实时的解，该参数是可选的，如果没有此参数，integrate函数会从t0计算到t1，不产生任何输出就返回

给定初始状态x0，integrate函数从初始时间t0到结束时间t1不断地调用给定的stepper，计算微分方程在不同时刻的解，用户还可以提供observer以分析某个时刻的状态值。具体选择哪个integrate函数取决于你想要什么类型的结果，也就是调用observer的频率。

integrate_const每过相等的时间间隔dt会调用一次observer，语法为：

1	integrate_const(stepper, system, x0, t0, t1, dt, observer)

integrate_n_steps和前面的类似，但它不需要知道结束时间，它只需要知道要计算的步数，语法为：

1	integrate_n_steps(stepper, system, x0, t0, dt, n, observer)

integrate_times计算用户给定时间点上的值，语法为：

1 2	integrate_times(stepper, system, x0, times_start, times_end, dt, observer) integrate_times(stepper, system, x0, time_range, dt, observer)

integrate_adaptive用于需要在每个时间间隔都调用observer的场合，语法为：

1	integrate_adaptive(stepper, system, x0, t0, t1, dt, observer)

integrate最方便, 不需要指定stepper，简单快捷，语法为：

1	integrate(system, x0, t0, t1, dt, observer)

求解器stepper的选择（比如自适应方式会根据误差修改时间间隔）会改变计算的具体实现方式, 但是observer的调用（也就是输出结果）依然遵循上述规则。

求解单摆模型

微分方程标准化

现在求单摆系统微分方程的解，以得到单摆角度随时间变化的规律。其微分方程 \[ \ddot{\theta}(t) = -\mu \dot{\theta}(t) - \frac{g}{L} \sin \theta(t) \]

\(\theta(t)\)是摆角随时间变化的规律，\(\mu\)是空气阻力系数，\(g\)是重力加速度，\(L\)是摆球到悬挂点的长度。首先将它转化成一阶微分方程组，新增一个变量角速度\(\omega(t) = \dot{\theta}(t)\)，方程变为 \[ \begin{align*} \dot{\theta}(t) & = \omega(t) \\ \dot{\omega}(t) & = -\mu \omega(t) - g \sin \theta(t) / L \end{align*} \]

令状态变量 \[ \mathbf{x} = \begin{bmatrix} x_1(t)\\ x_2(t) \end{bmatrix} = \begin{bmatrix} \theta(t)\\ \omega(t) \end{bmatrix} \]

微分方程组变为 \[ \dot{\mathbf{x}}= \begin{bmatrix} \dot{x}_1(t)\\ \dot{x}_2(t) \end{bmatrix} =\mathbf{f}(\mathbf{x}, t)= \begin{bmatrix} x_2(t)\\ -\mu x_2(t) - g \sin x_1(t) / L \end{bmatrix} \]

代码实现

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

要定义状态变量的类型，将state_type简单定义为std::vector<double>即可
要用方程表示微分方程模型，和MATLAB中模型方程的写法非常类似
要写一个Observer以打印出计算结果，Observer函数也可以直接将数据写入文件中
要选择合适的求解器stepper，各种stepper的特点总结可以看这里
要根据需要选择合适的integrate函数，一般选择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()