【问题标题】:Numerical computation of the Spherical harmonics expansion球面谐波展开的数值计算
【发布时间】:2016-08-31 12:24:37
【问题描述】:

我正在尝试理解球面谐波展开以解决更复杂的问题,但我期望从一个非常简单的计算得到的结果是不正确的。我不知道为什么会这样。


一点理论:众所周知,球面上的函数()可以定义为一些常数系数和球面的无限和谐波:

球谐函数定义为:

其中 是关联的勒让德多项式。

最后,常数系数 可以计算如下(类似于傅里叶变换):


问题:假设我们有一个以 为中心的球体,其中表面上的函数对于所有点 都等于。我们要计算常数系数,然后通过近似计算返回表面函数。由于 常数系数的计算减少到:

在数值上(在 Python 中)可以使用:

def Ylm(l,m,theta,phi):
    return scipy.special.sph_harm(m,l,theta,phi)

def flm(l,m):
    phi, theta = np.mgrid[0:pi:101j, 0:2*pi:101j]
    return Ylm(l,m,theta,phi).sum()

然后,通过计算 上的带限总和,我希望在 时看到 对于任何给定点。

L = 20
f = 0
theta0, phi0 = 0.0, 0.0
for l in xrange(0,L+1):
    for m in xrange(-l,l+1):
        f += flm(l,m)*Ylm(l,m,theta0,phi0)
print f

但对于,它给了我 而不是。对于,它给了我

我知道这似乎更像是一个数学问题,但公式应该是正确的。问题似乎出在我的计算上。这可能是一个非常愚蠢的错误,但我无法发现它。有什么建议吗?

谢谢

【问题讨论】:

  • 我同意你的理论部分,但你所说的问题奇迹般地使用了 f^l_m 的总和而不是积分。您确定不想将总和的每个元素与元素的大小相乘吗?否则积分将取决于您用于离散化的元素数量。
  • 也许我错了,但据我了解,连续空间中的积分减少为离散空间中的总和(如离散傅里叶变换)。实际上,缩放因子可能是解决方案,但我不知道如何检索该值。我对 L=20 做的总和是 441。根据我得到的结果,该值应该在 1/12860 左右。
  • 你说得对,你可以使用和而不是积分,也许我不明白你在做什么。我的观点是,当您将积分转换为总和时,您不能只删除 d Omega 而不用某种形式的 delta Omega_ij 替换它。在您的情况下,总和的每个元素都是一个高度为 Y^l_m(...) 且底面积为 (pi/100)*(pi/2/100) 的 kuboid,因此 dOmega 约为 1/2000,即偏离 2 pi。正如 dmuir 所写,我猜想使用 sin(theta) 与球坐标进行积分会解决这个问题。
  • 这个积分似乎是我不可或缺的路径。你是这样编码的吗?我希望看到该路径的离散化以及类似高斯正交的东西来评估每个部分。
  • 是的,正如 mars 和 dmuir 指出的那样,球面积分不能转换为简单的总和。我试图解决的问题来自这篇文章lmb.informatik.uni-freiburg.de/papers/download/fe_bu_icpr08.pdf,尤其是公式 8。要么文章是错误的,要么我对它的理解是错误的。是的,正如在这个线程math.stackexchange.com/questions/1337806/… 中指出的那样,积分可以使用高斯正交进行数值计算,但我不确定如何找到系数。欧米茄

标签: python math


【解决方案1】:

球谐函数与内积正交

<f|g> = Integral( f(theta,phi)*g(theta,phi)*sin(theta)*dphi*dtheta)

所以你应该计算系数

  clm = Integral( Ylm( theta, phi) * sin(theta)*dphi*dtheta)

【讨论】:

  • 所以我应该能够使用如下总和计算系数 f:flm = Sum( Ylm( theta, phi) * sin(theta) ),它转换为 Python 代码,如:def flm(l,m): phi, theta = np.mgrid[0:np.pi:101j, 0:2*np.pi:101j] return np.multiply(Ylm(l,m,theta, phi),np.sin(theta)).sum() ?
  • 对不起,我不知道 python,但我认为你需要对其进行规范化。如果你“积分”常数 1,你应该得到 4*pi(球体的面积),所以归一化常数是 4*pi/whatever-you-get
猜你喜欢
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 2021-04-26
  • 2011-12-22
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
相关资源
最近更新 更多