Having trouble with odeint in python

Question

I am trying to write a script to solve an ODE in python and graph the result. I am using scipy.integrate.odeint for this task. I followed a simple tutorial and modified my code to work with the ODE that I want to solve

import numpy as np
import matplotlib.pyplot as plt
import pylab as plb
from scipy.integrate import ode as odeint
from math import pi

t = np.arange(2, 67., 0.01)

#these are just some values I need for my function
p=2.7 * (10**3)     #kg/m^3   density
V=pi * 0.3042 * ((0.0251 / 2)**2)          #Volume m^3
C=904        #heat capacity J/kgK
A=pi * ((0.0251 / 2)**2)         #Cross sectional area m^2
e= 3490000      #emmisitivy 
b=5.6704 * (10** -8)     #boltzmans
D=0.251         #diameter m
T= 21.2        #ambient temp
# initial condition
x0 = 76

plt.clf()

def f(x, t, p, V, C, A, e, b, D, T):
    y=(1/(p*V*C)) * (A*e*b*np.power(T, 4) - (A*e*b*np.power(x, 4)) - ((1.32*A)/(D**0.25))*(np.power((x-T), 1.25)))
    return y

x = odeint(f, x0, t, (p, V, C, A, e, b, D, T))

# plot the numerical solution
plt.plot(t, x)
plt.xlabel('Time (min)') 
plt.ylabel('Temperature (Celsius)') 
plt.title('Temperature vs Time for a Rough Rod')

# plot empirical data
data = plb.loadtxt('rough.csv', skiprows=2)

x = data[:,0]
y = data[:,1]
sigma = data[:,2]

plt.errorbar(x, y, sigma, linestyle='', fmt='.')

plt.legend(['Numerical Solution', 'Data Points'], loc='best') 

plt.show()

This code runs perfectly for simple function like -x, but the problem is it does not work for the function I am using there. I can get a plot out of this method, but I get the same plot whether I use

y=(1/(p*V*C)) * (A*e*b*np.power(T, 4) - (A*e*b*np.power(x, 4)) - ((1.32*A)/(D**0.25))*(np.power((x-T), 1.25)))

or

y=(1/(p*V*C)) * (A*e*b*pow(T, 4) - (A*e*b*pow(x, 4)) - ((1.32*A)/(D**0.25))*(pow((x-T), 1.25)))

or

y=(1/(p*V*C)) * (A*e*b*np.power(T, 4) - (A*e*b*np.power(x, 4))

I think even this gives me the same plot

y=(1/(p*V*C)) * ((-A*e*b*np.power(x, 4))

also, my value of e should be between 0 and 1 but I have to use HUGE numbers to get something that looks like an exponential decay. I am basically trying to do this: sfu.ca/~rca10/rants/convection.pdf. The ODE in question is at the top of page 4. Where the guy in the link can use normal emissivities, I can't, even though were solving the same ODE's with the same values for everything (I am doing the exact same lab)

What the heck am I doing wrong?

Answer 1

I think that there is an error in your import statement:

from scipy.integrate import ode as odeint

should be:

from scipy.integrate import odeint

The rest should work. Since I do not have your input file ('rough.csv') I cannot reproduce your graph.