使用solve_ivp时出现ValueError: need at least one array to concatenate错误求助

问题代码 #

xa_max = (0.7 * k1 + (0.3-gamma) * k2) / (k1 + k2) #(0.5343)
npt = 50
xau = np.linspace(0,xa_max,npt)
V = np.linspace(0.,0.,npt)
Xa = np.zeros([npt,npt])

def batch(xa , t):
    z = 1 / (k1 * (1 - xa) - k2 * (gamma + xa))
    return z

def volumemin(V,xau):
    res = Pb_min - (c0[0] * xau ) / (t_reaz/ V + tc + ts)
    return res

for j in range(npt):
    Xa[j,:] = np.linspace(0, xau[j], npt)
    t0 = np.array([0.])
    xaeval = Xa[j,:] 
    SOL = solve_ivp(batch, (0, Xa[j,npt-1]), t0 , t_eval=xaeval)
    t = SOL.y
    t_reaz = max(t) 
    V0 = 3. #m^3
    V[j] = fsolve(volumemin, V0, args = (xau [j])) #m^3

Vmin = min(V)   

报错信息 #

File "C:\Users\Lenovo\Desktop\ex2_esame.py", line 36, in <module>
    SOL = solve_ivp(batch, (0, Xa[j,npt-1]), t0 , t_eval=xaeval)

  File "C:\Users\Lenovo\anaconda3\lib\site-packages\scipy\integrate\_ivp\ivp.py", line 650, in solve_ivp
    ts = np.hstack(ts)

  File "<__array_function__ internals>", line 5, in hstack

  File "C:\Users\Lenovo\anaconda3\lib\site-packages\numpy\core\shape_base.py", line 346, in hstack
    return _nx.concatenate(arrs, 1)

  File "<__array_function__ internals>", line 5, in concatenate

ValueError: need at least one array to concatenate

错误原因 #

1. 微分方程函数参数顺序错误：solve_ivp要求传入的微分方程函数第一个参数是自变量t，第二个是状态变量y，但你的batch函数把xa（状态变量）放在了第一个参数位置，导致积分逻辑混乱。
2. 无效积分区间：当j=0时，xau[0] = 0，所以Xa[j,npt-1] = 0，积分区间变成(0,0)，此时solve_ivp无法生成有效的积分点，导致内部执行np.hstack时没有可拼接的数组，触发报错。

修复方案 #

1. 修正batch函数的参数顺序，将t放在第一个位置：

   def batch(t, xa):
       z = 1 / (k1 * (1 - xa) - k2 * (gamma + xa))
       return z
   

2. 跳过无效的积分区间（j=0的情况），或者给xau的起始值设置一个极小的非零值，避免积分区间长度为0：

   # 方案1：跳过j=0
   for j in range(1, npt):
       # 原循环内的代码
   # 方案2：修改xau的起始值
   xau = np.linspace(1e-6, xa_max, npt)
   

3. 增加积分有效性判断，避免max(t)报错：

   if SOL.success:
       t_reaz = np.max(t)
   else:
       # 处理积分失败的情况，比如设默认值或跳过
       continue
   

修改后的完整代码示例 #

import numpy as np
from scipy.integrate import solve_ivp
from scipy.optimize import fsolve

# 假设这些变量已提前定义
k1 = 1.0
k2 = 0.5
gamma = 0.1
Pb_min = 0.8
c0 = [1.0]
tc = 0.2
ts = 0.1

xa_max = (0.7 * k1 + (0.3-gamma) * k2) / (k1 + k2) #(0.5343)
npt = 50
xau = np.linspace(1e-6, xa_max, npt)  # 用极小值替代0，避免无效区间
V = np.linspace(0.,0.,npt)
Xa = np.zeros([npt,npt])

def batch(t, xa):
    z = 1 / (k1 * (1 - xa) - k2 * (gamma + xa))
    return z

def volumemin(V_val, xau_val):
    res = Pb_min - (c0[0] * xau_val ) / (t_reaz/ V_val + tc + ts)
    return res

for j in range(npt):
    Xa[j,:] = np.linspace(0, xau[j], npt)
    t0 = np.array([0.])
    xaeval = Xa[j,:] 
    # 确保积分区间有效
    if Xa[j,npt-1] <= 0:
        continue
    SOL = solve_ivp(batch, (0, Xa[j,npt-1]), t0 , t_eval=xaeval)
    if not SOL.success:
        continue
    t = SOL.y
    t_reaz = np.max(t) 
    V0 = 3. #m^3
    V[j] = fsolve(volumemin, V0, args = (xau[j])) #m^3

Vmin = min(V[V != 0])  # 排除未赋值的0值

内容的提问来源于stack exchange，提问作者Lorenzo Giannetti