Plotting in spherical coordinates given the radial distance?

Question

For a homework problem we are asked to

---

Write a program that computes all the eigenvalues of the matrix A, using the Rayleigh Quotient iteration to find each eigenvalue.

i. Report the initial guess you used for each eigenvalue. Plot the error at the n iteration against the error at the n+1 iteration for each eigenvalue (in log log).

ii. Compute the Rayleigh quotient for x sampled on a discretization of the unit sphere (use spherical coordinates). Plot the result, for example with mesh (Θ,ϕ,r(Θ,ϕ)). Explain the number and location of the critical points.

---

I have completed the first part of the problem, but I am not sure I understand how to complete the second. After some research I have found various ways (Here and Here) to plot spherical coordinates in Python which has provided me some clues. Borrowing from those links and a few other sources I have

# 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)])

In the first chunk of code (Θ,ϕ) are created and translated to Cartesian coordinates, effectively a discritization of the unit circle in the R^3. This allows me then to create the x vectors (xyz above) used to calculate the Rayleigh Quotients (rqs above) at each discritization point, the r(Θ,ϕ) of the spherical coordinates.

Though, now that I have the radial distance I am not sure how to again recreate x, y, z properly for a meshgrid to plot as a surface. This might be out of the realm of StackOverflow and more for Math.Stack, but I am also not sure if this plot should end up being a "warped plane" or a "warped sphere".

In this SO answer, linked above too, I think the answer lies though. In the chunk of code

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 here I assume refers to the radial distance, though, this R is in a mesh when x, y, z are calculated. I have tried to reshape rqs above to the same dimension, but the values of the rqs do not line up with the subsequent grid and as a result produces obviously wrong plots.

I almost need a way to tie in the creation of the meshgrid with the calculation of the x. But the calculation seems to complex to be applied directly to the meshgrid..

How can I produce the spherical coordinates based plot given the radial distances?

---

EDIT: After some more searching I found this MatLab code which produces the desired plot, though I would still like to reproduce this in Python. I would like to say this MatLab code provides an overview of how to implement this in Python, but it seems to be some very old and esoteric code. Here is the plot(s) it produces

Answer 1

I don't know about the math or the Rayleigh Quotient but from what I gather, you want to calculate rqs as a function of the points of the unit sphere. For that, I would recommend to use meshgrid to generate value pairs for all Θ and ϕ. Then, since your matrix formula is defined for cartesian rather than spherical coordinates, I would convert my mesh to that and insert into the formula.

Then, finally, the result can be illustrated, using plot_surface on the unit sphere, where the (scaled) RQS data are used for 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()

This gives the following result, which seems to agree with theMatlab contours plot in the OP.

Note that there may be better ways (vectorized) to convert the X, Y, and Z meshes to vectors, but the main execution time seems to be in the plottting anyway, so I won't spend the time trying to figure out exactly how.