Chork-Niem 2.2 Kin-Kat

cn_2.2.py

from scipy import integrate
from numpy import array, linspace
from matplotlib import pyplot as plt
%matplotlib inline
xi = -0.1+1/3
k1=0.1
k_ggw=2**2*xi**3/((1/3-xi)*(2/3-2*xi)**2)
k2=k1/k_ggw
print('xi={:f}'.format(xi))
print('K={:f}'.format(k_ggw))
print('k2={:f}'.format(k2))
def df_dt(y, t):
    [ca,cb,cc,cd]=y
    output = 'r1={:f} r2={:f}'
    output += ' Q={:f}'
    output += ' ca/sum(c)={:f}'
    output = output.format(
        k1*ca*cb**2, 
        k2*cc**2*cd,
        cc**2*cd/ca/cb**2,
        ca/sum([ca,cb,cc,cd])
    )
    #print(output)
    dca=-k1*ca*cb**2+k2*cc**2*cd
    dcb=2*(-k1*ca*cb**2+k2*cc**2*cd)
    dcc=-2*(-k1*ca*cb**2+k2*cc**2*cd)
    dcd=-1*(-k1*ca*cb**2+k2*cc**2*cd)
    return [dca, dcb, dcc, dcd]
t=linspace(0,500,50)
c0=array([0.5,1,0,0])
soln = integrate.odeint(df_dt, c0,t)
plt.plot(t, soln)
plt.legend( labels=['a','b','c','d']);
c=soln[-1,:]
print('c='+str(c))
print('y='+str(c/sum(c)))