Возникли проблемы с моим кодом химической реакции ODE

Question

Итак, я пишу магистерскую диссертацию, и мне нужно смоделировать реакцию BR (Йод-йодид-осцилляция Бриггса-Раушера), и у меня есть некоторые проблемы с этим кодом. Каждый раз, когда я запускаю (тот же код), я получаю другой сюжет и ошибку. Я пробовал разные точки и разные варианты, но все равно.

lsoda--  at current t (=r1), mxstep (=i1) steps   
       taken on this call before reaching tout     
      in above message,  i1 =       500
******** ODEintWarning: Excess work done on this call (perhaps wrong Dfun type). 
  Run with full_output = 1 to get quantitative information.
  warnings.warn(warning_msg, ODEintWarning)

Кто-нибудь может указать мне правильное направление? Я не знаю, где я ошибся.

Спасибо заранее и извините за все ваши неприятности. Спасибо за ответы на все вопросы.

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

k1 = 3.1*(10**12) #M**-2 s**-1 - k11
k2 = 2.2  #s**-1 - k1-1
k3 = 5.0*10**9 #M**-2 s**-1 k12
k4 = 1.4*10**3 #M**-3 s**-1  k13
k5 = 3.0*10**9 #M**-1 s**-1  k14
k6 = 2.6*10**5 #M**-2 s**-1 k15
k7 = 3.494     #M**-1 s**-1   kc5
k8 = 2.0*10**3 #M**-1 s**-1  kd1

#L(0) = #[I-]
#L(1) = #[I2]
#L(2) = #[HOI]
#L(3) = #[H+]
#L(4) = #[H2O]
#L(5) = #[HOIO]
#L(6) = #[IO3-]
#L(7) = #[O2]
#L(8) = #[CH2(COOH)2]
#L(9) = #[CHI(COOH)2]
#L(10) = #[H2O2]


def BR(L,t):
    H1=k1*L[2]*L[0]*L[3]
    H2=k2*L[1]*L[4]
    H3=k3*L[3]*L[5]*L[0]
    H4=k4*(L[3]**2)*L[6]*L[1]
    H5=k5*L[5]**2
    H6=k6*L[7]*L[6]*L[5]
    H7=k7*L[8]*L[1]
    H8=k8*L[2]*L[10]

    dCdt = -H1 + H2 + (2*H3) + H4 + H5 - H8
    dAdt = -H1 + H2 - H3 - H4 + H7 + H8 
    dDdt = -H1 + H2 - H3 - (2*H4) + H5 - H6 + H7 + H8 
    dBdt = H1-H2-H7
    dEdt = H1-H2+H8 
    dFdt = -H3 + H4 - (2*H5) + H6
    dGdt = -H4 + H5 + H6 
    dHdt = H6  
    dIdt = -H7 
    dJdt = H7
    dKdt = H8 
    return(dCdt,dAdt,dBdt,dDdt,dEdt,dFdt,dGdt,dHdt,dIdt,dJdt,dKdt)
t = np.arange(0,100,1)
L0 = [0,0.02,0.04,0,1,0.04,0,0,10.06,0,0.5]
Conc= odeint(BR,L0,t)
cC= Conc[:,0]
cA= Conc[:,1]
cD= Conc[:,2]
cB= Conc[:,3]
cE= Conc[:,4]
cF= Conc[:,5]
cG= Conc[:,6]
cH= Conc[:,7]
cI= Conc[:,8]
cJ= Conc[:,9]
cK= Conc[:,10]

plt.plot(t,cA)
plt.plot(t,cB)
plt.plot(t,cC)
plt.plot(t,cD)
plt.plot(t,cE)
plt.plot(t,cF)
plt.plot(t,cG)
plt.plot(t,cH)
plt.plot(t,cI)
plt.plot(t,cJ)
plt.plot(t,cK)

Answer 1

Помимо несоответствий, упомянутых в комментариях, вы можете обойти внутренний предел шага, либо уменьшив временной шаг, либо увеличив

mxstep

параметр. В то же время вы должны подумать, достаточны ли для ваших целей допуски по умолчанию около 1e-6 или если вам нужна более высокая точность,

Conc= odeint(BR,L0,t, atol=1e-7, rtol=1e-11, mxstep=5000)  

---

Восстанавливая порядок увеличения букв, назначение реагентов должно быть

A: HOI, B: I-, C: H+, D: I2, E: H2O

и то же самое, что и в вашем коде. Химические уравнения:

A + B + C <--> E + F,     ->: k1 = 3.1e12, <-: k2 = 2.2
B + C + F --> 2A          k3 = 5e9
B + 2C + G --> A + F      k4 = 1.4e3
2F --> A + C + G          k5 = 3e9
C + F + G -- 2F + 0.5H    k6 = 2.6e5
D + I --> B + C + J       k7 = 3.494
A + K --> B + C + E + H   k8 = 2e3

Это должно привести к функции

def BR(L,t):
    A,B,C,D,E,F,G,H,I,J,K = L
    H1 = k1*A*B*C
    H2 = k2*D*E
    H3 = k3*B*C*F
    H4 = k4*B*C**2*F
    H5 = k5*F**2
    H6 = k6*C*F*G
    H7 = k7*D*I
    H8 = k8*A*K 

    dAdt = -H1 + H2 + 2*H3 + H4 + H5 - H8
    dBdt = -H1 + H2 -H3 - H4 + H7 + H8
    dCdt = -H1 + H2 - H3 -2*H4 + H5 - H6 + H7 + H8
    dDdt = H1 - H2 + H7 
    dEdt = H1 - H2 + H8
    dFdt = -H3 + H4 - (2*H5) + H6
    dGdt = -H4 + H5 + H6
    dHdt = 0.5*H6 + H8
    dIdt = -H7 #[CH2(COOH)2] se porablja
    dJdt = H7 #[CHI(COOH)2]nastaja
    dKdt = -H8 #[H2O2] nastaja
    return(dAdt,dBdt,dCdt,dDdt,dEdt,dFdt,dGdt,dHdt,dIdt,dJdt,dKdt)

L0 = [ 8e-11, 1e-10, 0.056, 8e-8, 1, 9e-11, 0.01, 2.5e3, 0.0015, 0, 0.33 ]
t = np.arange(0,2000,1)
Conc= odeint(BR,L0,t)
names = [ '[HOI]', '[I-]', '[H+]', '[I2]', '[H2O]', '[HOIO]', '[IO3-]', '[O2]','[CH2(COOH)2]','[CHI(COOH)2]','[H2O2]']

plt.subplots_adjust(hspace=0.4,wspace=0.55)
for k, name in enumerate(names):
    plt.subplot(6,2,k+1);
    plt.plot(t,Conc[:,k])
    plt.title(name)
plt.show()

(начальные значения от http://www.math.udel.edu/~rossi/Math512/2003/br5.pdf), дающие результаты