在给定径向距离的情况下在球坐标中绘图？答案

【问题标题】：Plotting in spherical coordinates given the radial distance?在给定径向距离的情况下在球坐标中绘图？
【发布时间】：2018-08-05 12:35:36
【问题描述】：

对于家庭作业问题，我们被要求

编写一个程序，计算矩阵A 的所有特征值，使用瑞利商迭代来找到每个特征值。

我。报告您对每个特征值使用的初始猜测。针对每个特征值（在日志中）绘制n 迭代中的错误与n+1 迭代中的错误。

二。计算在单位球体的离散化上采样的x 的瑞利商（使用球坐标）。绘制结果，例如使用网格 (Θ,ϕ,r(Θ,ϕ))。解释关键点的数量和位置。

我已经完成了问题的第一部分，但我不确定我是否理解如何完成第二部分。经过一些研究，我找到了各种方法（Here 和Here）在 Python 中绘制球坐标，这为我提供了一些线索。从这些链接和我拥有的其他一些来源中借用

# The matrix for completness
A = np.matrix([[4,3,4], [3,1,3], [4,3,4]])
A = scipy.linalg.hessenberg(A)

num_points = 50
theta, phi = np.linspace(0, 2 * np.pi, num_points), np.linspace(0, np.pi, num_points)
x = np.sin(phi) * np.cos(theta)
y = np.sin(phi) * np.sin(theta)
z = np.cos(phi)

xyz = np.stack((x, y, z), axis=1)
rqs = np.array([np.dot(x, np.dot(A, np.conj(x).T)).item(0) for _,x in enumerate(xyz)])

在第一个代码块中，(Θ,ϕ) 被创建并转换为笛卡尔坐标，有效地区分了R^3 中的单位圆。这允许我创建x 向量（上面的xyz），用于计算每个分类点的瑞利商（上面的rqs），即球坐标的r(Θ,ϕ)。

不过，现在我有了径向距离，我不确定如何再次正确地重新创建 x, y, z 以将 meshgrid 绘制为表面。这可能超出了 StackOverflow 的范围以及 Math.Stack 的范围，但我也不确定这个图应该最终成为“翘曲的平面”还是“翘曲的球体”。

在上面链接的SO answer 中，我认为答案就在这里。在代码块中

theta, phi = np.linspace(0, 2 * np.pi, 40), np.linspace(0, np.pi, 40)
THETA, PHI = np.meshgrid(theta, phi)
R = np.cos(PHI**2)
X = R * np.sin(PHI) * np.cos(THETA)
Y = R * np.sin(PHI) * np.sin(THETA)
Z = R * np.cos(PHI)

R 在这里我假设是指径向距离，但是，当计算 x, y, z 时，这个 R 在网格中。我尝试将上面的reshaperqs 设置为相同的维度，但rqs 的值与随后的网格不对齐，因此会产生明显错误的图。

我几乎需要一种方法来将meshgrid 的创建与x 的计算联系起来。但是直接应用于meshgrid..

如何在给定径向距离的情况下生成基于球坐标的图？

编辑： 经过更多搜索后，我发现this MatLab code 产生了所需的绘图，但我仍然想在 Python 中重现它。我想说这个 MatLab 代码提供了如何在 Python 中实现它的概述，但它似乎是一些非常古老和深奥的代码。这是它产生的情节

【问题讨论】：

你想绘制哪个国家？半径是角度的函数吗？
轮廓是矩阵(M, x) = x* M x / x* x (en.wikipedia.org/wiki/Rayleigh_quotient) 的瑞利商的轮廓。它是矩阵和向量的函数。其中向量是角度的函数；也就是说，向量是R^3 中单位球体的区分。所以函数应该看起来像r(M, x(Θ,ϕ)) = x(Θ,ϕ)* M x(Θ,ϕ) / x(Θ,ϕ)* x(Θ,ϕ)。那么在x(Θ,ϕ) = xyz 和r(M, x(Θ,ϕ)) = rqs 上面的代码中。问题是我无法弄清楚如何将rqs 与角度的创建以及它们随后转换为笛卡尔坐标相结合。
但是r是一个标量值，是Θ和ϕ的函数，而M可以被认为是常数（在这个特定的计算中）？
是的 r 将是一个标量（瑞利商是一个单一的标量值）。而且我认为这是我不确定的地方，也许这里需要一些 Math.Stack。根据我的理解，我基本上需要绘制一个单位球体，作为曲面图，在球体上的每个点都有一个标量值。可以说是一种“扭曲的球体”。 r（嗯rps）表示单位球面上每个点到原点的距离。
实际上.. 它应该绘制一个单位球体，但作为热图。单位球面上给定点的颜色更强烈应该代表r 的更大值。实际上，我刚刚找到了这个 MatLab 代码 (people.eecs.berkeley.edu/~demmel/ma221_Fall04/Matlab/…)，它创建了所需的绘图。我将在 OP 中上传一张图片。

标签： python matplotlib linear-algebra

【解决方案1】：

我不知道数学或瑞利商，但根据我的收集，您想计算 rqs 作为单位球体点的函数。为此，我建议使用meshgrid 为所有Θ 和ϕ 生成值对。然后，由于您的矩阵公式是为笛卡尔坐标而不是球坐标定义的，我会将我的网格转换为该坐标并插入到公式中。

然后，最后，可以在单位球体上使用plot_surface 来说明结果，其中（缩放的）RQS 数据用于facecolor：

import numpy as np
import scipy.linalg
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
import matplotlib.pyplot as plt

# The matrix for completness
A = np.matrix([[4,3,4], [3,1,3], [4,3,4]])
A = scipy.linalg.hessenberg(A)

num_points = 25
theta = np.linspace(0, 2 * np.pi, num_points)
phi = np.linspace(0, np.pi, num_points)

THETA, PHI = np.meshgrid(theta, phi)

X = np.sin(PHI) * np.cos(THETA)
Y = np.sin(PHI) * np.sin(THETA)
Z = np.cos(PHI)

# Calculate RQS for points on unit sphere:
RQS = np.empty(PHI.shape)
for theta_pos, phi_pos in itertools.product(range(num_points), range(num_points)):
    x = np.array([X[theta_pos, phi_pos],
                  Y[theta_pos, phi_pos],
                  Z[theta_pos, phi_pos]])
    RQS[theta_pos, phi_pos] = np.dot(x, np.dot(A, np.conj(x).T))

# normalize in range 0..1
maxRQS = abs(RQS).max()
N = (RQS+maxRQS)/(2*maxRQS)

fig = plt.figure()
ax = fig.gca(projection='3d')
surf = ax.plot_surface(
    X, Y, Z, rstride=1, cstride=1,
    facecolors=cm.jet(N),
    linewidth=0, antialiased=False, shade=False)
plt.show()

这给出了以下结果，这似乎与 OP 中的 Matlab 等高线图一致。

请注意，可能有更好的方法（矢量化）将X、Y 和 Z 网格转换为矢量，但主要执行时间似乎无论如何都在绘图中，所以我不会花时间试图弄清楚到底是怎么做的。

【讨论】：