0%

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

介绍

求解微分方程的数值解一般使用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求解。

集成方程

API中最重要的是集成函数(integrate functions),一共有5种,它们的调用接口很类似。integrate_const的函数调用方式为:

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

其中:

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

给定初始状态x0,集成函数从初始时间t0到结束时间t1不断地调用给定的stepper,计算微分方程在不同时刻的解,用户还可以提供observer以分析某个时刻的状态值。具体选择哪个集成函数取决于你想要什么类型的结果,也就是调用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/L \sin \theta(t) \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/L \sin x_1(t) \end{bmatrix} \]

代码实现

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

  1. 要定义状态变量的类型state_type,定义为std::vector<double>即可
  2. 要用方程表示微分方程模型,和MATLAB中模型方程的写法非常类似
  3. 要写一个Observer以打印出计算结果,Observer函数也可以直接将数据写入文件中
  4. 要选择合适的求解器stepper,各种stepper的特点总结可以看这里
  5. 要根据需要选择合适的集成函数,一般选择integrate_const即可满足要求

下面的代码可作为标准模板使用:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
#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中添加相应的内容。编译成功后运行,会得到如下的结果:

1
2
3
4
5
6
7
8
9
10
11
12
13
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等脚本语言提取数据并画图显示。下面是实现该功能的参考代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
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()

参考

  1. Overview
  2. Integrate funcitons
  3. Numerical integration in C++ with Boost odeint
  4. Odeint library