0

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?

0

1 Answer 1

2

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.

Sign up to request clarification or add additional context in comments.

Comments

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.