Python代码检查：函数(x⁶+2)/(8x²)绕x轴旋转曲面绘制及优化

需求 #

绘制函数$f(x)=\frac{x^6 + 2}{8x^2}$在区间$1 ≤ x ≤ 3$内绕x轴旋转的曲面，以下为原Python代码：

import matplotlib.pyplot as plt
import numpy as np

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')
x = np.linspace(1, 3, 3)
y = ((x ** 6) + 2) / (8 * x ** 2)
t = np.linspace(0, np.pi * 2, n)

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

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

# find the inverse of the function
y_inverse = x
x_inverse = ((y_inverse ** 6) + 2) / ( 8 * x ** 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$")

plt.tight_layout()
plt.show()

问题解答 #

1. 绕x轴旋转曲面的绘图是否正确？ #

不正确。核心问题有两点：

• 变量逻辑错误：绕x轴旋转时，应保持x为原取值，y和z由$f(x)$绕x轴旋转形成圆周（即$y = f(x)\cos(t)$，$z = f(x)\sin(t)$），原代码错误混淆了x与y的关系，误用逆函数逻辑构造坐标。
• 采样点过少：x仅取3个点，导致曲面呈现离散片状，完全没有连续曲面的平滑效果。

2. 更优的绘图方法 #

• 提升采样密度：将x的采样点从3个增加到100个以上，保证曲线和曲面的平滑度。
• 用numpy广播替代循环：利用numpy向量运算直接生成旋转坐标，避免手动循环赋值，代码更高效简洁。
• 正确构造旋转坐标：绕x轴旋转时，直接基于原函数值生成y、z方向的圆周坐标，无需额外求逆函数（仅绕y轴旋转时需要逆函数逻辑）。
• 优化可视化细节：给坐标轴添加清晰标签，使用颜色映射（cmap）增强曲面层次感，或用plot_wireframe展示曲面轮廓。

3. 对比Julia的平滑曲线，绘图及坐标轴范围调整方案 #

• 提升曲线平滑度：将x = np.linspace(1, 3, 3)改为x = np.linspace(1, 3, 100)，让2D函数曲线和3D曲面都足够平滑。
• 调整坐标轴范围：
  • 2D绘图：固定x轴范围为[1,3]，y轴范围可通过ax.axis('tight')自动适配函数值，或手动设置为[min(y), max(y)]，让曲线充满画布。
  • 3D绘图：添加x、y、z轴标签，设置z轴范围为[0, max(y)]，确保曲面完整形态被展示，避免空区域浪费空间。
• 优化显示效果：开启3D坐标轴网格，调整视角让曲面旋转形态更直观。

修正后的代码 #

import matplotlib.pyplot as plt
import numpy as np

# 提升采样点数量保证平滑
x = np.linspace(1, 3, 100)
t = np.linspace(0, 2 * np.pi, 100)
y = (x**6 + 2) / (8 * x**2)

# 构造绕x轴旋转的3D坐标：x不变，y和z由f(x)旋转得到
X, T = np.meshgrid(x, t)
Y = y * np.cos(T)
Z = y * np.sin(T)

# 绘图
fig = plt.figure(figsize=(12, 6))

# 2D函数曲线
ax1 = fig.add_subplot(121)
ax1.plot(x, y, linewidth=2)
ax1.set_title("$f(x) = \\frac{x^6 + 2}{8x^2}$")
ax1.set_xlabel("x")
ax1.set_ylabel("f(x)")
ax1.set_xlim(1, 3)
ax1.grid(True)

# 3D旋转曲面
ax2 = fig.add_subplot(122, projection='3d')
surf = ax2.plot_surface(X, Y, Z, cmap='viridis', edgecolor='none')
ax2.set_title("Revolution around x-axis")
ax2.set_xlabel("x")
ax2.set_ylabel("y")
ax2.set_zlabel("z")
ax2.set_zlim(0, np.max(y))
fig.colorbar(surf, shrink=0.5, aspect=5)

plt.tight_layout()
plt.show()

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