Assume I know that I can express an approximation of a function by $$ \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) $$ where $\theta$ is the inclination and $\phi$ azimuth, and $z_{l,m}^{1,2}$ are a sequence of indepedent standard normal random numbers, and $A_l$ is some pre-decided number so that it decays fast enough. The functions $L_{l,m}$ are given by $$ L_{l,m}=\sqrt{\frac{2l+1}{4\pi}\frac{(l-m)!}{(l+m)!}}P_{l,m}(\cos(\theta)) $$ where $P_{l,m}$ are the associated Lagrange polynomial

of the first type. This I know. I would like to evaluate this sum for

1) Large values of the truncation parameter $k$ (like 10^6)

2) many many values of $\theta$ and $\phi$.

I would like to do this in python, but i suppose I would be ok with writing a compiled program in Fortran for instance, as long as I can acess from python.

My attempt so far is to use scipy.special.lpmv and just blindly evaluate the sum, but that breaks down very fast for not very large values of l,m. There seems to be some sort of trouble with the scipy implementation, see.


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
def rand_sequence(noterms,seed):                                                                                                                  
def angular_spec(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!!  

def GRF(noterms,alpha,rands,x,y,z):
    for l in range(noterms):
        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)])

So, any suggestions? Will I have to give up my hope of being able to evaluate for large $k$? How can I speed this up. I will have a very large amount of $\theta$s and $\phi$s to evaluate in...

  • 2
    $\begingroup$ Do I understand correctly that you have an oscillatory sum with about 10^12 terms? And the terms include factorials of numbers like 10^6. If so, it is quite hopeless, I’m afraid. Usually there is an alternative representation of the function that converges much faster. A technique that is relevant to such sums is called the Watson transformation (jpier.org/PIER/pier75/08.07052502.Valagiannopoulos.pdf has an application of this method to sums of spherical harmonics). It is rather involved unfortunately. $\endgroup$ Oct 28 '19 at 18:19
  • 1
    $\begingroup$ Are the $\phi$ values equally spaced? $\endgroup$
    – vibe
    Oct 29 '19 at 4:27
  • $\begingroup$ The values of $\phi$ are not, sadly. I am essentially evaluating the sums at the vertices of a special tesselation of the sphere. $\endgroup$
    – edo
    Oct 29 '19 at 10:48
  • 1
    $\begingroup$ Your function is a linear combination of normal random variables, so it is a normal random variable itself (mean 0, variance given by a complicated sum of squares of spherical harmonics, probably without a closed form). That might simplify it, but in your code the random variables $z$'s appear to be fixed, being generated once, and not sampled randomly at every evaluation of the function, so maybe not. $\endgroup$
    – Kirill
    Oct 29 '19 at 16:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.