Python的科學：中斷微分方程的條件

我使用odeint模擬帶電粒子在電場（來源來源於此package）目前求解微分方程的制度下蟒蛇解決：Python的科學：中斷微分方程的條件

time = np.linspace(0, 5, 1000) 

def sm(x, t): 
    return np.array([x[1], eta*Ez0(x[0])]) 

traj = odeint(sm,[0,1.], time)

它工作正常，但我想盡快停止計算爲x [0] < 0。目前我只是阻止系統正演變：

def sm1(x, t): 
    if x[0] < 0: 
     return np.array([0, 0]) 
    else: 
     return np.array([x[1], eta*Ez0(x[0])]) 

traj = odeint(sm1,[0,1.],time)

但我GESS有更好的解決方案。我發現this但在我看來，它修復了步驟的數量，這是令人遺憾的。任何建議表示讚賞。

來源

2010-11-19 Mermoz

您的解決方案對我來說似乎很合理。你認爲什麼是錯的？ – 2010-11-19 17:09:47

你會得到一個警告：lsoda--在當前t（= r1），mxstep（= i1）這個呼叫在到達tout之前採取的步驟在上面的消息中，I1 = 500 在上面的消息中，R1 = 0.4223349048304E + 00 此次調用完成了過多的工作（可能是錯誤的Dfun類型）。以full_output = 1運行以獲取定量信息。 – Mermoz 2010-11-19 17:12:55

如果您編寫odeint函數的自定義擴展，您可以讓函數在完成時引發特定的異常。用Python做它可能會讓它慢得多，但我認爲你用C或Cython編寫了相同的東西。請注意，我沒有測試以下內容。

class ThatsEnoughOfThat(Exception): 
    pass 

def custom_odeint(func, y0, t): # + whatever parameters you need 
    for timestep in t: 
     try: 
      # Do stuff. Call odeint/other scipy functions? 
     except ThatsEnoughOfThat: 
      break 
    return completedstuff 

def sm2(x, t): 
    if x[0] < 0: 
     raise ThatsEnoughOfThat 
    return np.array([x[1], eta*Ez0(x[0])])

來源

2010-11-19 18:06:53

Python的科學：中斷微分方程的條件

回答

相關問題