2
我有一个由一个函数描述的方程组。用于ODE系统的Python编码
- 产品从反应物
- 部分产品打破
- 的细分产品的某一百分比被回收到初始反应物
- 系统继续提供更多的产品周期正在进行，直到所有的构建极限反应物在非循环产品内或不可用的“损失产品”内
鉴于产品组成上没有变化。我需要进入系统的反应物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()
问候。
欢迎设置您的大规模行动颂歌堆栈溢出。这里你真正的问题是什么？ – MrLeeh
感谢您的欢迎。我想知道为什么当我试图回收反应物时，会影响它们被利用的比例。 – Uaru427
另外我想知道为什么产品2，迷失1和迷失2不会随着时间积累。我希望系统耗尽反应物，使整个数量陷入这三组中。 – Uaru427