Python codant pour le système ODE

Question

J'ai un système d'équations décrites par une fonction dans laquelle.

1. produits sont construits à partir de réactifs
2. certains produits se décomposent
3. un certain pourcentage des produits ventilés sont recyclés dans réactifs initiaux
4. le système continue de cycles avec d'autres produits étant fabriqués jusqu'à ce que toutes le réactif limitant fait partie du produit sans recyclage ou des "produits perdus" inutilisables

Étant donné que les produits ne changent pas de composition. J'ai besoin que la quantité de réactif 1 entrant dans le système soit directement proportionnelle à la quantité de réactif 2 entrant dans le système. Ainsi, lorsque tout le réactif 1 est consommé, pas plus du réactif 2 n'est consommé. Actuellement, le rapport de la consommation de réactifs est constant lorsqu'il n'y a pas de recyclage des réactifs, cependant, lorsque les réactifs sont recyclés dans les lignes react1=-R1 - R5 et react2=-R2 - R6, le rapport des réactifs utilisés est modifié.

Le deuxième problème est que, pendant le cyclage du produit 2, les produits perdus n'augmentent pas continuellement. Au lieu de cela, ils semblent maintenir un ratio fixe avec le produit 1 et les produits recyclés respectivement.

Je suppose que les deux problèmes sont causés par la façon dont j'ai essayé de recycler les réactifs dans le système. Toute aide serait appréciée.

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() 

Cordialement.

Answer 1

OP Duffymo offre une observation sonore, à mon avis. Vous devriez avoir une EDO individuelle pour chaque réactif. Ayant moi-même traité quelques problèmes de cette nature, j'ai aussi quelques conseils:

-Si chaque ODE des lois d'action de masse est linéaire (c'est-à-dire qu'aucune concentration n'est multipliée), une approche implicite est facilement mise en œuvre. Je suggérerais de commencer par la «méthode trapézoïdale», puis de mettre en œuvre des schémas d'intégration d'ordre supérieur et des nausées.

-Si votre système n'est pas linéaire, ce problème est un peu plus difficile à résoudre. Si les magnitudes de toutes vos constantes d'équilibre/vitesse sont similaires à celles que vous pouvez probablement utiliser un schéma d'intégration explicite sans trop de problèmes. Explicit RK4 ne devrait pas être trop difficile à mettre en œuvre. Je suggérerais de vérifier tout système d'équations différentielles que vous développez contre les dérivations d'actions de masse publiées dans la littérature scientifique, surtout si la cinétique chimique n'est pas votre domaine d'expertise. S'il s'agit d'un système chimique couramment étudié, une telle dérivation ne devrait pas être trop difficile à trouver. -Si votre système d'action de masse a été configuré correctement et que votre schéma d'intégration fonctionne correctement, la conservation de la masse doit être maintenue de manière visible. Pour vérifier cela, je sélectionnerais un élément de ce système qui a des coefficients stœchiométriques égaux dans toutes les formes (c'est-à-dire le carbone dans le CO2 et le CO) et vérifie que son abondance dans le système est constante.

Voici un lien qui pourrait vous aider à démarrer la mise en place de vos EDO d'action de masse

Law of Mass Action