2017-03-01 14 views
1

私は、生体系における化学フラックスを表す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(おそらく今明らかに新しい)と私はこれが何らかの種類のコーディングエラーのために感じている。私は本当にこれを解決する助けに感謝します。

+0

あなたの微分方程式における用語の兆候を確認してください。ここに私にとって疑問に思われるものがあります:「r2」と「r3」は両方とも負であり、これらの反応は両方とも二酸化炭素を生成します。ではなぜ 'dCdt = -r2 + r3'ですか?それは 'dCdt = -r2 - r3'であってはなりませんか? (私は化学者ではなく、このコメントはコード内のコメントに基づいています。) –

答えて

0

完全な出力が必要なので、odeintも同様にインポータードを返します。サイドノートでは

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

、多分あなたは問題のあなたの種類のために構築されていchemreacの表情を持っていると思います。

1

解決策を探す時間間隔を短くしたところ、tspan=np.arange(0,5,0.1)の溶液は時間とともに非常に急速に成長し、~1e19に達することがわかりました。 tspan=np.arange(0,10,0.1)については、~1e38に達します。だから、あなたの解は指数関数的に無限大になり、問題はあなたのパラメータ、初期条件、または方程式にあるかもしれません。我々はyのみ必要としながら、odeint戻り、(y, infodict)をタプルので、私はちょうどConc=odeint(rxn,y0,tspan,full_output = 1,mxstep=10000)[0] を使用し、あなたの問題の2番目のエラーを回避するために

tspan=np.arange(0,10,0.1)ため

ソリューション: enter image description here

関連する問題