描述化學方程式的Scipy ODE錯誤

我有一大組ODE代表生物系統中的化學通量。分子正在反應，隔離和循環。我試圖讓這個功能以這樣一種方式發揮作用，讓我知道在一系列條件下可以生產出多少某些產品。描述化學方程式的Scipy ODE錯誤

我使用這些包

import numpy as np 
from scipy.integrate import odeint 
import matplotlib.pyplot as plt 
import math 
import pylab as p

這是我的職責

def rxn(y,t): 
    k1=1   #population coco * rate of photosynthesis 
    k2=0.5  #population methan * rate of reaction 
    k3=0.01 
    i=25 #day shape approximation 
    Pe=0.1 #approx photosynthetic efficiency 
    Ce=0.1 #approx calcium carbonate production efficiency 
    x =i*math.sin(math.pi*t/24)**2 
    # x= day night shift 

    #Sugar Production 
    r1= x * (y[0]*y[1]) -k1*(y[2]*y[3]) 
    # R1 ** 2cho3 + 2h+ <-> o2+ 2ch3cooh *** 

    # Anaerobic Respiration 
    r2= Pe* -k2*(y[4]*y[5]) 
    # R2 *** ch3cooh  -> co2 + ch4  *** 

    # Calcium carbonate production 
    r3= x* Ce * -k3*(y[6]*y[4]*y[7]) 
    # Ca + 2hco3  -> Caco3 + co2 + h2o 

    dWdt =    + r3 #Water 
    dCdt =  - r2 + r3 #Carbon Dioxide 
    dAdt = r1 + r2   #Acetate 
    dMdt =  - r2   #Methane 
    dOdt = 2*r1     #Oxygen 
    dCadt=    - r3 #Calcium 
    dCbdt= 2*r1   -2*r3 #Carbonate 
    dCacdt=    + r3 #Calcium carbonate 





    return [dWdt,dCdt,dAdt,dMdt,dOdt,dCadt,dCbdt,dCacdt]

然後我的代碼的其餘部分如下

# Timespan (0 - hours - increment) 
tspan=np.arange(0,50,0.1) 

#Starting Concentrations 
#y0 = H2o co2 chcooh ch4 o2 ca hco-3 caco3 
y0=[100,40,10,10,10,80,70,10] 

Conc=odeint(rxn,y0,tspan,full_output = 1,mxstep=5000) 

CW = Conc[:,0] 
CC = Conc[:,1] 
CA = Conc[:,2] 
CM = Conc[:,3] 
CO = Conc[:,4] 
CCa= Conc[:,5] 
CCo= Conc[:,6] 
CCc= Conc[:,7] 

plt.plot(tspan,CC,label='co2') 
plt.plot(tspan,CA,label='ch3cooh') 
plt.plot(tspan,CM,label='ch4') 
plt.plot(tspan,CO,label='o2') 
plt.plot(tspan,CCa,label='Ca') 
plt.plot(tspan,CCo,label='hco3-') 
plt.plot(tspan,CCc,label='caco3') 

plt.xlabel("time (hours)") 
plt.ylabel("moles") 
plt.title("Nutrient Flux?") 
plt.legend() 

p.show()

在運行此我得到一個對融合和類型都有很多錯誤。即;

lsoda-- at t (=r1) and step size h (=r2), the  
     corrector convergence failed repeatedly  
     or with abs(h) = hmin ls 
     in above, r1 = 0.3550854455646D+01 r2 = 0.2492601566412D-09

和

 File "Reaction.py", line 62, in <module> 
    CW = Conc[:,0] 
TypeError: tuple indices must be integers or slices, not tuple

我讀過每個錯誤意味着其他堆棧交換答案，不過，我實在看不出怎麼這麼簡單的一組微分方程的也可以這麼僵而且我也不太明白我從哪裏得到類型錯誤。我對python非常陌生（現在可能很明顯），我有一種感覺，這是由於某種編碼錯誤。我真的很感謝解決這個問題的任何幫助。

來源

2017-03-01 AW1991

檢查跡象。這裏有一些看起來有問題的東西：'r2'和'r3'都是負值，這兩種反應都會產生二氧化碳。那麼爲什麼'dCdt = -r2 + r3'？它不應該是'dCdt = -r2 - r3'嗎？（我不是化學家;這個評論是基於你在代碼中的評論。） –

既然你想要完整的輸出，odeint也會返回一個infodict。

Conc, info_dict = odeint(rxn, y0, tspan, full_output=1, mxstep=5000)

在一個側面說明，也許你想有這對你這樣的問題，建chemreac看看。

來源

2017-03-01 06:44:28 Ohjeah

我剛剛減少了尋求解決方案的時間間隔，發現對於tspan=np.arange(0,5,0.1)解決方案隨着時間增長非常迅速並達到~1e19。對於tspan=np.arange(0,10,0.1)它達到~1e38。所以可能你的解決方案指數趨於無窮大，問題可能出現在你的參數，初始條件或方程中。

爲避免你的問題我只是用Conc=odeint(rxn,y0,tspan,full_output = 1,mxstep=10000)[0]第二個錯誤因爲odeint返回元組(y, infodict)，而我們只需要y。

解決方案tspan=np.arange(0,10,0.1)：在您的微分方程條款

來源

2017-03-01 07:50:39

描述化學方程式的Scipy ODE錯誤

回答

相關問題