用Radau方法求解scipy.integrated在python的中子动力学方程

Question

在python中，我尝试用RADAU方法求解两种反馈(燃料温度反馈和冷却剂温度反馈)的中子动力学方程组。

import numpy as np
from scipy.integrate import Radau

def kin(x, t):
    beta = []
    lam = []
    lam = [0.001334, 0.032739, 0.12078, 0.30278, 0.84949, 2.853]
    beta = [0.000256, 0.00146, 0.001306, 0.002843, 0.000937, 0.000202]
    lifetime = 0.000015
    betasum = sum(beta)
    alfa_ttop = -0.000018
    alfa_ttn = -0.00026
    po0 = -1.0 * betasum
    n0 = 35.2 * 1000000
    Ttop0 = 377
    mtop = 1469.71
    ctop = 300
    kt = 11000
    Tvh = 271
    Gtn = 179.9
    ctn = 5500
    gamv = 900
    mtn = 500
    n = x[0]
    c1 = x[1]
    c2 = x[2]
    c3 = x[3]
    c4 = x[4]
    c5 = x[5]
    c6 = x[6]
    Ttop = x[7]
    Ttn = x[8]
    dndt = (po0 + alfa_ttop * (Ttop - Ttop0) + alfa_ttn * (Ttn - Tvh) - betasum) / lifetime * n + lam[0] * c1 + lam[1] * c2 + lam[2] * c3 + lam[3] * c4 + lam[4] * c5 + lam[5] * c6
    dc1dt = beta[0] / lifetime * n - lam[0] * c1
    dc2dt = beta[1] / lifetime * n - lam[1] * c2
    dc3dt = beta[2] / lifetime * n - lam[2] * c3
    dc4dt = beta[3] / lifetime * n - lam[3] * c4
    dc5dt = beta[4] / lifetime * n - lam[4] * c5
    dc6dt = beta[5] / lifetime * n - lam[5] * c6
    dTtopdt = 1.0 / (mtop * ctop) * (n - kt * (Ttop - Ttn))
    dTtndt = 1.0 / (mtn * ctn) * (kt * (Ttop - Ttn) - gamv * ctn * Gtn * (Ttn - Tvh))

   return (dndt, dc1dt, dc2dt, dc3dt, dc4dt, dc5dt, dc6dt, dTtopdt, dTtndt)


n0 = 35.2 * 1000000
beta = []
lam = []
lam = [0.001334, 0.032739, 0.12078, 0.30278, 0.84949, 2.853]
beta = [0.000256, 0.00146, 0.001306, 0.002843, 0.000937, 0.000202]
lifetime = 0.000015
Tvh = 271
Ttop0 = 377

x0 = np.array(
[n0, beta[0] * n0 / (lifetime * lam[0]), beta[1] * n0 / (lifetime * lam[1]), beta[2] * n0 / (lifetime * lam[2]),
 beta[3] * n0 / (lifetime * lam[3]), beta[4] * n0 / (lifetime * lam[4]), beta[5] * n0 / (lifetime * lam[5]), Ttop0,
 Tvh])
 t = np.linspace(0, 350, 700)
 t_bound = 700
 x = Radau(kin, t, x0, t_bound)

 n = x[:, 0]
 for i in range(0, len(n)):
    print(t[i], n[i] / 1000000)

我收到了下一个错误：

Traceback (most recent call last):
File "D:/Apps/untitled2/Scripts/RADAU.py", line 62, in <module>
    x = Radau(kin, t, x0, t_bound)
  File "D:\Apps\lib\site-packages\scipy\integrate\_ivp\radau.py", line 288, in __init__
    super(Radau, self).__init__(fun, t0, y0, t_bound, vectorized)
  File "D:\Apps\lib\site-packages\scipy\integrate\_ivp\base.py", line 145, in __init__
    self.direction = np.sign(t_bound - t0) if t_bound != t0 else 1
ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

我该怎么做才能纠正这个错误呢？首先，我在python中用odeint进行了求解，它成功了，但是我发现对于这个方程组，odeint是不合适的，因为它是刚性微分方程组。

Answer 1

在文档中没有明确说明这一点，但是Radau只是该方法的stepper类，成功调用这个构造函数的返回是一个stepper对象。要获得类似于odeint的结果，请使用通用接口方法

sol = solve_ivp(kin, (t[0],t[-1]), x0, t_eval=t, max_step=max_step, method='Radau', atol=atol, rtol=rtol)
x = sol.y

• kin有参数顺序t,x，
• max_step必须扮演t_bound的角色来限制ODE的计算范围，兼容您的值将是max_step=350，
• 建议始终明确地控制公差，特别是为了适应期望的状态向量的规模，
• 输出也转换到odeint的输出，sol.y[j,k]是组件j在时间索引k.

>