scipy.odeint求解时np.array无法与y值相乘报错是什么原因？

• 维度不匹配：scipy.integrate.odeint要求传入的初始状态y0为一维numpy数组，且微分方程函数返回的导数维度必须与y0完全一致。你代码中L、S1~S6都是长度为17的数组，计算得到的dhdt1、dhdt2也都是长度17的数组，二者组合后总维度为34，但你定义的y0为单值标量或单元素列表，维度完全不匹配，触发运算报错。
• 状态格式不符合要求：odeint不支持二维状态输入，所有状态必须拼接为一维向量传入，计算导数后也需返回同样长度的一维向量，不能返回嵌套列表。

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# 提前定义你的常量，此处为占位提示，请替换为实际数值
# Th231_LambdaAux = 长度为17的数组
# phi_1~phi_6为你的常量值
# Th231_S_1~Th231_S_6为长度为17的数组

def flux(y,t):
    # 将一维状态向量拆分为两个长度17的数组，对应两个状态量
    y1, y2 = y.reshape(2, 17)
    # 常量赋值
    L = Th231_LambdaAux
    k1 = phi_1
    k2 = phi_2
    k3 = phi_3
    k4 = phi_4
    k5 = phi_5
    k6 = phi_6
    S1 = Th231_S_1
    S2 = Th231_S_2
    S3 = Th231_S_3
    S4 = Th231_S_4
    S5 = Th231_S_5
    S6 = Th231_S_6 

    # 通量计算
    fout1 = (1-0.9998)*y1
    fout2 = 0.9998*y2

    # 公共系数项
    coef = L + k1*S1 + k2*S2 + k3*S3 + k4*S4 + k5*S5 + k6*S6
    dhdt1 = coef * fout1
    dhdt2 = coef * fout2
    # 合并为长度34的一维数组返回，和y0维度匹配
    return np.hstack([dhdt1, dhdt2])

# 时间序列
t = np.linspace(0,315360000) 

# 初始条件：17个点的y1初始值 + 17个点的y2初始值，拼接为长度34的一维数组
# 示例为所有点的y1、y2初始值都为1e22，可根据需求自行修改不同点的初始值
y0 = np.ones(2*17) * 1e22

# 积分求解
y = odeint(flux,y0,t)

# 绘图示例：此处以第0个网格点的y1、y2为例，可修改索引查看其他点的结果
plt.figure(1)
plt.loglog(t,y[:,0],'b-', label='y1')
plt.loglog(t,y[:,17],'r--', label='y2')
plt.xlabel('Time (s)')
plt.ylabel('N [nuclide/cm3]')
plt.legend()
plt.show()

如果两个状态量的物理逻辑和示例拆分方式不同，可以自行调整reshape和拼接的规则，只要保证初始值维度和导数返回维度一致即可。

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