计算科学 - 球谐展开的有效评估 - 吾爱随笔录

假设我知道我可以用其中是倾角和方位角，是一系列独立的标准正态随机数，是一些预-确定的数字，使其衰减得足够快。函数由其中是相关的拉格朗日多项式

\sum_{l = 0}^{k} (\sqrt{A_{l}} z_{l, 0}^{1} L_{l, 0} (θ) + \sqrt{2 A_{l}} \sum_{m = 1}^{l} L_{l, m} (θ) (z_{l, m}^{1} \cos (m ϕ) + z_{l, m}^{2} \sin (m ϕ)))

$\sum_{l=0}^{k}\left( \sqrt{A_l} z_{l,0}^1 L_{l,0}(\theta)+\sqrt{2A_l}\sum_{m=1}^l L_{l,m}(\theta)(z_{l,m}^1 \cos(m\phi)+z_{l,m}^2\sin(m\phi))\right)$

θ

$\theta$

ϕ

$\phi$

z_{l, m}^{1, 2}

$z_{l,m}^{1,2}$

A_{l}

$A_l$

L_{l, m}

$L_{l,m}$

L_{l, m} = \sqrt{\frac{2 l + 1}{4 π} \frac{(l - m)!}{(l + m)!}} P_{l, m} (\cos (θ))

$L_{l,m}=\sqrt{\frac{2l+1}{4\pi}\frac{(l-m)!}{(l+m)!}}P_{l,m}(\cos(\theta))$

P_{l, m}

$P_{l,m}$

第一类。这个我知道。我想评估这笔款项

的大值（如10^6） $k$

和的许多值。 $\theta$ $\phi$

我想在 python 中执行此操作，但我想我可以在 Fortran 中编写一个编译程序，只要我可以从 python 访问。

到目前为止，我的尝试是使用 scipy.special.lpmv 并且只是盲目地评估总和，但是对于不是很大的 l，m 值，它很快就会崩溃。scipy 实现似乎存在某种问题，请参阅。

https://stackoverflow.com/questions/39249594/numerical-computation-of-the-spherical-harmonics-expansion

import numpy as np
from scipy.special import lpmv sph_harm                                                                                                                                                                        
from math import factorial
from dolfin import *
from mshr import *
from numpy import sqrt,cos,sin,pi
import numpy as np
no_terms=50
def rand_sequence(noterms,seed):                                                                                                                  
    np.random.seed(seed)
    return(np.random.normal(0,1,[noterms,noterms,2]))
def angular_spec(l,alpha):
    return((1+l)**(-alpha))
def plm(l,m,theta):
    return(sqrt((2*l+1)*factorial(l-m)/4/pi/factorial(l+m))*lpmv(m,l,cos(theta)))    #it seems to break here!!  
rands=rand_sequence(no_terms,123)

def GRF(noterms,alpha,rands,x,y,z):
    phi=np.arctan2(y,x)
    theta=np.arccos(z)
    out=0.0                                                                                                                                                                      
    for l in range(noterms):
        a=sqrt(angular_spec(l,alpha))*rands[l,0,0]*plm(l,0,theta)
        b=sqrt(2*angular_spec(l,alpha))*sum([plm(l,m,theta)(rands[l,m,0]*cos(m*phi)+rands[l,m,1]*sin(m*phi)) for m in range(1,l+1)])
        out=np.add(out,np.add(a,b))                                                                                                                                        
    return(out)

那么，有什么建议吗？我是否必须放弃能够评估大的希望？我怎样才能加快速度。我将有大量的和来评估...... $k$ $\theta$ $\phi$