Sympy与Numpy结合特殊函数出错：DiracDelta未定义

问题代码 #

我尝试创建螺线管磁场数值解的可视化图形，编写了如下代码：

import numpy as np
import sympy as smp
from scipy.integrate import quad_vec

t, x, y, z = smp.symbols("t x y z")

l = smp.Matrix([0.5 * smp.cos(t), 
                0.5 * smp.sin(t), 
                (t / (600 * smp.pi)) * smp.Heaviside(t - 1) * smp.Heaviside(t)])

r = smp.Matrix([x, y, z])
sep = r-l

integrand = smp.diff(l, t).cross(sep) / sep.norm()**3

dBxdt = smp.lambdify([t, x, y, z], integrand[0])
dBydt = smp.lambdify([t, x, y, z], integrand[1])
dBzdt = smp.lambdify([t, x, y, z], integrand[2])

def B(x, y, z):
    return np.array([quad_vec(dBxdt, 0, 2*np.pi, args=(x, y, z))[0],
                     quad_vec(dBydt, 0, 2*np.pi, args=(x, y, z))[0],
                     quad_vec(dBzdt, 0, 2*np.pi, args=(x, y, z))[0]])

x = np.linspace(-2, 2, 20)
xv, yv, zv = np.meshgrid(x, x, x)

B_field = B(xv, yv, zv)
Bx, By, Bz = B_field

错误信息 #

运行后出现如下错误：

File <lambdifygenerated-8>:2, in _lambdifygenerated(t, x, y, z)
      1 def _lambdifygenerated(t, x, y, z):
----> 2     return (-(y - 0.5*sin(t))*((1/600)*t*DiracDelta(t)*select([less(t, 1),equal(t, 1),True], [0,1/2,1], default=nan)/pi + (1/600)*t*DiracDelta(t - 1)*select([less(t, 0),equal(t, 0),True], [0,1/2,1], default=nan)/pi + (1/600)*select([less(t, 0),equal(t, 0),True], [0,1/2,1], default=nan)*select([less(t, 1),equal(t, 1),True], [0,1/2,1], default=nan)/pi) + 0.5*(-1/600*t*select([less(t, 0),equal(t, 0),True], [0,1/2,1], default=nan)*select([less(t, 1),equal(t, 1),True], [0,1/2,1], default=nan)/pi + z)*cos(t))/(abs(x - 0.5*cos(t))**2 + abs(y - 0.5*sin(t))**2 + abs((1/600)*t*select([less(t, 0),equal(t, 0),True], [0,1/2,1], default=nan)*select([less(t, 1),equal(t, 1),True], [0,1/2,1], default=nan)/pi - z)**2)**(3/2)

NameError: name 'DiracDelta' is not defined

解决方案 #

错误原因 #

对Heaviside阶跃函数求导后会生成DiracDelta狄拉克δ函数，但lambdify默认不会导入该函数，且数值积分场景下不需要符号化的δ函数——原代码用两个Heaviside函数限制t的范围，完全可以替换为分段函数避免引入δ函数。

修改步骤 #

1. 替换Heaviside为分段表达式
   用smp.Piecewise替代Heaviside组合，明确定义t在不同区间的取值，求导后不会产生δ函数：

   l = smp.Matrix([
       0.5 * smp.cos(t), 
       0.5 * smp.sin(t), 
       smp.Piecewise(
           (0, t < 0),
           (t / (600 * smp.pi), (t >= 0) & (t <= 1)),
           (0, t > 1)
       )
   ])
   

2. 指定lambdify的映射模块
   调用lambdify时添加modules参数，确保符号函数正确映射到numpy/scipy的数值实现：

   dBxdt = smp.lambdify([t, x, y, z], integrand[0], modules=['numpy', 'scipy'])
   dBydt = smp.lambdify([t, x, y, z], integrand[1], modules=['numpy', 'scipy'])
   dBzdt = smp.lambdify([t, x, y, z], integrand[2], modules=['numpy', 'scipy'])
   

完整可运行代码 #

import numpy as np
import sympy as smp
from scipy.integrate import quad_vec

t, x, y, z = smp.symbols("t x y z")

# 用Piecewise替代Heaviside，避免求导产生DiracDelta
l = smp.Matrix([
    0.5 * smp.cos(t), 
    0.5 * smp.sin(t), 
    smp.Piecewise(
        (0, t < 0),
        (t / (600 * smp.pi), (t >= 0) & (t <= 1)),
        (0, t > 1)
    )
])

r = smp.Matrix([x, y, z])
sep = r - l

integrand = smp.diff(l, t).cross(sep) / sep.norm()**3

# 指定modules确保函数正确映射到数值实现
dBxdt = smp.lambdify([t, x, y, z], integrand[0], modules=['numpy', 'scipy'])
dBydt = smp.lambdify([t, x, y, z], integrand[1], modules=['numpy', 'scipy'])
dBzdt = smp.lambdify([t, x, y, z], integrand[2], modules=['numpy', 'scipy'])

def B(x, y, z):
    return np.array([
        quad_vec(dBxdt, 0, 2*np.pi, args=(x, y, z))[0],
        quad_vec(dBydt, 0, 2*np.pi, args=(x, y, z))[0],
        quad_vec(dBzdt, 0, 2*np.pi, args=(x, y, z))[0]
    ])

# 生成网格并计算磁场
x = np.linspace(-2, 2, 20)
xv, yv, zv = np.meshgrid(x, x, x)

B_field = B(xv, yv, zv)
Bx, By, Bz = B_field

可选：添加可视化代码 #

如果需要绘制3D磁场矢量图，可追加以下代码：

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')

# 绘制归一化后的磁场矢量
ax.quiver(xv, yv, zv, Bx, By, Bz, length=0.2, normalize=True)

ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
ax.set_title('Solenoid Magnetic Field')
plt.show()

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