用C++的odeint庫求解微分方程
微分方程的標(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 分別是初始時間和結(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 微分方程標(biāo)準(zhǔn)化
現(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{\mathrmwppm3vysvbp\boldsymbol{x}}{\mathrmwppm3vysvbpt}= \frac{\mathrmwppm3vysvbp}{\mathrmwppm3vysvbpt} \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即可滿足要求
下面的代碼可作為標(biāo)準(zhǔn)模板使用：
#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)文章
這篇文章主要介紹了編輯器寫C語言輸出中文亂碼問題及解決方案,具有很好的參考價值,希望對大家有所幫助,如有錯誤或未考慮完全的地方,望不吝賜教2023-09-09
這篇文章主要為大家詳細(xì)介紹了在Qt中如何實現(xiàn)讀取和寫入配置文件的數(shù)據(jù)，文中的示例代碼講解詳細(xì)，具有一定的學(xué)習(xí)價值，感興趣的小伙伴可以了解一下2023-03-03
linux內(nèi)核select/poll，epoll實現(xiàn)與區(qū)別
這篇文章主要介紹了linux內(nèi)核select/poll，epoll實現(xiàn)與區(qū)別,需要的朋友可以參考下2016-11-11
這篇文章主要為大家介紹了C語言軟件spi虛擬總線中間層設(shè)計詳解，有需要的朋友可以借鑒參考下，希望能夠有所幫助，祝大家多多進步，早日升職加薪2023-01-01
大家好，本篇文章主要講的是C語言實現(xiàn)二叉樹層次遍歷介紹，感興趣的同學(xué)趕快來看一看吧，對你有幫助的話記得收藏一下2022-01-01
C++構(gòu)造函數(shù)初始化列表的實現(xiàn)詳解
構(gòu)造函數(shù)主要作用在于創(chuàng)建對象時為對象的成員屬性賦值，構(gòu)造函數(shù)由編譯器自動調(diào)用,無須手動調(diào)用；析構(gòu)函數(shù)主要作用在于對象銷毀前系統(tǒng)自動調(diào)用，執(zhí)行一 些清理工作2022-09-09
C語言數(shù)據(jù)結(jié)構(gòu)之二分法查找詳解
二分查找算法是在有序數(shù)組中用到的較為頻繁的一種算法，在未接觸二分查找算法時，最通用的一種做法是，對數(shù)組進行遍歷，跟每個元素進行比較，其時間為O(n)，但二分查找算法更優(yōu)2022-02-02
這篇文章主要介紹了C++?std::thread?如何使用,C++中的std::thread類提供了一種方便的多線程編程方式，在使用std::thread類時，我們需要注意線程間的同步和通信問題，以確保多個線程之間的正確協(xié)同工作需要的朋友可以參考下2023-03-03
FFmpeg實戰(zhàn)之分離出PCM數(shù)據(jù)
PCM(Pulse?Code?Modulation，脈沖編碼調(diào)制)音頻數(shù)據(jù)是未經(jīng)壓縮的音頻采樣數(shù)據(jù)裸流，它是由模擬信號經(jīng)過采樣、量化、編碼轉(zhuǎn)換成的標(biāo)準(zhǔn)數(shù)字音頻數(shù)據(jù)。本文將通過FFmpeg實現(xiàn)分離PCM數(shù)據(jù)，感興趣的可以了解一下2023-02-02
這篇文章主要介紹了使用單鏈表實現(xiàn)多項式計算示例,需要的朋友可以參考下2014-03-03