数组与‘非均匀形状’问题求助：ODE积分触发ValueError报错

原本用于三维图形绘制的代码，添加随时间t变化的变量s1后触发如下报错：

ValueError: setting an array element with a sequence. The requested array has an inhomogeneous shape after 1 dimensions. The detected shape was (6,) + inhomogeneous part.

已尝试数组/列表转换、内容检查（确认无虚部），回滚numpy版本失败，只能使用最新版numpy，求其他解决方法。

原代码如下：

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

b = 0.5182
c = 0.22439
d = 1.465
a1 = 30
a2 = 30
eps = 0.1
s = np.random.normal(0, 1, 100000)
n = s[1]
s1 = []
for i in range(1000):
    s1.append(s[i])
as1 = np.asarray (s1)
def f(x):
    if x<=-1.2:
        f = 0.528
    if -1.2<x<1.2:
        f = x*(1-x*x)
    else:
        f = -0.528
    return f
def dH_dt(H, t=0,*s1):
        return np.array([(H[1]+s1),
                     -H[0]-b*H[1]+H[2],
                     c*(a1*f(H[0])-H[2])-d*H[1],

                     (H[4]+eps*(H[0]+H[5]))+s1,
                     -H[3]-b*H[4]+H[5],
                     c*(a2*f(H[3])-H[5])-d*H[4]])

                        
t = np.linspace(0, 100, 1000)
H0 = [0.05, 0.07, 0.09, 0.09, 0.09, 0.09]
H, infodict = integrate.odeint(dH_dt, H0, t, args=(s1,), full_output=True)
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot(H[:,0], H[:,1], H[:,2])
plt.show()

问题根源 #

1. s1传递逻辑错误：odeint会逐个传入t数组中的时间点，但原代码把整个s1列表传给每个时间步，导致H[1]+s1是单个数值加整个列表，生成序列而非单个值，最终返回的数组形状不统一。
2. 参数定义冗余：*s1的可变参数写法无必要，且未实现“随时间取对应s1元素”的逻辑。

修正步骤 #

• 将s1转为numpy数组，方便按索引取值
• 修改dH_dt函数，根据当前t匹配对应的索引，取出该时间步的s1值
• 确保返回的每个元素都是单个数值，而非序列

修正后代码 #

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

b = 0.5182
c = 0.22439
d = 1.465
a1 = 30
a2 = 30
eps = 0.1
s = np.random.normal(0, 1, 100000)
# 直接生成与t点数匹配的numpy数组
s1 = s[:1000].astype(np.float64)

def f(x):
    if x <= -1.2:
        return 0.528
    elif -1.2 < x < 1.2:
        return x * (1 - x*x)
    else:
        return -0.528

def dH_dt(H, t, s1, t_array):
    # 根据当前t找到对应索引，取s1的对应值
    idx = np.argmin(np.abs(t_array - t))
    current_s1 = s1[idx]
    return np.array([
        H[1] + current_s1,
        -H[0] - b*H[1] + H[2],
        c*(a1*f(H[0]) - H[2]) - d*H[1],
        H[4] + eps*(H[0] + H[5]) + current_s1,
        -H[3] - b*H[4] + H[5],
        c*(a2*f(H[3]) - H[5]) - d*H[4]
    ])

t = np.linspace(0, 100, 1000)
H0 = [0.05, 0.07, 0.09, 0.09, 0.09, 0.09]
# 传入t数组和s1数组作为参数
H, infodict = odeint(dH_dt, H0, t, args=(s1, t), full_output=True)

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot(H[:,0], H[:,1], H[:,2])
plt.show()

额外优化说明 #

• 重构f函数分支逻辑，用elif避免多分支同时触发的潜在问题
• 直接生成s1为numpy数组，省去列表转数组的冗余步骤
• 通过np.argmin(np.abs(t_array - t))精准匹配当前时间步对应的s1元素，保证随机值随时间正确变化

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