运行欧拉-理查森法代码触发ValueError，请求问题排查

问题场景 #

使用Python结合NumPy实现欧拉-理查森法求解洛伦兹方程组时，执行x[i+1] = x[i] + k2_x等数组赋值操作时，触发ValueError: setting an array element with a sequence错误。此前类似数组分配代码未出现问题，源代码及报错信息如下：

源代码 #

import numpy as np
import matplotlib.pyplot as plt
#parameters
sigma=10.
rho=28.
beta=8/3.
ti=0.
tf=100
dt=0.01
#pre-allocation
x = np.zeros(tf)
y = np.zeros(tf)
z = np.zeros(tf)
#initial conditions
x[0]=1.
y[0]=1.
z[0]=1.
#functions
fx= lambda x: sigma*(y-x)  #y too?
fy= lambda y: x*(rho-z)-y
fz= lambda z: x*y-(beta*z)
#euler-richardson
for i in np.arange(0,tf-1):
    k1_x = fx(x[i])
    k1_y = fy(y[i])
    k1_z = fz(z[i])
    
    k2_x = fx((x[i]+(0.5*k1_x))*dt)  #maybe just dt?
    k2_y = fy((y[i]+(0.5*k1_y))*dt)
    k2_z = fz((z[i]+(0.5*k1_z))*dt)
    
    x[i+1] = x[i] + k2_x
    y[i+1] = y[i] + k2_y
    z[i+1] = z[i] + k2_z

报错信息 #

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
TypeError: only size-1 arrays can be converted to Python scalars

The above exception was the direct cause of the following exception:

ValueError                                Traceback (most recent call last)
Input In [10], in <cell line: 2>()
      8 k2_y = fy((y[i]+(0.5*k1_y))*dt)
      9 k2_z = fz((z[i]+(0.5*k1_z))*dt)
---> 11 x[i+1] = x[i] + k2_x
     12 y[i+1] = y[i] + k2_y
     13 z[i+1] = z[i] + k2_z

ValueError: setting an array element with a sequence.

错误原因分析 #

1. 微分方程函数定义错误：
   洛伦兹方程组的每个方程都依赖x、y、z三个变量，但定义的fx、fy、fz仅接收单个参数，且内部引用了全局的x/y/z数组（而非当前迭代的标量值）。例如fx=lambda x: sigma*(y-x)中的y是整个NumPy数组，当传入标量x[i]时，计算结果是一个数组，而非单个标量，导致后续k1_x、k2_x均为数组，无法赋值给x[i+1]这个标量位置。

2. 欧拉-理查森法公式实现错误：
   欧拉-理查森法中，中间状态点的计算应为x_mid = x[i] + 0.5 * dt * k1_x，而非(x[i]+0.5*k1_x)*dt，颠倒了运算顺序，导致中间状态值完全错误。

3. 数组预分配长度错误：
   tf=100是终止时间，dt=0.01，实际迭代步数应为int((tf - ti)/dt) = 10000，但用np.zeros(tf)只分配了100个元素的数组，后续循环会出现越界问题（当前错误出现在更早步骤，但这是潜在隐患）。

修正后的代码 #

import numpy as np
import matplotlib.pyplot as plt

# 参数定义
sigma = 10.
rho = 28.
beta = 8/3.
ti = 0.
tf = 100
dt = 0.01

# 计算总步数，预分配数组
n_steps = int((tf - ti) / dt)
x = np.zeros(n_steps)
y = np.zeros(n_steps)
z = np.zeros(n_steps)

# 初始条件
x[0] = 1.
y[0] = 1.
z[0] = 1.

# 重构微分方程函数，接收x,y,z三个参数
def fx(x_val, y_val, z_val):
    return sigma * (y_val - x_val)

def fy(x_val, y_val, z_val):
    return x_val * (rho - z_val) - y_val

def fz(x_val, y_val, z_val):
    return x_val * y_val - beta * z_val

# 欧拉-理查森法迭代
for i in range(n_steps - 1):
    # 计算k1：当前点的导数
    k1_x = fx(x[i], y[i], z[i])
    k1_y = fy(x[i], y[i], z[i])
    k1_z = fz(x[i], y[i], z[i])
    
    # 计算中间状态点
    x_mid = x[i] + 0.5 * dt * k1_x
    y_mid = y[i] + 0.5 * dt * k1_y
    z_mid = z[i] + 0.5 * dt * k1_z
    
    # 计算k2：中间点的导数
    k2_x = fx(x_mid, y_mid, z_mid)
    k2_y = fy(x_mid, y_mid, z_mid)
    k2_z = fz(x_mid, y_mid, z_mid)
    
    # 更新下一个状态
    x[i+1] = x[i] + dt * k2_x
    y[i+1] = y[i] + dt * k2_y
    z[i+1] = z[i] + dt * k2_z

# 绘制洛伦兹吸引子
fig = plt.figure()
ax = fig.add_subplot(projection='3d')
ax.plot(x, y, z)
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
plt.show()

修正说明 #

• 重构微分方程函数，明确传入当前迭代的x_val、y_val、z_val标量值，避免引用全局数组。
• 修正欧拉-理查森法的中间状态计算逻辑，符合方法的数学定义。
• 正确计算迭代步数并预分配数组，避免长度不足的问题。
• 补充状态更新时的dt因子，这是原代码遗漏的关键步骤（欧拉类方法需乘以时间步长）。

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