我试图用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是不合适的,因为它是刚性微分方程组
文档中没有明确说明,但是
Radau
只是方法的stepper类,成功调用此构造函数的返回是stepper对象。要获得类似于odeint
的结果,请使用通用接口方法kin
具有参数顺序t,x
max_step
必须扮演t_bound
的角色,以限制ODE在未来的评估时间,与您的值兼容的是max_step=350
odeint
的输出,sol.y[j,k]
是组件j
在时间索引k
李>相关问题 更多 >
编程相关推荐