微分方程的標(biāo)準(zhǔn)形式為：

即：\dot{\boldsymbol{x}} = \boldsymbol{f}(\boldsymbol{x}, t),\, \boldsymbol{x}(0) = \boldsymbol{x_0}

這是一階微分方程組， \boldsymbol{x} 和 \boldsymbol{f}(\boldsymbol{x}, t) 均為向量。如果要求解高階微分方程，需要先轉(zhuǎn)換成一階微分方程組后再用odeint求解。

1、集成方程

API中最重要的是集成函數(shù)（integrate functions），一共有5種，它們的調(diào)用接口很類似。 integrate_const 的函數(shù)調(diào)用方式為：

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

其中：

• stepper 是求解器，也就是所使用的數(shù)值算法（例如Runge-Kutta或Euler法）
• system 是待求解的微分方程
• x0 是初始條件
• t0 和 t1 分別是初始時(shí)間和結(jié)束時(shí)間
• dt 是時(shí)間間隔，它重要與否取決于求解器的類型
• observer 是每N個(gè)時(shí)間間隔調(diào)用一次的函數(shù)，可用來打印實(shí)時(shí)的解，該參數(shù)是可選的，如果沒有此參數(shù)，集成函數(shù)會(huì)從 t0 計(jì)算到 t1 ，不產(chǎn)生任何輸出就返回
  

給定初始狀態(tài) x0 ，集成函數(shù)從初始時(shí)間 t0 到結(jié)束時(shí)間 t1 不斷地調(diào)用給定的 stepper ，計(jì)算微分方程在不同時(shí)刻的解，用戶還可以提供 observer 以分析某個(gè)時(shí)刻的狀態(tài)值。具體選擇哪個(gè)集成函數(shù)取決于你想要什么類型的結(jié)果，也就是調(diào)用 observer 的頻次。

integrate_const 每過相等的時(shí)間間隔 dt 會(huì)調(diào)用一次 observer ，語法為：

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

integrate_n_steps 和前面的類似，但它不需要知道結(jié)束時(shí)間，它只需要知道要計(jì)算的步數(shù)，語法為：

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

integrate_times 計(jì)算在用戶給定時(shí)間點(diǎn)的值，語法為：

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

integrate_adaptive 用于需要在每個(gè)時(shí)間間隔調(diào)用 observer 的場合，語法為：

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

integrate 是最方便的集成函數(shù), 不需要指定 stepper ，簡單快捷，語法為：

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

求解器stepper的選擇（比如自適應(yīng)方式會(huì)根據(jù)誤差修改時(shí)間間隔）會(huì)改變計(jì)算的具體實(shí)現(xiàn)方式, 但是observer的調(diào)用（也就是你的輸出結(jié)果）依然遵循上述規(guī)則。

2、求解單擺模型

2.1 微分方程標(biāo)準(zhǔn)化

現(xiàn)在求單擺系統(tǒng)微分方程的解，以得出單擺角度隨時(shí)間變化的規(guī)律。其微分方程

即：\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}


令狀態(tài)變量

即：\boldsymbol{x} = \begin{bmatrix} x_1(t)\\ x_2(t) \end{bmatrix} = \begin{bmatrix} \theta(t)\\ \omega(t) \end{bmatrix}

微分方程組變?yōu)?/strong>

#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; // 阻力系數(shù)

// 定義狀態(tài)變量的類型
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打印狀態(tài)值
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);
}

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    ----

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()