已设置x、r非负，绘制圆及旋转体（球体）仍报ufunc循环错误求解

核心错误原因 #

• SymPy符号与NumPy数值环境混用：代码中r始终是SymPy的符号对象，未赋予具体数值，但绘图部分使用NumPy的数值函数（如np.sqrt）对其运算，NumPy无法处理符号类型，直接触发报错。
• 变量名冲突：先定义了SymPy符号x，后续又将x赋值为NumPy数组，虽不会直接报错，但会导致逻辑混淆，增加调试难度。
• 积分结果无法数值化：vxn = vxf.evalf()因r是符号，无具体值，无法得到有效数值，后续标题格式化会出现异常。

修正后的代码 #

# Calculate the surface area of y = sqrt(r^2 - x^2)
# revolved about the x-axis

import matplotlib.pyplot as plt
import numpy as np
import sympy as sy

# SymPy符号计算部分，使用独立变量名避免冲突
x_sym = sy.Symbol("x", nonnegative=True)
r_sym = sy.Symbol("r", nonnegative=True)

def f(x):
    return sy.sqrt(r_sym**2 - x**2)

def fd(x):
    return sy.simplify(sy.diff(f(x), x))

def vx(x):
    return 2*sy.pi*(f(x)*sy.sqrt(1 + (fd(x) ** 2)))
  
vxi = sy.Integral(vx(x_sym), (x_sym, -r_sym, r_sym))
vxf = vxi.simplify().doit()

# 给r赋予具体数值，用于数值计算和绘图
r_val = 2
# 计算表面积的数值结果
surface_area = vxf.subs(r_sym, r_val).evalf()

n = 100

fig = plt.figure(figsize=(14, 7))
ax1 = fig.add_subplot(221)
ax2 = fig.add_subplot(222, projection='3d')
ax3 = fig.add_subplot(223)
ax4 = fig.add_subplot(224, projection='3d')

# 绘图部分使用NumPy数组，变量名与SymPy符号区分
x_np = np.linspace(-r_val, r_val, 100)  # 修正x的范围，符合圆的定义域
# Plot the circle
y_np = np.sqrt(r_val ** 2 - x_np ** 2)
t = np.linspace(0, np.pi * 2, n)

xn = np.outer(x_np, np.cos(t))
yn = np.outer(x_np, np.sin(t))
zn = np.zeros_like(xn)
for i in range(len(x_np)):
    zn[i:i + 1, :] = np.full_like(zn[0, :], y_np[i])

ax1.plot(x_np, y_np)
ax1.set_title("$f(x)$")
ax2.plot_surface(xn, yn, zn)
ax2.set_title("$f(x)$: Revolution around $y$")

# 修正反函数计算逻辑
y_inverse = np.linspace(0, r_val, 100)
x_inverse = np.sqrt(r_val ** 2 - y_inverse ** 2)
xn_inverse = np.outer(x_inverse, np.cos(t))
yn_inverse = np.outer(x_inverse, np.sin(t))
zn_inverse = np.zeros_like(xn_inverse)
for i in range(len(x_inverse)):
    zn_inverse[i:i + 1, :] = np.full_like(zn_inverse[0, :], y_inverse[i])

ax3.plot(x_inverse, y_inverse)
ax3.set_title("Inverse of $f(x)$")
ax4.plot_surface(xn_inverse, yn_inverse, zn_inverse)
ax4.set_title("$f(x)$: Revolution around $x$ \n Surface Area = {:.2f}".format(surface_area))

plt.tight_layout()
plt.show()

关键修正点说明 #

• 把SymPy符号变量命名为x_sym、r_sym，绘图用的NumPy数组命名为x_np、r_val，彻底避免变量冲突。
• 给r赋予具体数值r_val=2，让NumPy能进行正常的数值运算。
• 修正x_np的范围为[-r_val, r_val]，符合圆的定义域（原代码[1,3]在r=2时会出现根号内负数的问题）。
• 反函数部分重新定义y_inverse的范围，确保计算合法。
• 用subs方法替换符号为具体数值后，再调用evalf()得到有效的表面积数值。

内容的提问来源于stack exchange，提问作者Freya the Goddess