Zach Thomas - 1 year ago 120

Python Question

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)/(V*p) + d_HR*k*CA/(p*C) + UA*(Tc - T)/(V*p*C)

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:

Recommended for you: Get network issues from **WhatsUp Gold**. **Not end users.**

Answer Source

Apparently there is a sign error in the formula `k = k0*np.exp(8750*1/T)`

. If you change that to `k = k0*np.exp(-8750*1/T)`

, you'll get plots like those you expect.

Recommended from our users: **Dynamic Network Monitoring from WhatsUp Gold from IPSwitch**. ** Free Download**