Zach Thomas - 1 year ago 128
Python Question

# Error in scipy odeint use for system of ODE's

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!

The two equations are as follows:

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

dT/dt = w*(Ti - T)/(Vp) + d_HRkCA/(pC) + UA*(Tc - T)/(VpC)

The code is:

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

Here is what the graphs are supposed to produce:

And here is the output of my code above: