我正在尝试在Python中使用SciPy中的solve_ivp和Randau方法求解一阶ODE(单ODE,而不是系统)。 ODE 非常僵硬,即 ,我正在采用更新的 SciPY 的solve_ivp 方法。对于刚性问题,SciPY 推荐 Randau 积分方法。附上代码片段:
for gamma in gamma_array:
for q_nondim in q_nondim_array:
zs = []
us = []
r = solve_ivp(model, (0,1),[0], method = 'Radau',jac = jacobian,t_eval= [0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8,0.9, 1.0])
if r.status <0 :
print(r.message)
print(r.success)
us = np.array(r.y)
zs = np.array(r.t)
pressure = pressure_generator(us)
pressure_array[i,j]= pressure[-1]
j = j+1
i = i+1
j =0
这里,model是计算du/dz的可调用函数:
def model(t,y):
A1 = (1+t_by_a)*(1+t_by_a)
A2 = A1 + y*y +2*y
A3 = 2*gamma*(np.power((3+1/n)*q_nondim,n))
Dr1 = 2*A1*(1+y)/(A2*A2)
Dr2 = 2*(1+y)/(A2)
Dr3 = 2/(1+y)
Dr4 = 2/((1+y)*(1+y)*(1+y))
A = np.power((1+y),(1+3*n))
dudz = -A3/(A*(Dr1+Dr2-Dr3-Dr4))
return dudz
jacobian 是一个函数,它将 ODE 的雅可比矩阵返回为形状 (1,1) 的数组
def jacobian(t,y):
C = 2*gamma*pow((3+1/n)*q_nondim,n)
a= t_by_a
b=((C*(a**2 + 2*a + (y + 1)**2))*((y + 1)**(1 - 3*n))* ((a**4)* (3*n*(y**2 + 2*y + 2) - 2) + 4*(a**3)*(3*n*(y**2 + 2*y + 2) - 2) + (a**2)* (3*n*(y**4 + 4*y**3 + 13*y**2 + 18*y + 12) - 2*(2*y**4 + 8*y**3 + 15*y**2 + 14*y + 9)) + 2*a*((y + 1)**2)*(3*n*(y**2 + 2*y + 4) - 2*(2*y**2 + 4*y + 5)) + 2*(3*n - 4)*(y + 1)**4))/(2*a*(a + 2)*((a**2)*(y**2 + 2*y + 2) + 2*a*(y**2 + 2*y + 2) + 2*(y + 1)**2)**2)
c = np.array([b]) #convert a list into a 1D array
d = np.reshape(c,(1,1))
return c
压力生成器是一个函数,它获取 ODE 的解并将其转换为所需的变量(在本例中为压力)。
def pressure_generator(u):
A1 = (1+t_by_a)*(1+t_by_a)
A2 = A1 + u*u +2*u
A3 = A1/A2
F1 = A3 + np.log(A3)
A4 = 1/((1+u)*(1+u))
F2 = A4 + np.log(A4)
p = (1/gamma)*(F1 - F2)
return p
但是,这里的集成似乎失败了。 r.success 为 False,r.message 为:
所需步长小于数字之间的间距。
我再次尝试将 t_eval 数组更改为:
t_eval= [0.0,0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1.0]
但即使如此,集成也失败并显示相同的消息。此外,我收到另一条错误消息:
Value error: setting array element with a sequence.
关于:
pressure_array[i,j]= pressure[-1]
代码行。
请帮助我完成此集成并解决 ODE。由于solve_ivp是SciPy的新增功能,因此与其相关的文献相当有限,官方文档也没有多大帮助。
你找到解决办法了吗?我使用 radau 也遇到类似的问题