Error in scipy odeint use for system of ODE's
<p>I am trying to reproduce the two graphs shown below, which are temperature and concentration profiles as a function of time. I have checked my method and code a million times but cannot seem to find an error in it but I cannot reproduce those graphs. All values are constants except CA and T. Could it be an issue with the accuracy of odeint from scipy? Any help would be much appreciated! </p>

<p>The two equations are as follows:</p>

<p>dCA/dt = q*(CAi - CA)/V - k*CA</p>

<p>dT/dt = w*(Ti - T)/(V<em>p) + d_HR</em>k<em>CA/(p</em>C) + UA*(Tc - T)/(V<em>p</em>C)</p>

<p>The code is: </p>

<pre><code>import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

def ODESolve(y, t, q, CAi, V, k0, w, Ti, p, dH_R, C, UA, Tc):

    CA, T = y
    k = k0*np.exp(8750*1/T)
    dydt = [q*(CAi - CA)/V - k*CA, w*(Ti - T)/(V*p) + \
            dH_R*k*CA/(p*C) + UA*(Tc - T)/(V*p*C)]

    return dydt

q = 100
CAi = 1.0
V = 100
p = 1000
C = .239
dH_R = 5*(10**4)
k0 = 7.2*(10**10)
UA = 5*10**4
CA0 = .5
T0 = 350
Ti = T0 
w = p*q 

y0 = [CA0, T0]
t = np.linspace(0, 20, 100)

Tc = 305
sol1 = odeint(ODESolve, y0, t, args = (q, CAi, V, k0, w, Ti, p, dH_R, C, UA, Tc))

Tc = 300
sol2 = odeint(ODESolve, y0, t, args = (q, CAi, V, k0, w, Ti, p, dH_R, C, UA, Tc))

Tc = 290
sol3 = odeint(ODESolve, y0, t, args = (q, CAi, V, k0, w, Ti, p, dH_R, C, UA, Tc))

plt.figure(1)
plt.plot(t, sol1[:,0], label = 'Tc = 305')
plt.plot(t, sol2[:,0], label = 'Tc = 300')
plt.plot(t, sol3[:,0], label = 'Tc = 290')
plt.ylim(ymax = 1, ymin = 0)
plt.title ('CA(t)')
plt.legend(loc = 'best')

plt.figure(2)
plt.plot(t, sol1[:,1], label = 'Tc = 305')
plt.plot(t, sol2[:,1], label = 'Tc = 300')
plt.plot(t, sol3[:,1], label = 'Tc = 290')
plt.ylim(ymax = 450, ymin = 300)
plt.legend(loc = 'best')
plt.title ('T(t)')

plt.show()
</code></pre>

<p>Here is what the graphs are supposed to produce: </p>

<p><a href="http://i.stack.imgur.com/XDrLp.png" rel="nofollow"><img src="http://i.stack.imgur.com/XDrLp.png" alt="enter image description here"></a></p>

<p>And here is the output of my code above:</p>

<p><a href="http://i.stack.imgur.com/nq6uS.png" rel="nofollow"><img src="http://i.stack.imgur.com/nq6uS.png" alt="enter image description here"></a>
<a href="http://i.stack.imgur.com/924MS.png" rel="nofollow"><img src="http://i.stack.imgur.com/924MS.png" alt="enter image description here"></a></p>

<p>Apparently there is a sign error in the formula <code>k = k0*np.exp(8750*1/T)</code>.  If you change that to <code>k = k0*np.exp(-8750*1/T)</code>, you'll get plots like those you expect.</p>

python scipy ode odeint