python - 描述化学方程的 Scipy ODE 错误

Question

我有一大组代表生物系统中化学通量的常微分方程。分子发生 react 、隔离和循环的地方。我试图让它发挥作用，让我了解在一系列条件下可以生产多少特定产品。

我正在使用这些软件包

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

我已经阅读了其他堆栈交换答案中每个错误的含义，但是，我真的不明白这样一组简单的 ODE 为何会如此僵硬，而且我真的不明白我在哪里得到类型错误来自任一。我对 python 非常陌生(现在可能很明显)，我有一种感觉这是由于某种编码错误造成的。我真的很感激任何帮助解决这个问题的人。

Answer 1

我只是减少了寻求解决方案的时间间隔，发现对于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)的解决方案: