求解器在不调用派生回调函数的情况下进行积分
Solver does the integration without calling the derivative callback function
我有一个 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, t, 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], P, solver.y[1:]
states.append(gas.state, t=solver.t)
print("t outside ODE function", solver.t);
print("\n")
求解器具有自适应步长。这意味着它按照适合给定误差容限的内部步骤进行。在从一个步进点到下一个步进点的分段中,解值被插值。因此,时间循环的一系列外部步骤可能会落入同一内部段。如果您将错误容忍度设置为较小的默认水平,情况可能会发生逆转,每个外部值请求都需要几个内部步骤。
我有一个 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, t, 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], P, solver.y[1:]
states.append(gas.state, t=solver.t)
print("t outside ODE function", solver.t);
print("\n")
求解器具有自适应步长。这意味着它按照适合给定误差容限的内部步骤进行。在从一个步进点到下一个步进点的分段中,解值被插值。因此,时间循环的一系列外部步骤可能会落入同一内部段。如果您将错误容忍度设置为较小的默认水平,情况可能会发生逆转,每个外部值请求都需要几个内部步骤。