Python中ODE求解器的错误

Question

我正在使用一个Spark内燃机模型，由于某些原因，我正在使用python对燃烧进行建模。我试着使用颂歌的解算器，但结果完全不符合实际。我发现圆柱体体积的积分是错误的。我已经尝试使用“odeint”和“ode”解算器，但结果是一样的。代码显示体积与θ的导数，并进行积分以找到体积。我把分析方程式拿来比较。

OBS：我在使用Matlab时也遇到过类似的问题，但当我尝试在三角函数中使用度时就遇到了。当我改变弧度时，问题就解决了。守则如下：

from scipy.integrate import odeint
from scipy.integrate import ode
from scipy import integrate
import math
import sympy
from sympy import sqrt, sin, cos, tan, atan
from pylab import *
from RatesComp import *
V_real=np.zeros((100))

def Volume(V,theta):
   V_sol = V[0]
   dVdtheta = Vtdc*(r-1)/2 *( sin(theta) + eps/2*sin(2*theta)/sqrt(1-(eps**2)*sin(theta)**2))
   return [dVdtheta]

#Geometry
eps = 0.25;        #  half stroke to rod ratio, s/2l
r = 10;            #  compression ratio
Vtdc = 6.9813e-05  # volume at TDC

# Initial Conditions
theta0 = - pi
V_init =  0.0006283
theta = linspace(-pi,pi,100)
solve = odeint( Volume, V_init, theta)

# Analytical Result
Size = len(theta)

for i in range(0, Size,1):
    V_real[i] = Vtdc*(1+(r-1)/2*(1-cos(theta[i])+ 1/eps*(1-(1-(eps**2)*sin(theta[i])**2)**0.5)))

figure(1)
plot(theta, solve[:,0],label="Comput")
plot(theta, V_real[0:Size],label="Real")
ylabel('Volume [m^3]')
xlabel('CA [Rad]')
legend()
grid(True)
show()

我展示的图是圆柱体的体积。计算结果与实际相符

有人能帮助了解为什么会发生这个问题吗？

Answer 1

如果使用python2，请在第一行添加：from\uuuuuu future\uuuuu import division以根据文档使用Python3中的division：https://mail.python.org/pipermail/tutor/2008-March/060886.html

在python2中，当您除以两个整数值时，将得到一个整数结果，而不是浮点值。它可以解决您的问题，而无需对代码进行大量更改。

Answer 2

显然你用的是python2在这里，r=10的声明给出了r类型整数，这导致了在'real'解决方案中的（r-1）/2中不需要的整数除法。在导数函数中，有一个浮点值Vtdc作为产品中的第一个因子，之后整个产品的评估都是浮点的。

因此，更改为r=10.0或使用（r-1.0）/2或0.5*（r-1）。

您还应该设置V_init=r*Vtdc，因为这是V_real（-pi）的值。