球面調和関数の値をPython3/matplotlibで３D表示する

２次元の調和関数の値を、r=1の円週上で均等にスキャンして合計するとゼロになるらしい。
３次元の調和関数の値を、r=1の球面上で均等にスキャンして合計するとゼロになるか、Pythonで試してみた。
その前に、値を３次元でmatplotlibを使い、赤青の濃淡で表示させてみた。

３次元調和関数

{{\displaystyle {\partial ^{2} \over \partial x^{2}}\phi (x,y,z)+{\partial ^{2} \over \partial y^{2}}\phi (x,y,z)+{\partial ^{2} \over \partial z^{2}}\phi (x,y,z)=0.}
}

これを満たす式で、一番簡単なのをネットで見つける。

f(x,y,z) = xyz

この式の値を　r=1　の球面上のサンプル点で値を取り、色付けをしてスキャター表示する。
まずは、球面上の点を均等にサンプリングする方法を考える。
普通に、緯度、経度で分割すると極部付近での点が多くなるので、緯線の長さに応じて点の数を調整している。

濃い赤が大きなプラス、濃い青が大きなマイナスを示しており、球面全体の平均を取るとゼロになりそうな配置である。
Python3で、実際に合計値を取るとゼロでした。
Python3でプログラムを実行してもらうと、スクロールして見れます。真ん中の白い部分は、極です。

プログラムは、こちら

harmony.py

import math
from mpl_toolkits import mplot3d
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm

#この式を、入れ替える
def f(x,y,z): 
  return x*y*z

fig = plt.figure()
ax = plt.axes(projection='3d')

c = 0
N = 50
zdata = np.linspace(0,0,2*N*N)
xdata = np.linspace(0,0,2*N*N)
ydata = np.linspace(0,0,2*N*N)
udata = np.linspace(0,0,2*N*N)

for j in range(0,N): #緯度のループ
  ps = math.pi*(j+0.5)/N
  w = int(N*math.sin(ps))
  for i in range(-1*w,w): #経度のループ
    th = math.pi*(i+0.5)/w
    # 極座標からxyz座標へ変換
    z = math.cos(ps)
    x = math.sin(ps)*math.cos(th)
    y = math.sin(ps)*math.sin(th)
    # ドット表示用データを記録
    xdata[c] = x
    ydata[c] = y
    zdata[c] = z
    u = f(x,y,z) 
    udata[c] = u #uが調和関数の値
    c = c + 1

# 3D スキャター表示
ax.set_aspect('equal')
ax.scatter3D(xdata, ydata, zdata, c=udata, cmap=cm.coolwarm);
plt.show()

関数f()を、いろいろ変えると、きれいな調和関数を３D表示できるので、次回。