旋转体体积积分计算：绘图验证、负结果解析及工具对比问询

1. 绘图正确性分析 #

• 二维图ax1：正确。绘制了直线$x=2y$（对应$y=x/2$）和曲线$x=(y/2)^{2$（对应$y}2=4x$）在$y∈[0,8]$区间的图像，两者交点为$(0,0)$和$(16,8)$，围成区域符合要求。
• 二维图ax3：标题与标注错误。当前绘制的是$x$作为$y$的函数（即原函数的横轴纵轴互换），并非严格意义上的反函数，建议将标题改为“$x$关于$y$的函数”，label对应修改为“$x=2y$”和“$x=y^2/4$”。
• 三维图ax4（绕x轴旋转）：正确。通过参数方程绘制了$Y^2+Z2=(2\sqrt{x})^2$和$Y2+Z^2=(x/2)2$两个旋转曲面，准确呈现了绕x轴旋转的空心旋转体。
• 三维图ax2（绕y轴旋转）：错误。绕y轴旋转应以$y$为参数，构造$X=x(y)\cosθ$、$Z=x(y)\sinθ$、$Y=y$的曲面，但当前代码误用$x$作为参数，混淆了坐标轴变量，导致旋转体形状完全错误。

2. 体积计算结果验证 #

• 绕x轴体积Vx：正确。手动计算得：
  $$
  V_x = \pi\int_{0}^{16} \left[(2\sqrt{x})^2 - (\frac{x}{2})^2\right]dx = \pi\int_{0}^{16}(4x - \frac{x^2}{4})dx = \frac{512\pi}{3} \approx 536.17
  $$
  与代码中sympy计算结果一致。
• 绕y轴体积Vy：错误。代码中误用$x$作为积分变量，而绕y轴旋转应针对$y$积分，正确计算式为：
  $$
  V_y = \pi\int_{0}^{8} \left[(2y)^2 - (\frac{y^2}{4})2\right]dy = \pi\int_{0}^{8}(4y2 - \frac{y^4}{16})dy = \frac{4096\pi}{15} \approx 857.86
  $$
  需修改代码中积分变量为$y$，重新定义对应的被积函数。

3. 负体积相关问题解答 #

• 负体积成因：当被积函数写成$r1(x)^2 - r2(x)^2$时，在$x∈[0,16]$区间内，$r2(x)=2\sqrt{x} ≥ r1(x)=x/2$，导致被积函数为负，积分结果自然为负。
• 是否使用绝对值：是。体积是正标量，若不小心搞反半径顺序，取绝对值可得到正确体积；更严谨的做法是确保被积函数为外半径平方减去内半径平方（大值减小值）。
• 负体积的实际意义：几何上无实际意义，仅代表积分时边界顺序搞反，结果与真实体积大小相等、符号相反。
• 积分计算体积的要点：
  • 明确旋转轴，选择圆盘/ washer/壳层法等合适方法。
  • 确定积分变量与上下限（由曲线交点确定）。
  • 区分内外半径：保证被积函数为外边界到旋转轴距离的平方减去内边界的平方，确保非负。
  • 验证区间内函数大小：确认哪个边界离旋转轴更远，避免顺序错误。
  • 绕y轴旋转时，若用x作为函数，需使用$x$关于$y$的表达式，积分变量为$y$。

4. Sympy与Numpy/Scipy在积分计算的优劣 #

• Sympy优势：
  • 支持符号积分，能输出精确的解析解（如$\frac{512\pi}{3}$），而非数值近似值。
  • 适合数学推导，可展示积分过程、化简表达式，便于理解数学原理。
  • 支持符号变量运算，无需提前代入数值。
• Numpy/Scipy优势：
  • 数值积分效率高，适合复杂无解析解的函数、大规模数据或实时计算场景。
  • 天然支持数组运算，适配工程领域的数值计算需求。
• 总结：追求精确解析结果或数学推导时选Sympy；需要快速数值近似或处理实际数据时选Numpy/Scipy，二者可结合使用（如用Sympy求解析解，转数值用于可视化）。

原代码（MWE） #

# Compare the plot at xy axis with the solid of revolution toward x and y axis
# For region bounded by the line x - 2y = 0 and y^2 = 4x
# Plotting the revolution of the bounded region 
# can be done by limiting the np.linspace of the y, u, and x_inverse values 
# You can determine the limits by finding the intersection points of the two functions.

import matplotlib.pyplot as plt
import numpy as np
import sympy as sy
  
def r1(x):
    return x/2

def r2(x):
    return 2*(x**(1/2))

def r3(x):
    return 2*x

def r4(x):
    return (x/2)**(2)

def vx(x):
    return np.pi*(r2(x)**2 - r1(x)**2)

def vy(x):
    return np.pi*(r3(x)**2 - r4(x)**2)
  
x = sy.Symbol("x")
vx = sy.integrate(vx(x), (x, 0, 16))
vy = sy.integrate(vy(x), (x, 0, 8))

n = 200

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')

y = np.linspace(0, 8, n)
x1 = (2*y)
x2 = (y / 2) ** (2)
t = np.linspace(0, np.pi * 2, n)

u = np.linspace(0, 16, n)
v = np.linspace(0, 2 * np.pi, n)
U, V = np.meshgrid(u, v)

X = U
Y1 = (2 * U ** (1/2)) * np.cos(V)
Z1 = (2 * U ** (1/2)) * np.sin(V)

Y2 = (U / 2) * np.cos(V)
Z2 = (U / 2) * np.sin(V)

Y3 = ((U / 2) ** (2)) * np.cos(V)
Z3 = ((U / 2) ** (2)) * np.sin(V)

Y4 = (2*U) * np.cos(V)
Z4 = (2*U) * np.sin(V)

ax1.plot(x1, y, label='$y=x/2$')
ax1.plot(x2, y, label='$y=2 \sqrt{x}$')
ax1.legend()

ax1.set_title('$f(x)$')
ax2.plot_surface(X, Y3, Z3, alpha=0.3, color='red', rstride=6, cstride=12)
ax2.plot_surface(X, Y4, Z4, alpha=0.3, color='blue', rstride=6, cstride=12)
ax2.set_title("$f(x)$: Revolution around $y$ \n Volume = {}".format(vy))

# find the inverse of the function
x_inverse = np.linspace(0, 8, n)
y1_inverse = np.power(2*x_inverse, 1)
y2_inverse = np.power(x_inverse / 2, 2)

ax3.plot(x_inverse, y1_inverse, label='Inverse of $y=x/2$')
ax3.plot(x_inverse, y2_inverse, label='Inverse of $y=2 \sqrt{x}$')
ax3.set_title('Inverse of $f(x)$')
ax3.legend()

ax4.plot_surface(X, Y1, Z1, alpha=0.3, color='red', rstride=6, cstride=12)
ax4.plot_surface(X, Y2, Z2, alpha=0.3, color='blue', rstride=6, cstride=12)
ax4.set_title("$f(x)$: Revolution around $x$ \n Volume = {}".format(vx))

plt.tight_layout()
plt.show()

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