用於ODE系統的Python編碼

Question

我有一個由一個函數描述的方程組。

1. 產品從反應物
2. 部分產品打破
3. 的細分產品的某一百分比被回收到初始反應物
4. 系統繼續提供更多的產品週期正在進行，直到所有的構建極限反應物在非循環產品內或不可用的「損失產品」內

鑑於產品組成上沒有變化。我需要進入系統的反應物1的量與進入系統的反應物2的量成正比。因此，當全部反應物1被消耗時，不再消耗反應物2。

目前反應物消耗的比例在沒有反應物再循環時是恆定的，但是當反應物在react1=-R1 - R5和react2=-R2 - R6的循環中循環時，所用反應物的比例改變。

第二個問題是，在循環產品2期間，丟失的產品不會持續增加。相反，他們似乎與產品1和回收產品分別保持固定比例。

我認爲這兩個問題都是由於我試圖回收系統內的反應物而引起的。任何幫助，將不勝感激。

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

def sample_func (y,t): 


    k1 = 10**-4 
    k2 = k1/4 
    k3 = 0.1 
    Recycle0=0.8 
    Recycle2=0.7 


    R1= -k1*y[0]*y[2] # Rate of substance 1 consumption 
    R2= -k2*y[0]*y[2] # Rate of substance 2 consumption 
    #These must be constantly proportional to one another 

    R3= 0.2*R1+0.7*R2 #Product 1 
    R4= 0.8*R1+0.3*R2 #Product 2 

    R5=R3*Recycle0  #Recycled substance 1 of product 1 
    R6=R3*Recycle2  #Recycled substance 2 of product 1 

    R7=R3*(1-Recycle0) 
    R8=R3*(1-Recycle2) 

    used1 =  R1 
    react1=  -R1 - R5 
    used2 =  R2 
    react2=  -R2 - R6 
    prod1=  -R3 
    prod2=  -R4 
    recycledr1 =-R5 
    recycledr2 =-R6 
    lost1  =-R7 
    lost2  =-R8 

    return [used1, react1, used2,react2,prod1,prod2, recycledr1,recycledr2,lost1,lost2] 





y0=(3,11,3,12,0.01,0.01,0.01,0.01,0.01,0.01) 
tspan=np.arange(0,15000,1) 
Conc= odeint(sample_func,y0,tspan) 



react1   = Conc[:,0] 
used1   = Conc[:,1] 
react2   = Conc[:,2] 
used2   = Conc[:,3] 
prod1   = Conc[:,4] 
prod2   = Conc[:,5] 
recycledr1  = Conc[:,6] 
recycledr2  = Conc[:,7] 

print("Consumed R1 & R2 RATIOS AT DIFFERENT TIME POINTS") 
print((Conc[1:2,1]-Conc[0:1,1])/(Conc[1:2,3]-Conc[0:1,3]), " 1 HOURS") 
print((Conc[50:51,1]-Conc[0:1,1])/(Conc[50:51,3]-Conc[0:1,3]), "50 HOURS") 
print((Conc[1000:1001,1]-Conc[0:1,1])/(Conc[1000:1001,3]-Conc[0:1,3]), "1000 HOURS") 

plt.plot(tspan,react1,label='react1') 
plt.plot(tspan,used1,label='used1') 
plt.plot(tspan,react2,label='react2') 
plt.plot(tspan,used2,label='used2') 
plt.plot(tspan,prod1,label='product1') 
plt.plot(tspan,prod2,label='product2') 
plt.plot(tspan,recycledr1,label='recycled react 1') 
plt.plot(tspan,recycledr2,label='recycled react 2') 


plt.xlabel("time (hours)") 
plt.ylabel("quantity") 
plt.title("production v time") 
plt.legend() 

p.show() 

問候。

Answer 1

OP Duffymo提供了一個合理的觀察，在我看來。每個反應物應該有一個單獨的ODE。我自己處理了這種性質的問題幾次，我也有幾條建議：

- 如果來自質量作用法則的每個ODE都是線性的（即沒有濃度相乘），隱式方法很容易實現。我建議先從「梯形方法」開始，然後再實施高階整合方案，讓人噁心。

- 如果你的系統是非線性的，這個問題有點難以解決。如果所有的均衡/速率常數的幅度相似，那麼可以使用明確的整合方案，而不會有太多麻煩。顯式RK4不應太難實現。

- 我會建議您驗證您在科學文獻中發佈的針對大規模行動派生的任何微分方程組，尤其是如果化學動力學不是您的專業領域。如果這是一個通常被調查的化學系統，這樣的推導不應該太難找到。

- 如果您的質量操作系統已正確設置並且您的集成方案正常工作，應該明顯地保持質量守恆。爲了證實這一點，我將選擇這個系統中所有形式（即二氧化碳和一氧化碳中的碳）具有相同化學計量係數的元素，並確認它在系統中的丰度是恆定的。

這裏是一個鏈接，可以幫助您開始

Law of Mass Action