如何在Python中绘制三元函数Q(x,y,z)的4D可视化图形？

我定义了一个接收x、y、z三个输入参数的函数Q，希望在Python中实现Q关于x、y、z的4D可视化绘图。以下是我的代码：

from scipy.special import ellipk, ellipe, ellipkm1
from numpy import pi, sqrt, linspace
from pylab import plot, xlabel, ylabel, suptitle, legend, show
from mpl_toolkits import mplot3d
import numpy as np

uo = 4*pi     # 磁导率常数 - 单位为毫高斯
u=uo        
N=18
i=0.4
a=0.413
#z=0.214
r=  lambda x,y: sqrt(x**2+y**2) 
k = lambda x, y, z, a: sqrt((4*sqrt(x**2+y**2)*a)/(a**2+(sqrt(x**2+y**2+z**2))**2+2*a*sqrt(x**2+y**2)))  # k = f(测量点距离, 线圈半径)
K = lambda k: ellipk(k**2.0)          
E = lambda k: ellipe(k**2.0)  
 
def Q(x,y,z):
    return np.where(((x==0)) & (y==0), 0, (((N*u*i)/(2.0*pi))*(x*z)/(((x**2+y**2)*sqrt(a**2+x**2+y**2+z**2+2*a*sqrt(x**2+y**2)))))*\
                    ((E(k(x,y,z,a))*(a**2+x**2+y**2+z**2)/(a**2+x**2+y**2+z**2-2*a*sqrt(x**2+y**2))) - K(k(x,y,z,a))))

4D可视化实现方案 #

4D可视化的核心是用颜色表示第四个维度（即Q的取值），以下是两种常用的实现方式：

方式一：3D散点图 #

生成x、y、z的网格点，计算每个点的Q值，通过颜色映射展示Q的大小：

import matplotlib.pyplot as plt

# 生成网格数据
x = np.linspace(-1, 1, 50)
y = np.linspace(-1, 1, 50)
z = np.linspace(-1, 1, 50)
X, Y, Z = np.meshgrid(x, y, z)

# 计算Q值
Q_vals = Q(X, Y, Z)

# 绘制3D散点图
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')

# 用scatter绘制，c参数绑定Q值作为颜色，cmap指定配色方案
scatter = ax.scatter(X.flatten(), Y.flatten(), Z.flatten(), c=Q_vals.flatten(), cmap='viridis')

# 添加标签和颜色条
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
plt.colorbar(scatter, label='Q(x,y,z)')
plt.title('Q(x,y,z)的4D可视化（3D散点图）')
plt.show()

方式二：多子图3D曲面图（固定单一变量） #

固定z的不同取值，绘制x-y平面上的Q值曲面，通过多个子图展示z变化时Q的分布：

fig = plt.figure(figsize=(15, 10))
z_values = [-0.8, -0.4, 0, 0.4, 0.8]  # 选取几个固定的z值

for idx, z_fixed in enumerate(z_values, 1):
    # 生成x-y网格
    X, Y = np.meshgrid(x, y)
    # 计算固定z时的Q值
    Q_fixed = Q(X, Y, z_fixed)
    
    # 创建子图
    ax = fig.add_subplot(2, 3, idx, projection='3d')
    surf = ax.plot_surface(X, Y, Q_fixed, cmap='viridis')
    
    ax.set_xlabel('X')
    ax.set_ylabel('Y')
    ax.set_zlabel('Q')
    ax.set_title(f'Q(x,y,z={z_fixed})')
    plt.colorbar(surf, ax=ax, shrink=0.5)

plt.tight_layout()
plt.show()

优化建议 #

原函数Q中x==0和y==0的判断可能因浮点数精度问题失效，建议改为判断x**2 + y**2 < 1e-8避免除以零错误：

def Q(x,y,z):
    return np.where(x**2 + y**2 < 1e-8, 0, (((N*u*i)/(2.0*pi))*(x*z)/(((x**2+y**2)*sqrt(a**2+x**2+y**2+z**2+2*a*sqrt(x**2+y**2)))))*\
                    ((E(k(x,y,z,a))*(a**2+x**2+y**2+z**2)/(a**2+x**2+y**2+z**2-2*a*sqrt(x**2+y**2))) - K(k(x,y,z,a))))

内容的提问来源于stack exchange，提问作者SAHAR GOHARSHENASAN