解决scipy odeint报错：func返回数组须为一维但得ndim=2

报错场景 #

使用Python scipy库的odeint函数求解太阳能集热器吸热管内二苯醚传热流体沿管长温度分布的常微分方程时，触发报错：

The array return by func must be one-dimensional, but got ndim=2
即自定义微分方程右端函数返回的数组必须为一维数组，但实际返回了ndim=2的二维数组。

原始问题代码如下：

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
import math
def scp(y,z):
    dh=  0.07 #internal diameter(m)
    tk = 0.002 #thickness of the absorber tube(m)
    do = dh+ 2*tk #outer diamater of the absorber tube
    kat = 45 #thermal conductivity of the absorber tube 
    usDPE =0.5 #velocity of the diphenyl-ether (m/s)
    I = 800 #global insolation (W/m2K)
    Cr = 25 #concentration ratio
    nuo = 0.785 #optical efficiency
    epsr = 0.095 #emissivity of the solar selcective coating
    sigm =5.67*10**-8 #Stefan Boltzann constant
    Pat= math.pi*dh #perimeter of the tube 
    Tsurr = 293
    Tht=y[0]  
    B1=0.2336 
    B2=120.57
    B3=-147.33
    B4=7E-07 
    CpDPE=(1/(B1+(B2/Tht)+(B3/Tht**2)+(B4*Tht)))*1000 #Specific heat coefficient J/kgK
    B5=1.33718
    B6=-0.9233E-03
    B7=-0.08992E-06 
    rhoDPE=(B5+(B6*Tht)+(B7*(Tht)**2))*1000  #Specific heat coefficient kg/m3
    B8=0.0451
    B9=0.0034
    kDPE=B8*math.exp(B9*Tht) #Thermal conductivity of DPE (W/m.K)
    D=4.832
    mu0=0.063
    T0=136.7
    muDPE=math.exp(math.log(mu0)+((D*T0)/(Tht-T0)))/1000 #viscosity of DPE (kg/m.s)
    ReDPE=(dh*usDPE*rhoDPE)/muDPE #Reynolds number of DPE
    PrDPE=(CpDPE*muDPE)/kDPE #Prandtl number of DPE
    NuDPE=0.023*ReDPE**0.8*PrDPE**0.4 #Nusselt number of DPE
    hj=(NuDPE*kDPE)/dh #heat transfer coefficent of fluid
    Uo= (hj*kat)/(kat+(tk*hj)) #Overall heat transfer coefficent of fluid
    At= ((math.pi)/4)*dh**2 #Area of the tube
    mDPE = rhoDPE*usDPE*At #mass flow rate of DPE
    ppar = [(epsr*sigm),0,0,Uo,-(I*Cr*nuo+(Uo*Tht)+(epsr*sigm*Tsurr**4))]
    Tw = np.roots(ppar)
    Tw1=Tw[3:4]
    Tw2=abs(Tw1)
    print(Tw2)
    dThtdz= (-(4/do)*Uo*(Tht-Tw2))/(rhoDPE*usDPE*CpDPE)
    return [dThtdz]
y0 = 563
# length points
z=np.linspace(0,50,1000)
# solve ODE
y = odeint(scp,y0,z,rtol=None,atol=None)
Tht=y[:,0]
plt.figure(1)
plt.plot(z,Tht)
plt.xlabel('length')
plt.ylabel('Heat transfer fluid temperature (K)')

错误根因 #

• 核心问题出在壁温Tw的取值逻辑：np.roots()返回的四次方程根是一维numpy数组，使用切片Tw[3:4]取到的是形状为(1,)的一维数组而非标量，后续计算得到的dThtdz也是同形状的数组。此时用[dThtdz]包裹返回值，相当于把一个长度为1的数组嵌套进列表，最终返回给odeint的是形状为(1,1)的二维数组，不符合接口要求。
• 冗余问题：微分方程右端函数会在积分的每一步被调用，函数内的print(Tw2)会产生上千行冗余打印，拖慢运行速度。

修正步骤 #

• 将切片取根的代码Tw1=Tw[3:4]改为直接取标量值，避免数组嵌套。
• 调整返回语句，直接返回标量导数值或一维数组，不要做多余的列表嵌套。
• 删除函数内的冗余打印语句，避免无效输出。
• （可选优化）四次方程求解可能返回复数根，增加实根筛选逻辑，确保取到物理合理的正实数壁温。

修正后完整可运行代码 #

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
import math
def scp(y,z):
    dh=  0.07 # 管内径(m)
    tk = 0.002 # 吸热管壁厚(m)
    do = dh+ 2*tk # 吸热管外径(m)
    kat = 45 # 吸热管导热系数 W/(m·K)
    usDPE =0.5 # 二苯醚流速(m/s)
    I = 800 # 太阳辐照度 W/m2
    Cr = 25 # 聚光比
    nuo = 0.785 # 光学效率
    epsr = 0.095 # 选择性涂层发射率
    sigm =5.67*10**-8 # 斯忒藩-玻尔兹曼常数
    Tsurr = 293 # 环境温度 K
    Tht=y[0]  
    # 二苯醚定压比热容 J/(kg·K)
    B1=0.2336 
    B2=120.57
    B3=-147.33
    B4=7E-07 
    CpDPE=(1/(B1+(B2/Tht)+(B3/Tht**2)+(B4*Tht)))*1000
    # 二苯醚密度 kg/m3
    B5=1.33718
    B6=-0.9233E-03
    B7=-0.08992E-06 
    rhoDPE=(B5+(B6*Tht)+(B7*(Tht)**2))*1000
    # 二苯醚导热系数 W/(m·K)
    B8=0.0451
    B9=0.0034
    kDPE=B8*math.exp(B9*Tht)
    # 二苯醚动力粘度 kg/(m·s)
    D=4.832
    mu0=0.063
    T0=136.7
    muDPE=math.exp(math.log(mu0)+((D*T0)/(Tht-T0)))/1000
    # 对流换热系数计算
    ReDPE=(dh*usDPE*rhoDPE)/muDPE
    PrDPE=(CpDPE*muDPE)/kDPE
    NuDPE=0.023*ReDPE**0.8*PrDPE**0.4
    hj=(NuDPE*kDPE)/dh
    Uo= (hj*kat)/(kat+(tk*hj))
    At= ((math.pi)/4)*dh**2
    mDPE = rhoDPE*usDPE*At
    # 求解壁温四次方程，筛选正实根
    ppar = [(epsr*sigm),0,0,Uo,-(I*Cr*nuo+(Uo*Tht)+(epsr*sigm*Tsurr**4))]
    Tw = np.roots(ppar)
    # 取物理合理的、高于流体温度的正实根
    real_roots = Tw[np.isreal(Tw)].real
    Tw2 = real_roots[real_roots > Tht][0]
    # 计算温度沿管长的导数
    dThtdz= (-(4/do)*Uo*(Tht-Tw2))/(rhoDPE*usDPE*CpDPE)
    return dThtdz

y0 = 563 # 入口温度 K
z=np.linspace(0,50,1000) # 管长离散点
y = odeint(scp,y0,z,rtol=1e-6,atol=1e-8)
Tht=y[:,0]
# 绘制温度分布曲线
plt.figure(1)
plt.plot(z,Tht)
plt.xlabel('管长 (m)')
plt.ylabel('传热流体温度 (K)')
plt.show()

注：修正后的代码增加了实根筛选逻辑，避免复数根导致的计算异常，同时调整了求解容差提升计算精度。

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