循环求解刚性ODE系统时遇Could not broadcast input array错误求助

问题根源分析 #

1. 循环变量覆盖错误：原代码中for y in T: y=T完全破坏了循环逻辑，把单个温度值的循环变量强行替换成整个温度数组，导致后续计算的速率常数k1~k5都是数组，最终返回的导数数组形状变为(6,6)，与solve_ivp要求的(6,)一维数组不匹配，触发广播错误。
2. 函数定义位置不当：将ODE右侧函数f放在循环内部，每次循环重复定义函数，不仅冗余，还容易引发变量作用域混淆。
3. 参数传递错误：调用solve_ivp时传递整个温度数组T，而非当前循环的单个温度值，导致函数内部温度参数混乱。
4. 平方根警告：当PB（B组分分压）变为负数时，(KB*PB)**0.5会产生无效值，需添加保护逻辑避免负数开平方。

修正后的代码 #

import scipy as sc
from scipy.integrate import solve_ivp
import numpy as np
import matplotlib.pyplot as plt

# 定义ODE右侧函数，放在循环外部
def f(t, x, current_T):
    FA, FB, FC, FD, FE, FF = x
    e = np.e  # 直接使用numpy的自然常数更准确
    R = 8.314
    Tm = 723.15
    A1 = 5.5
    A2 = 0.686
    A3 = 1.58
    A4 = 2.6
    A5 = 0.787
    E1 = 90500
    E2 = 165000
    E3 = 150000
    E4 = 139000
    E5 = 132000
    KB = 6.54e-12
    KD = 1.19
    m2 = 0.922
    m3 = 0.906
    m4 = 1.23
    m5 = 0.905
    Patm = 0.8 * 101325
    FT = 178.47

    # 计算分压，添加保护避免负数
    PA = max((FA / FT) * Patm, 0)
    PB = max((FB / FT) * Patm, 0)
    PC = max((FC / FT) * Patm, 0)
    PD = max((FD / FT) * Patm, 0)

    # 计算速率常数，使用当前循环的单个温度值
    k1 = e ** (A1 - (E1 / R) * ((1 / current_T) - (1 / Tm)))
    k2 = e ** (A2 - (E2 / R) * ((1 / current_T) - (1 / Tm)))
    k3 = e ** (A3 - (E3 / R) * ((1 / current_T) - (1 / Tm)))
    k4 = e ** (A4 - (E4 / R) * ((1 / current_T) - (1 / Tm)))
    k5 = e ** (A5 - (E5 / R) * ((1 / current_T) - (1 / Tm)))

    # 避免负数开平方
    sqrt_term = np.sqrt(KB * PB) if PB >= 0 else 0
    Tast = 1 / (1 + sqrt_term + KD * PD)
    TB = sqrt_term * Tast
    TD = KD * PD * Tast

    r1 = (k1 / 1000) * TB * PA
    r2 = (k2 / 1000) * (TB ** m2) * PA
    r3 = (k3 / 1000) * (TB ** m3) * PA
    r4 = (k4 / 1000) * (TB ** m4) * PC
    r5 = (k5 / 1000) * (TB ** m5) * PC

    rA = -r1 - r2 - r3
    rB = -r1 - 7 * r2 - 5 * r3 - 6 * r4 - 4 * r5
    rC = r1 - r4 - r5
    rD = r1 + 3 * r2 + 3 * r3 + 2 * r4 + 2 * r5
    rE = 2 * r2 + 2 * r4
    rF = 2 * r3 + 2 * r5

    return [rA, rB, rC, rD, rE, rF]

# 温度数组
T = np.array([250, 300, 350, 400, 450, 500])
x0 = (5, 5, 0, 0, 0, 0)
t0 = 0
t1 = 40

# 遍历每个温度求解
for current_T in T:
    print(f"===== 求解温度 {current_T} K =====")
    soln = solve_ivp(f, (t0, t1), x0, method="Radau", args=(current_T,))
    print("各组分最终流量：")
    print(soln.y[:, -1])  # 打印最终时刻的流量值

关键修改说明 #

• 修正循环逻辑：使用current_T遍历温度数组，确保每次循环处理单个温度值，而非整个数组。
• 函数移至循环外：将ODE函数f定义在循环外部，避免重复定义，同时明确参数current_T用于接收当前循环的温度值。
• 参数传递修正：调用solve_ivp时，args=(current_T,)传递当前循环的单个温度，而非整个数组。
• 添加数值保护：通过max(..., 0)确保分压非负，避免负数开平方引发的警告；同时对平方根项添加条件判断，进一步避免无效值。
• 优化常数使用：将手动定义的e=2.711828替换为np.e，提高计算精度。

内容的提问来源于stack exchange，提问作者Axel Flores