SciPy的：實施微分方程的兩種方法：兩種不同的解決方案：回答

我試圖解決我的論文在化學微分方程，有我絆了關於SciPy的的微分方程求解「odeint」的問題。SciPy的：實施微分方程的兩種方法：兩種不同的解決方案：回答

首先我由函數CIDNP_1實現差分（CIDNP是一個化學現象，這解釋了不尋常的變量）根據SciPy的網站上的例子。但是，即使正確的方向，解決方案也是如此。

import numpy as np 
import matplotlib.pyplot as plt 
import scipy.integrate 

R0 = 5e+5 
kt = 5e5/R0 
beta = 3/R0 

def CIDNP_1(y, t): 
    dP_dt, dQ_dt = y 

    def R(t): 
     return R0/(1 + kt*R0*t) 

    dP_dt = -kt*dP_dt*R(t) - kt*beta*(R(t))**2 
    dQ_dt = +kt*dP_dt*R(t) + kt*beta*(R(t))**2 
    return [dP_dt, dQ_dt] 


def CIDNP_2(y, t):  
    dP_dt, dQ_dt = y 

    def R(t): 
     return R0/(1 + kt*R0*t) 

    return [-kt*dP_dt*R(t) - kt*beta*(R(t))**2, \ 
      +kt*dP_dt*R(t) + kt*beta*(R(t))**2] 

y0 = [-1, +1] 
t = np.linspace(1e-9, 100e-6, 1e3) 
sol_1 = scipy.integrate.odeint(CIDNP_1, y0, t) 
sol_2 = scipy.integrate.odeint(CIDNP_2, y0, t)

然後，我改變了我的解決方案CIDNP_2什麼都給了正確的結果，但在我眼裏沒有在執行無差異的變量dP_dt和dQ_dt不實施CIDNP_1改變。

那麼爲什麼作爲實施CIDNP_1給出錯誤的結果任何人都可以給我一個提示，我會爲那麼至少最後兩個小時內未完全喪失很幸運。

問候，

雅各布

來源

2017-09-06 JakobJakobson

您確認輸入'y'是語義上的時間導數，而不僅僅是國家回答'y = [P，Q]'RESP。 'P，Q = y'？這種重新貼標籤也會避免發生混淆。 – LutzL

Semantivally你是對的。微分方程描述了系統的狀態並且沒有導數。但是，因爲我最初只是想在五分鐘內編碼，所以我沒有太多考慮。在下一個版本中，我會改變它。 – JakobJakobson

在你的第一個版本，你執行更新不是在同一時間，當你執行到線

dP_dt = -kt*dP_dt*R(t) - kt*beta*(R(t))**2 
dQ_dt = +kt*dP_dt*R(t) + kt*beta*(R(t))**2

不simulatnously;因此，您使用已經更新dP_dt更新dQ_dt。這是ODE系統的錯誤實施。你的第二種方法更好，因爲它沒有這種錯誤。你要麼必須直接返回更新後的值，或者您有計算新dQ_dt之前在另一個變量保存的dP_dt新計算出的值。

new_dP_dt = -kt*dP_dt*R(t) - kt*beta*(R(t))**2 
new_dQ_dt = +kt*dP_dt*R(t) + kt*beta*(R(t))**2 

return [new_dP_dt, new_dQ_dt]

這會解決您的問題。

來源

2017-09-06 20:06:14 FlashTek

非常感謝。下次我會盡量不在同一行中覆蓋我的變量。 Regards – JakobJakobson

在CIDNP_1您使用新的值來計算dQ_dt之前更改dP_dt值：

dP_dt = -kt*dP_dt*R(t) - kt*beta*(R(t))**2 # changed dP_dt! 
dQ_dt = +kt*dP_dt*R(t) + kt*beta*(R(t))**2 # use the changed value!

在CIDNP_2您可以同時計算它們，即dQ_dt用的dP_dt原始值，而不是改變一個計算。你可以認爲它像

a = -kt*dP_dt*R(t) - kt*beta*(R(t))**2  # no overwriting - 
b = +kt*dP_dt*R(t) + kt*beta*(R(t))**2  # uses original value of dP_dt 
return [a, b]

你也很可能加快東西有點像

def CIDNP_3(y, t): 
    dP_dt, dQ_dt = y 
    R_t = R0/(1 + kt * R0 * t) 
    res = kt * R_t * (dP_dt + beta * R_t) 
    return [-res, res]

來源

2017-09-06 20:12:19

非常感謝。下次我會盡量不在同一行中覆蓋我的變量。問候 – JakobJakobson

SciPy的：實施微分方程的兩種方法：兩種不同的解決方案：回答

回答

相關問題