Hi everyone,
I am trying to integrate polynomials on an n-sphere and implemented (see
attached) the algorithm published in Alan Genz "Fully Symmetric
Interpolatory Rules for Multiple Integrals over Hyper-Spherical Surfaces"
Journal of Computational and Applied Mathematics 157: 187-195, 2003. On
double precision, everything is fine; however, I lose precision in some
integrals using more than double prec.
I am using randomized samples on the sphere which (by construction) should
integrate the polynomials (that I tested) exactly. However, integrating all
polynomials up to degree 6 on a 5-dimensional sphere with points that
should integrate all polynomials up to degree 7 exactly, I obtain the
following test results :
mp.prec = 1024 (aka mp.dps = 307)
polynomial | degree of polynomial | log_10 ( difference of two
randomized integrations ) | log_10 ( relative error of weights ) |
analytic value of integral
[0, 0, 0, 0, 0, 0] 0 -inf -306.962262
31.0062767
[1, 0, 0, 0, 0, 0] 1 -308.111506 -306.962262 0.0
[2, 0, 0, 0, 0, 0] 2 -14.519354 -306.962262 5.16771278
[1, 1, 0, 0, 0, 0] 2 -14.8787358 -306.962262 0.0
[3, 0, 0, 0, 0, 0] 3 -308.583548 -306.962262 0.0
[2, 1, 0, 0, 0, 0] 3 -309.157806 -306.962262 0.0
[1, 1, 1, 0, 0, 0] 3 -309.759866 -306.962262 0.0
[4, 0, 0, 0, 0, 0] 4 -14.6442927 -306.962262 1.93789229
[3, 1, 0, 0, 0, 0] 4 -15.3047045 -306.962262 0.0
[2, 2, 0, 0, 0, 0] 4 -15.0673319 -306.962262 0.645964098
[2, 1, 1, 0, 0, 0] 4 -16.0594155 -306.962262 0.0
[1, 1, 1, 1, 0, 0] 4 -31.8146065 -306.962262 0.0
[5, 0, 0, 0, 0, 0] 5 -309.37189 -306.962262 0.0
[4, 1, 0, 0, 0, 0] 5 -309.481766 -306.962262 0.0
[3, 2, 0, 0, 0, 0] 5 -309.739262 -306.962262 0.0
[3, 1, 1, 0, 0, 0] 5 -310.095075 -306.962262 0.0
[2, 2, 1, 0, 0, 0] 5 -310.474598 -306.962262 0.0
[2, 1, 1, 1, 0, 0] 5 -310.714154 -306.962262 0.0
[1, 1, 1, 1, 1, 0] 5 -311.46643 -306.962262 0.0
[6, 0, 0, 0, 0, 0] 6 -14.7692314 -306.962262 0.968946146
[5, 1, 0, 0, 0, 0] 6 -15.6057345 -306.962262 0.0
[4, 2, 0, 0, 0, 0] 6 -15.4314091 -306.962262 0.193789229
[3, 3, 0, 0, 0, 0] 6 -15.8275833 -306.962262 0.0
[4, 1, 1, 0, 0, 0] 6 -16.5822942 -306.962262 0.0
[3, 2, 1, 0, 0, 0] 6 -16.4753651 -306.962262 0.0
[2, 2, 2, 0, 0, 0] 6 -16.0997045 -306.962262
0.0645964098
[3, 1, 1, 1, 0, 0] 6 -32.3374853 -306.962262 0.0
[2, 2, 1, 1, 0, 0] 6 -17.3552924 -306.962262 0.0
[2, 1, 1, 1, 1, 0] 6 -33.4100353 -306.962262 0.0
[1, 1, 1, 1, 1, 1] 6 -49.1751847 -306.962262 0.0
It should be noted that the polynomial are written in the following form: a
list p of length n represents the polynomial
z_0^{p[0]} z_1^{p[1]} z_2^{p[2]} ... z_{n-1}^{p[n-1]}.
If we look at the second to last column, we can see that the weights add
are exact with respect to mp.prec=1024 but do not add up to 1.0 exactly.
Hence, the loss in precision is unlikely to be caused by the weights being
only on double prec (and even if it is, then I don't understand why because
they are calculated with mp.prec=1024, as well).
The really interesting column is the middle one. Here, I took two
randomized set of integration points (both should integrate all polynomials
exactly up to machine error) and printed the decadic logarithm of the
absolute value of their difference. In other words, we are expecting the
middle column to be populated with numbers in the vicinity of -307 (just
like the fourth column). However, we only get these values for odd degree
polynomials (and the constant 1 but that is just the sum of weights, i.e.,
the randomization doesn't do anything). For even degree polynomials, we
have significantly reduced precision.
Does anyone have an idea what might be the reason for such behavior?
Thank you very much,
Tobias
--
You received this message because you are subscribed to the Google Groups
"sympy" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
To post to this group, send email to [email protected].
Visit this group at https://groups.google.com/group/sympy.
To view this discussion on the web visit
https://groups.google.com/d/msgid/sympy/cf76f223-96a0-4789-8530-db224622f6bb%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
import sympy as smp
import numpy
from copy import deepcopy
mp = smp.mpmath
prec = 1024 # set precision of computation in bit
mp.mp.prec = prec # do not set this parameter directly as it will be changed for output purposes and reset to prec above
# define function to cycle partitions
def next_partition(old_part):
new_part = deepcopy(old_part)
part_sum = new_part[0]
for i in range(1,len(old_part)):
part_sum += new_part[i]
if new_part[0] <= new_part[i]+1:
new_part[i] = 0
else:
new_part[0] = part_sum - i * (new_part[i]+1)
for k in range(1,i+1):
new_part[k] = new_part[i]+1
return new_part
new_part[0] = part_sum + 1
return new_part
# find all partitions
def get_all_partitions(m,dim):
part = [m]+[0]*(dim-1)
while sum(part) == m:
yield part
part = deepcopy(next_partition(part))
# volume of S^{dim-1} - 2 \Gamma(1/2)^{dim} / \Gamma(dim/2)
def vol_factor(dim):
return mp.mpmathify(2.)*mp.power(mp.gamma(.5),mp.mpmathify(dim))/mp.gamma(mp.mpmathify(dim)/2.)
# normalized weights
def get_weight(s,g,d):
w_s = mp.mpmathify(0.)
k = [0]*d
while True:
i_s = 0
k_s = 0
m_c = 0
k_n = mp.mpmathify(1.)
p_d = mp.mpmathify(1.)
for i in range(d):
i_s += 1
if k[i] == 0:
p_d *= mp.mpmathify(1.-i_s)
else:
p_d *= mp.mpmathify(d*k_n)/mp.mpmathify(s+k_s)
k_s += 2
k_n += 2
if i_s == g[m_c]:
m_c += 1
i_s = 0
k_n = 1
w_s += p_d
for i in range(d):
k[i] += 1
if k[i] < 2:
break
k[i] = 0
if i == d-1 and k[d-1] == 0:
break
for i in range(s):
for j in range(g[i]):
w_s = w_s/mp.mpmathify( g[i]-j )
c = 0
j = 0
while j < len(g) and g[j] > 0:
c += 1
j += 1
return mp.power(mp.mpmathify(2.),mp.mpmathify(-c))*w_s
# calculate all permutations (including identical permutations, i.e., [1,1,0] will appear twice)
def permute(l,low=0):
if low+1 >= len(l):
yield l
else:
for l1 in permute(l,low+1):
yield l1
for i in range(low+1,len(l)):
l[low],l[i]=l[i],l[low]
for l1 in permute(l,low+1):
yield l1
l[low],l[i]=l[i],l[low]
# calculate all sign combinations of non-zero elements
def signs(p,low=0):
if low >= len(p):
yield p
else:
for p1 in signs(p,low+1):
yield p1
if not (p[low]==0.):
p[low] *= -1.0
for p1 in signs(p,low+1):
yield p1
p[low] *= -1.0
# calculate a random orthogonal Matrix
def random_orthogonal(n):
O = numpy.zeros((n,n))
while numpy.linalg.det(O) == 0:
O = numpy.random.normal(size=(n,n)) #generate nxn matrix with normally distributed entries and det =/= 0
Q,R = numpy.linalg.qr(O) # QR decomposition: Q is obtained by Gram Schmidt on the columns of O
return Q
class spherical_rule: # Genz's method for S^k with k > 2; for k = 2, use t-design
def __init__(self,dim,deg,randomized=True,normalized = False):
self.__dim = dim # dimension of space
self.__m = max(int(deg/2),1) # degree parameter m
self.__deg = 2*self.__m+1 # exact up to degree self.__deg
self.__random = randomized # shall points be randomized? default = yes
self.__normalized = normalized # shall weights be normalized? default = no
self.__vol_fact = vol_factor(dim) # get volume of sphere
self.__random_orthogonal = mp.matrix(random_orthogonal(dim)) # get a randomized orthogonal matrix
self.points = [] # calculate all points of the spherical rule
for p in self.__generate_points():
if p not in self.points:
self.points.append(p)
def __generate_points(self):
if self.__dim == 2: # for S^2 use the k-th unit roots as (k-1)-design
phase = mp.mpmathify(1.)
if self.__random: # if randomizied, add a randomly chosen phase shift
phase = mp.exp(mp.mpmathify(2.j)*mp.pi*mp.rand())
w = mp.mpmathify(1.)/mp.mpmathify(self.__deg + 1.) # weights of k-point t-design: 1/k
if not self.__normalized:
w *= self.__vol_fact # normalized weights to be multiplied by volume of sphere for spherical integration
for root in mp.unitroots(self.__deg + 1): # cycle all unit roots
x = phase * root # add phase shift
yield [ mp.matrix([[x.real],[x.imag]]) , w ] # return point (vector) and weight in a list
else: # for all other dimensions, use Genz
for partition in get_all_partitions(self.__m, self.__dim): # generate all partitions of m in dim summands -> generates points
w = get_weight(self.__dim, partition, self.__m) # calculate corresponding weight (normalized)
if not self.__normalized: # if spherical integration, multiply by volume of sphere
w *= self.__vol_fact
# let g be the generator, then the component x[i] of the raw point associated with g is given by sqrt( g[i]/m )
p_raw = []
for i in range(len(partition)):
p_raw.append(mp.sqrt( mp.mpmathify(partition[i])/mp.mpmathify(self.__m) ))
for p_permuted in permute(p_raw): # raw point needs to be permuted
for p in signs(p_permuted): # and signs need to be swapped
if self.__random: # randomize if wanted by multiplying point with orthogonal matrix we initialized
yield [ self.__random_orthogonal*mp.matrix(p) , w ] # return randomized point and weight
else:
yield [ mp.matrix(p) , w ] # return non-randomized point and weight otherwise
# forcing the program to terminate
#####################################
##### test polynomial exactness #####
#####################################
# calculate \int_{S^{n-1}} z_1^{a_1} ... z_n^{a_n} d\vol_{S^{n-1}}(z)
# result is zero if at least one a_j is odd
# otherwise let b_j := (a_j+1)/2. then \int = 2 \Gamma(b_1) ... \Gamma(b_n) / \Gamma(b_1+...+b_n) = 2 / \Gamma(\sum_{j=1}^n b_j) * \prod_{j=1}^n\Gamma(b_j)
def analytical(powers):
r = mp.mpmathify(2.)
beta = mp.mpmathify(0.)
for j in range(len(powers)):
r *= mp.mpmathify((powers[j]+1)%2) * mp.gamma(mp.mpmathify(powers[j]+1.)/mp.mpmathify(2.))
beta += mp.mpmathify(powers[j]+1.)/mp.mpmathify(2.)
return r/mp.gamma(beta)
# run a test integrating all monomials of degree <= deg over S^{dim-1}
def test(dim,deg,rand):
points = spherical_rule(dim,deg,rand,False).points # calculate points and weights
for degree in range(deg+1): # run over all monomial degrees
for p in get_all_partitions(degree,dim): # get all monomials (ordered by highest power since everything is symmetric)
integral = mp.mpmathify(0.)
weights = mp.mpmathify(0.)
for point in points:
int_add = mp.mpmathify(1.)
for i in range(dim):
int_add *= mp.power(point[0][i,0],p[i]) # calculate product z_i^{p_i}
integral += int_add * point[1] # multiply with weight and sum up
weights += point[1] # sum up weights
"""
if analytical(p) == 0.:
print "abs",p, mp.log(abs(integral-analytical(p)),10) #integral,analytical(p),abs(integral-analytical(p)) # divide by sum of weights and multiply with volume of the sphere - absolute error (analytic value is zero)
else:
print "rel",p, mp.log(abs(integral-analytical(p))/analytical(p),10) #integral,analytical(p),(integral-analytical(p))/analytical(p) # divide by sum of weights and multiply with volume of the sphere - relative error (analytic non-zero)
"""
yield p,integral,weights,analytical(p) # return the list of powers, the integral, the sum of weights, and the analytical integral
if __name__ == "__main__":
print mp.mp
dim = 6
deg = 6
rand = True
y = test(dim,deg,rand)
for x in test(dim,deg,rand): # getting two randomized quadrature rules and run the test
z = next(y) # cycling the generator y
a,b = mp.log(abs(x[1]-z[1]),10),mp.log(abs(x[2]-vol_factor(dim))/vol_factor(dim),10) # get log (decadic) of difference of the two integations and of the relative error of the weights
with mp.workprec(32): # reduce precision to single (otherwise output in next line will be hard to read)
print x[0],sum(x[0]),a,b,x[3] # print list of powers, degree of monomial, the two logarithms above, and the analytic value of the integral
