將剛性ODE與Python集成

我正在尋找一個好的庫，它將在Python中集成剛性ODE。問題是，scipy的odeint給我提供了很好的解決方案有時候是，但是初始條件稍有變化就會導致它掉下來並放棄。 MATLAB的僵硬求解器（ode15s和ode23s）很快樂地解決了同樣的問題，但我無法使用它（即使是從Python開始，因爲沒有一個用於MATLAB C API的Python綁定實現了回調函數，而且我需要傳遞一個函數到ODE求解器）。我正在嘗試PyGSL，但它非常複雜。任何建議將不勝感激。將剛性ODE與Python集成

編輯：我與PyGSL的具體問題是選擇正確的步驟功能。其中有幾個，但沒有直接類似於ode15s或ode23s（bdf公式和修改的Rosenbrock如果有意義的話）。那麼選擇僵硬系統有什麼好的步驟功能？爲了確保它達到穩定狀態，我必須解決這個系統很長時間，GSL解算器要麼選擇一個微小的時間步長，要麼選擇一個太大的時間步長。

來源

2010-01-18 Chinmay Kanchi

我想用PyGSL來幫助你。我從來沒有用過它，但我有GSL的經驗。我只是看了pygsl（odeiv.py）中提供的示例，它看起來與C中的幾乎相同。您是否認爲PyGSL由於python接口或GSL本身而非常複雜？ – YuppieNetworking 2010-01-18 20:04:09

好吧，可怕的複雜也許是多餘的:)。它比MATLAB或scipy的複雜度要高出一個數量級。爲了澄清，python接口與C接口幾乎相同，所以它的庫本身很複雜。另外，PyGSL不會記錄odeiv，所以我必須使用C文檔來找出在Python中要做什麼。不好玩。 – 2010-01-18 20:07:14

編輯該問題。 – 2010-01-18 20:09:29

Python可以調用C.ODEPACK中的行業標準是LSODE。它是公共領域。您可以下載C version。這些求解器非常棘手，所以最好使用一些經過良好測試的代碼。

補充：確保你確實有一個僵硬的系統，即如果比率（特徵值）相差超過2或3個數量級。另外，如果系統是僵硬的，但你只是尋找一個穩態解決方案，這些求解器給你選擇以代數方式求解一些方程。否則，像DVERK這樣的Runge-Kutta求解器將會是一個很好且更簡單的解決方案。

在這裏添加，因爲它不會在評論適合：這是從DLSODE頭DOC：

C  T  :INOUT Value of the independent variable. On return it 
C     will be the current value of t (normally TOUT). 
C 
C  TOUT :IN  Next point where output is desired (.NE. T).

而且，是米氏動力學是非線性的。儘管如此，艾特肯加速技術也適用於它。（如果你想做一個簡短的解釋，首先考慮Y是一個標量的簡單情況，你運行系統得到3個Y（T）點，通過它們擬合一個指數曲線（簡單代數），然後將Y設置爲漸近線，重複一次，現在推廣到Y是一個向量，假設3個點在一個平面上 - 如果它們不是的話，這是可以的）。除此之外，除非你有一個強制函數（如常量IV滴），否則MM消除會衰減遠離系統將接近線性。希望有所幫助。

來源

2010-01-20 19:58:01

感謝LSODE提示。我發現有人已經在http://www.cs.uiuc.edu/homes/mrgates2/ode/上爲它編寫了一個Python接口。明天試試看，如果有效，請接受你的答案。 – 2010-01-20 22:32:23

@Chinmay：賓果！希望對你有效。 – 2010-01-20 23:03:11

我知道我確實有一個僵硬的系統，但是看着LSODE源代碼，我意識到我的程序中有一個可能的錯誤。 scipy.integrate.odeint求解器基於LSODA，和所有LS *求解器一樣，它需要一個時間向量來計算解。原來，我計算這個時間向量太粗糙了。我從MATLAB切換到Python，並認爲numpy.arange會在MATLAB中類似於i = t0：tn。它沒有。所以我一直在整個時間只求解整數... – 2010-01-21 00:33:16

目前我正在一點ODE及其求解的，所以你的問題對我來說很有趣......

從我所聽到和讀到，對於剛性問題的正確方法去是選擇一個隱含的方法作爲一個階梯函數（糾正我，如果我錯了，我仍然在學習ODE求解器的魔法）。我不能引用你在哪裏閱讀這篇文章，因爲我不記得，但這裏有一個來自gsl-help的thread，那裏有類似的問題。

因此，總之，似乎bsimp方法值得一試，儘管它需要一個雅可比函數。如果您無法計算雅可比行列式，我會嘗試使用rk2imp,rk4imp或任何齒輪方法。

來源

2010-01-19 11:20:03 YuppieNetworking

如果你可以通過Matlab的ode15s解決你的問題，你應該可以通過求解scipy的vode求解。爲了模擬ode15s，我使用以下設置：

ode15s = scipy.integrate.ode(f) 
ode15s.set_integrator('vode', method='bdf', order=15, nsteps=3000) 
ode15s.set_initial_value(u0, t0)

，然後你可以愉快地ode15s.integrate(t_final)解決您的問題。它應該在一個棘手的問題上工作得很好。

（也參見http://www.scipy.org/NumPy_for_Matlab_Users）

來源

2010-02-11 11:01:28

好！感謝那。最後我用'vode'或多或少地得出了一個類似的解決方案。 – 2010-02-11 11:07:43

'''bdf'''的最大訂單是5.無論如何，沒有意義設置高於12的訂單。因爲12是亞當斯的最大值。 – 2016-07-22 15:05:05

PyDSTool包裹10期Radau解算器，它是一個極好的隱式僵硬積分器。這比odeint有更多的設置，但比PyGSL少很多。最大的好處是你的RHS函數被指定爲一個字符串（通常，雖然你可以使用符號操作來構建一個系統）並且被轉換爲C，所以沒有緩慢的python回調，整個事情會非常快。

來源

2011-07-26 19:38:41 RHC

將剛性ODE與Python集成

回答

相關問題