我正在尝试编写一个脚本来解决ODE in python和绘制结果。我正在使用 scipy.integrate.odeint 来完成这项任务。我遵循了一个简单的教程并修改了我的代码以使用我想要解决的 ODE
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()
这段代码对于像 -x 这样的简单函数运行完美，但问题是它不适用于我在那里使用的函数。我可以从这种方法中得到一个情节，但无论我使用，我都会得到相同的情节
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)))
或者
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)))
或者
y=(1/(p*V*C)) * (A*e*b*np.power(T, 4) - (A*e*b*np.power(x, 4))
我想即使这也给了我同样的情节
y=(1/(p*V*C)) * ((-A*e*b*np.power(x, 4))
此外，我的 e 值应该在 0 和 1 之间，但我必须使用巨大的数字来获得看起来像指数衰减的东西。我基本上是在尝试这样做：sfu.ca/~rca10/rants/convection.pdf。有问题的 ODE 位于第 4 页的顶部。链接中的人可以使用正常的发射率，但我不能，即使正在用相同的值解决相同的 ODE（我正在做完全相同的实验室）
我到底做错了什么？