求解器在不调用导数回调函数的情况下进行积分

问题描述

我有一个 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")
    

解决方法

求解器具有自适应步长。这意味着它在适应给定容错的内部步骤中进行。在从一个步骤点到下一步骤点的部分中,解值被插值。因此,可能会发生时间循环的外部步骤序列落入同一内部段的情况。如果您将容错设置为默认值,则可能会发生相反的情况,即每个外部值请求都需要多个内部步骤。