如何创建 4D 复杂曲面图？

Question

我有以下 Matlab 代码，我想将其转换为 Python 3 代码。

r = (0:1:15)';                           % create a matrix of complex inputs
theta = pi*(-2:0.05:2);
z = r*exp(1i*theta);
%w = z.^(1/2)  ;                          % calculate the complex outputs
w = sqrt(r)*exp(1i*theta/2);

figure('Name','Graphique complexe','units','normalized','outerposition',[ 0.08 0.1 0.8 0.55]);
subplot(121)

surf(real(z),imag(z),real(w),imag(w))    % visualize the complex function using surf
xlabel('Real(z)')
ylabel('Imag(z)')
zlabel('Real(u)')
cb = colorbar;
colormap jet;                            % gradient from blue to red
cb.Label.String = 'Imag(v)';

subplot(122)
surf(real(z),imag(z),imag(w),real(w))    % visualize the complex function using surf
xlabel('Real(z)')
ylabel('Imag(z)')
zlabel('Imag(v)')
cb = colorbar;
colormap jet;                            % gradient from blue to red
cb.Label.String = 'Real(u)';

结果和原始讨论可以在这里找到。在这个SO page上还有一个讨论。但是，我未能运行和重现这些代码。我接下来可以尝试什么？

Answer 1

如果您花时间学习 matplotlib（尤其是 3d 轴）的工作原理，这将非常简单：

import numpy as np  
import matplotlib.pyplot as plt  
import matplotlib.cm as cm 
from mpl_toolkits.mplot3d import Axes3D 
 
# compute data to plot 
r, theta = np.mgrid[1:16, -2*np.pi:2*np.pi:50j] 
z = r * np.exp(1j*theta)  
w = np.sqrt(r) * np.exp(1j*theta/2)  
 
# plot data  
fig = plt.figure()  
for plot_index in [1, 2]: 
    if plot_index == 1: 
        z_data, c_data = w.real, w.imag 
        z_comp, c_comp = 'Re', 'Im' 
    else: 
        z_data, c_data = w.imag, w.real 
        z_comp, c_comp = 'Im', 'Re' 
    c_data = (c_data - c_data.min()) / c_data.ptp() 
    colors = cm.viridis(c_data) 
 
    ax = fig.add_subplot(f'12{plot_index}', projection='3d') 
    surf = ax.plot_surface(z.real, z.imag, z_data, facecolors=colors,
                           clim=[z_data.min(), z_data.max()])
    ax.set_xlabel('$Re z$')  
    ax.set_ylabel('$Im z$')   
    ax.set_zlabel(f'${z_comp} w$')  
    cb = plt.colorbar(surf, ax=ax)  
    cb.set_label(f'${c_comp} w$')  
 
plt.show()

结果：

应该注意的一些事项：

• Viridis 颜色图很好，jet 很差。
• 一般来说，复杂（互锁）3d 几何图形可能存在渲染问题，因为 matplotlib 有一个 2d 渲染器。幸运的是，在这种情况下，数据集耦合得足够紧密，即使您以交互方式围绕图形旋转，这似乎也不会发生。（但如果您将两个相交的曲面绘制在一起，情况可能会有所不同。）
• 人们可能希望启用标签的乳胶渲染，以使结果更加酥脆。
• 如果根据数据的 z 分量使用着色的默认选项，阴影看起来会好很多。

如果我们还想移植我的 MATLAB 答案的第二部分，您将不得不使用一种技巧将函数的两个分支缝合在一起（正如我所说，这对于正确渲染互锁表面是必要的）。对于上面代码中的具体示例，这不会给您完美的结果，因为两个分支本身在虚部中都包含不连续性，因此无论我们努力很好地渲染两个表面，结果看起来都会有点糟糕：

import numpy as np 
import matplotlib.pyplot as plt 
import matplotlib.cm as cm 
from mpl_toolkits.mplot3d import Axes3D 
 
# compute data to plot 
r0 = 15 
re, im = np.mgrid[-r0:r0:31j, -r0:r0:31j] 
z = re + 1j*im 
r, theta = abs(z), np.angle(z) 
w1 = np.sqrt(r) * np.exp(1j*theta/2)  # first branch 
w2 = np.sqrt(r) * np.exp(1j*(theta + 2*np.pi)/2)  # second branch 
 
# plot data 
fig = plt.figure() 
for plot_index in [1, 2]: 
    # construct transparent bridge 
    re_bridge = np.vstack([re[-1, :], re[0, :]]) 
    im_bridge = np.vstack([im[-1, :], im[0, :]]) 
    c_bridge = np.full(re_bridge.shape + (4,), [1, 1, 1, 0])  # 0% opacity
 
    re_surf = np.vstack([re, re_bridge, re]) 
    im_surf = np.vstack([im, im_bridge, im]) 
    w12 = np.array([w1, w2]) 
    if plot_index == 1: 
        z_comp, c_comp = 'Re', 'Im' 
        z12, c12 = w12.real, w12.imag 
    else: 
        z_comp, c_comp = 'Im', 'Re' 
        z12, c12 = w12.imag, w12.real 
         
    color_arrays = cm.viridis((c12 - c12.min()) / c12.ptp()) 
    z1,z2 = z12 
    c1,c2 = color_arrays 
     
    z_bridge = np.vstack([z1[-1, :], z2[0, :]]) 
    z_surf = np.vstack([z1, z_bridge, z2]) 
    c_surf = np.vstack([c1, c_bridge, c2]) 
     
    ax = fig.add_subplot(f'12{plot_index}', projection='3d') 
    surf = ax.plot_surface(re_surf, im_surf, z_surf, facecolors=c_surf, 
                           clim=[c12.min(), c12.max()], 
                           rstride=1, cstride=1) 
    ax.set_xlabel('$Re z$') 
    ax.set_ylabel('$Im z$') 
    ax.set_zlabel(f'${z_comp} w$') 
    cb = plt.colorbar(surf, ax=ax) 
    cb.set_label(f'${c_comp} w$') 
  
plt.show()

右图中丑陋的跳跃可能需要大量工作才能解决，但这并不容易：这是两个表面数据集中在负实参数下发生的实际不连续性。由于您的实际问题可能更像这样，您可能不需要面对这个问题，您可以使用上述拼接（桥接）技巧来组合您的曲面。