用C++的odeint庫求解微分方程
微分方程的標準形式為：
即：\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 分別是初始時間和結(jié)束時間dt是時間間隔，它重要與否取決于求解器的類型observer是每N個時間間隔調(diào)用一次的函數(shù)，可用來打印實時的解，該參數(shù)是可選的，如果沒有此參數(shù)，集成函數(shù)會從 t0 計算到 t1 ，不產(chǎn)生任何輸出就返回
給定初始狀態(tài) x0 ，集成函數(shù)從初始時間 t0 到結(jié)束時間 t1 不斷地調(diào)用給定的 stepper ，計算微分方程在不同時刻的解，用戶還可以提供 observer 以分析某個時刻的狀態(tài)值。具體選擇哪個集成函數(shù)取決于你想要什么類型的結(jié)果，也就是調(diào)用 observer 的頻次。
integrate_const 每過相等的時間間隔 dt 會調(diào)用一次 observer ，語法為：
integrate_const(stepper, system, x0, t0, t1, dt, observer)
integrate_n_steps 和前面的類似，但它不需要知道結(jié)束時間，它只需要知道要計算的步數(shù)，語法為：
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 用于需要在每個時間間隔調(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)方式會根據(jù)誤差修改時間間隔）會改變計算的具體實現(xiàn)方式, 但是observer的調(diào)用（也就是你的輸出結(jié)果）依然遵循上述規(guī)則。
2、求解單擺模型
2.1 微分方程標準化
現(xiàn)在求單擺系統(tǒng)微分方程的解，以得出單擺角度隨時間變化的規(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>
即：\frac{\mathrmvvxyksv9kd\boldsymbol{x}}{\mathrmvvxyksv9kdt}= \frac{\mathrmvvxyksv9kd}{\mathrmvvxyksv9kdt} \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 代碼實現(xiàn)
代碼中有如下幾個關(guān)鍵點：
- 要定義狀態(tài)變量的類型
state_type，定義為std::vector<double>即可 - 要用方程表示微分方程模型，和
MATLAB中模型方程的寫法非常類似 - 要寫一個
Observer以打印出計算結(jié)果，Observer函數(shù)也可以直接將數(shù)據(jù)寫入文件中 - 要選擇合適的求解器
stepper，各種stepper的特點總結(jié)可以看 這里 - 要根據(jù)需要選擇合適的集成函數(shù)，一般選擇
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; // 阻力系數(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);
}
編譯該程序依賴boost庫，要在 CMakeLists.txt 中添加相應(yīng)的內(nèi)容。編譯成功后運行，會得到如下的結(jié)果：
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 ----
可以將輸出數(shù)據(jù)重定向到文本文件 data.txt 中，然后使用Python等腳本語言提取數(shù)據(jù)并畫圖顯示。下面是實現(xiàn)該功能的參考代碼：
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()
到此這篇關(guān)于用C++的odeint庫求解微分方程的文章就介紹到這了,更多相關(guān)C++的odeint庫求解微分方程內(nèi)容請搜索腳本之家以前的文章或繼續(xù)瀏覽下面的相關(guān)文章希望大家以后多多支持腳本之家！
相關(guān)文章
這篇文章主要為大家詳細介紹了opencv2實現(xiàn)10張圖像上下左右拼接融合，文中示例代碼介紹的非常詳細，具有一定的參考價值，感興趣的小伙伴們可以參考一下2020-03-03
這篇文章主要介紹了c語言 字符串轉(zhuǎn)大寫的簡單實例，有需要的朋友可以參考一下2013-12-12
淺析C/C++中動態(tài)鏈接庫的創(chuàng)建和調(diào)用
下面小編就為大家?guī)硪黄獪\析C/C++中動態(tài)鏈接庫的創(chuàng)建和調(diào)用。小編覺得挺不錯的，現(xiàn)在分享給大家，也給大家做個參考，一起跟隨小編過來看看吧2016-05-05
所謂嵌套(Nest),就是一條語句里面還有另一條語句,例如 for 里面還有 for,while 里 面還有 while,或者 for 里面有 while,while 里面有 if-else,這都是允許的2022-05-05
內(nèi)存對齊的目的是為了提高CPU讀寫內(nèi)存里數(shù)據(jù)的速度?，F(xiàn)代的CPU讀取內(nèi)存并不是一個一個字節(jié)挨著讀取，這樣做的效率非常低?，F(xiàn)代的CPU一般以4個字節(jié)(32bit數(shù)據(jù)總線）或者8個字節(jié)（64bit數(shù)據(jù)總線）為一組，一組一組地讀寫內(nèi)存里的數(shù)據(jù)2021-08-08
這篇文章主要為大家詳細介紹了C語言實現(xiàn)通訊錄程序，文中示例代碼介紹的非常詳細，具有一定的參考價值，感興趣的小伙伴們可以參考一下2021-10-10
這篇文章主要介紹了C語言怎么獲得進程的PE文件信息的相關(guān)代碼,需要的朋友可以參考下2016-01-01
這篇文章主要介紹了史上最強C語言分支和循環(huán)教程詳解,本文通過代碼演示給大家介紹的非常詳細，對大家的學(xué)習(xí)或工作具有一定的參考借鑒價值，需要的朋友可以參考下2021-11-11
本篇文章是對C語言實現(xiàn)shell指令的方法進行了詳細的分析介紹，需要的朋友參考下2013-05-05
這篇文章主要介紹了C++如何實現(xiàn)字符串的部分復(fù)制問題,具有很好的參考價值,希望對大家有所幫助,如有錯誤或未考慮完全的地方,望不吝賜教2023-08-08