问题描述
我有一个 python 代码 (example from Cantera.org),它使用 scipy.integrate.ode 来解决 ODE 系统。代码工作正常,结果是合理的。但是,我注意到有关 ode 求解器的某些内容对我来说没有意义。
我在计算导数向量 (print("t inside ODE function",t)
) 的函数内部 (__call__(self,t,y)
) 和 while 循环 (print("t outside ODE function",solver.t);
) 中的函数外部放置了一个打印函数。
我希望在求解器进行时间积分时必须调用内部打印,然后调用外部打印。换句话说,两个 "t outside ODE function"
不能在没有 "t inside ODE function"
的情况下接连出现。然而,这发生在 while 循环中的一些迭代中,这意味着求解器在不计算导数的情况下进行积分。
我想知道这怎么可能
import cantera as ct
import numpy as np
import scipy.integrate
class ReactorOde:
def __init__(self,gas):
# Parameters of the ODE system and auxiliary data are stored in the
# ReactorOde object.
self.gas = gas
self.P = gas.P
def __call__(self,y):
"""the ODE function,y' = f(t,y) """
# State vector is [T,Y_1,Y_2,... Y_K]
self.gas.set_unnormalized_mass_fractions(y[1:])
self.gas.TP = y[0],self.P
rho = self.gas.density
print("t inside ODE function",t)
wdot = self.gas.net_production_rates
dTdt = - (np.dot(self.gas.partial_molar_enthalpies,wdot) /
(rho * self.gas.cp))
dYdt = wdot * self.gas.molecular_weights / rho
return np.hstack((dTdt,dYdt))
gas = ct.solution('gri30.yaml')
# Initial condition
P = ct.one_atm
gas.TPX = 1001,P,'H2:2,O2:1,N2:4'
y0 = np.hstack((gas.T,gas.Y))
# Set up objects representing the ODE and the solver
ode = ReactorOde(gas)
solver = scipy.integrate.ode(ode)
solver.set_integrator('vode',method='bdf',with_jacobian=True)
solver.set_initial_value(y0,0.0)
# Integrate the equations,keeping T(t) and Y(k,t)
t_end = 1e-3
states = ct.solutionArray(gas,1,extra={'t': [0.0]})
dt = 1e-5
while solver.successful() and solver.t < t_end:
solver.integrate(solver.t + dt)
gas.TPY = solver.y[0],solver.y[1:]
states.append(gas.state,t=solver.t)
print("t outside ODE function",solver.t);
print("\n")
解决方法
求解器具有自适应步长。这意味着它在适应给定容错的内部步骤中进行。在从一个步骤点到下一步骤点的部分中,解值被插值。因此,可能会发生时间循环的外部步骤序列落入同一内部段的情况。如果您将容错设置为默认值,则可能会发生相反的情况,即每个外部值请求都需要多个内部步骤。