Package hedge :: Module quadrature
[hide private]
[frames] | no frames]

Source Code for Module hedge.quadrature

  1  """1D quadrature for Jacobi polynomials. Grundmann-Moeller cubature on the simplex.""" 
  2   
  3  from __future__ import division 
  4   
  5  __copyright__ = "Copyright (C) 2007 Andreas Kloeckner" 
  6   
  7  __license__ = """ 
  8  This program is free software: you can redistribute it and/or modify 
  9  it under the terms of the GNU General Public License as published by 
 10  the Free Software Foundation, either version 3 of the License, or 
 11  (at your option) any later version. 
 12   
 13  This program is distributed in the hope that it will be useful, 
 14  but WITHOUT ANY WARRANTY; without even the implied warranty of 
 15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the 
 16  GNU General Public License for more details. 
 17   
 18  You should have received a copy of the GNU General Public License 
 19  along with this program.  If not, see U{http://www.gnu.org/licenses/}. 
 20  """ 
 21   
 22   
 23   
 24   
 25  import numpy 
26 27 28 29 30 -def jacobi_gauss_lobatto_points(alpha, beta, N):
31 """Compute the M{N}th order Gauss-Lobatto quadrature 32 points, x, associated with the Jacobi polynomial, 33 of type (alpha,beta) > -1 ( <> -0.5). 34 """ 35 36 x = numpy.zeros((N+1,)) 37 x[0] = -1 38 x[-1] = 1 39 40 if N == 1: 41 return x 42 43 x[1:-1] = numpy.array( 44 JacobiGaussQuadrature(alpha+1, beta+1, N-2).points 45 ).real 46 return x
47
48 49 50 51 -def legendre_gauss_lobatto_points(N):
52 """Compute the M{N}th order Gauss-Lobatto quadrature 53 points, x, associated with the Legendre polynomials. 54 """ 55 return jacobi_gauss_lobatto_points(0, 0, N)
56
57 58 59 60 -class Quadrature(object):
61 """An abstract quadrature rule."""
62 - def __init__(self, points, weights):
63 self.weights = weights 64 self.points = points 65 self.data = zip(points, weights)
66
67 - def __call__(self, f):
68 """Integrate the callable f with respect to the given quadrature rule. 69 """ 70 return sum(w*f(x) for x, w in self.data)
71
72 73 74 -class JacobiGaussQuadrature(Quadrature):
75 """An M{N}th order Gauss quadrature associated with the Jacobi 76 polynomials of type M{(alpha,beta) > -1} 77 78 C{alpha} and C{beta} may not be -0.5. 79 80 Integrates on the interval (-1,1). 81 """
82 - def __init__(self, alpha, beta, N):
83 x, w = self.compute_weights_and_nodes(N, alpha, beta) 84 Quadrature.__init__(self, x, w)
85 86 @staticmethod
87 - def compute_weights_and_nodes(N, alpha, beta):
88 """Return (nodes, weights) for an n-th order Gauss quadrature 89 with the Jacobi polynomials of type (alpha, beta). 90 """ 91 # follows 92 # Gene H. Golub, John H. Welsch, Calculation of Gauss Quadrature Rules, 93 # Mathematics of Computation, Vol. 23, No. 106 (Apr., 1969), pp. 221-230 94 # doi:10.2307/2004418 95 96 # see also doc/hedge-notes.tm for correspondence with the Jacobi 97 # recursion from Hesthaven/Warburton's book 98 99 from math import sqrt 100 101 apb = alpha+beta 102 103 # see Appendix A of Hesthaven/Warburton for these formulas 104 def a(n): 105 return ( 106 2/(2*n+apb) 107 * 108 sqrt( 109 (n*(n+apb)*(n+alpha)*(n+beta)) 110 / 111 ((2*n+apb-1)*(2*n+apb+1)) 112 ) 113 )
114 115 def b(n): 116 if n == 0: 117 return ( 118 -(alpha-beta) 119 / 120 (apb+2) 121 ) 122 else: 123 return ( 124 -(alpha**2-beta**2) 125 / 126 ((2*n+apb)*(2*n+apb+2)) 127 )
128 129 T = numpy.zeros((N+1, N+1)) 130 131 for n in range(N+1): 132 T[n,n] = b(n) 133 if n > 0: 134 T[n,n-1] = current_a 135 if n < N: 136 next_a = a(n+1) 137 T[n,n+1] = next_a 138 current_a = next_a 139 140 assert numpy.linalg.norm(T-T.T) < 1e-12 141 eigval, eigvec = numpy.linalg.eigh(T) 142 143 from numpy import dot, diag 144 assert numpy.linalg.norm(dot(T, eigvec) - dot(eigvec, diag(eigval))) < 1e-12 145 146 from hedge.polynomial import JacobiFunction 147 p0 = JacobiFunction(alpha, beta, 0) 148 nodes = eigval 149 weights = [eigvec[0,i]**2 / p0(nodes[i])**2 for i in range(N+1)] 150 151 return nodes, weights 152
153 154 155 156 -class LegendreGaussQuadrature(JacobiGaussQuadrature):
157 """An M{N}th order Gauss quadrature associated with the Legendre polynomials. 158 """
159 - def __init__(self, N):
160 JacobiGaussQuadrature.__init__(self, 0, 0, N)
161
162 163 164 165 -class TransformedQuadrature(Quadrature):
166 """A quadrature rule on an arbitrary interval M{(a,b)}. """ 167
168 - def __init__(self, quad, left, right):
169 """Transform a given quadrature rule `quad' onto an arbitrary 170 interval (left, right). 171 """ 172 self.left = left 173 self.right = right 174 175 length = right-left 176 assert length > 0 177 half_length = length / 2 178 Quadrature.__init__(self, 179 [left + (p+1)/2*length for p in quad.points], 180 [w*half_length for w in quad.weights])
181
182 183 184 185 -def _extended_euclidean(q, r):
186 """Return a tuple (p, a, b) such that p = aq + br, 187 where p is the greatest common divisor. 188 """ 189 190 # see [Davenport], Appendix, p. 214 191 192 if abs(q) < abs(r): 193 p, a, b = _extended_euclidean(r, q) 194 return p, b, a 195 196 Q = 1, 0 197 R = 0, 1 198 199 while r: 200 quot, t = divmod(q, r) 201 T = Q[0] - quot*R[0], Q[1] - quot*R[1] 202 q, r = r, t 203 Q, R = R, T 204 205 return q, Q[0], Q[1]
206
207 208 209 210 -def _gcd(q, r):
211 return _extended_euclidean(q, r)[0]
212
213 214 215 216 -def _simplify_fraction((a, b)):
217 gcd = _gcd(a,b) 218 return (a//gcd, b//gcd)
219
220 221 222 223 -class SimplexCubature(object):
224 """Cubature on an M{n}-simplex. 225 226 cf. 227 A. Grundmann and H.M. Moeller, 228 Invariant integration formulas for the n-simplex by combinatorial methods, 229 SIAM J. Numer. Anal. 15 (1978), 282--290. 230 231 This cubature rule has both negative and positive weights. 232 """
233 - def __init__(self, order, dimension):
234 s = order 235 n = dimension 236 d = 2*s+1 237 238 from pytools import \ 239 generate_decreasing_nonnegative_tuples_summing_to, \ 240 generate_unique_permutations, \ 241 factorial, \ 242 wandering_element 243 244 points_to_weights = {} 245 246 for i in xrange(s+1): 247 weight = (-1)**i * 2**(-2*s) \ 248 * (d + n-2*i)**d \ 249 / factorial(i) \ 250 / factorial(d+n-i) 251 252 for t in generate_decreasing_nonnegative_tuples_summing_to(s-i, n+1): 253 for beta in generate_unique_permutations(t): 254 denominator = d+n-2*i 255 point = tuple( 256 _simplify_fraction((2*beta_i+1, denominator)) 257 for beta_i in beta) 258 259 points_to_weights[point] = points_to_weights.get(point, 0) + weight 260 261 from operator import add 262 263 vertices = [-1 * numpy.ones((n,))] \ 264 + [numpy.array(x) for x in wandering_element(n, landscape=-1, wanderer=1)] 265 266 self.pos_points = [] 267 self.pos_weights = [] 268 self.neg_points = [] 269 self.neg_weights = [] 270 271 dim_factor = 2**n 272 for p, w in points_to_weights.iteritems(): 273 real_p = reduce(add, (a/b*v for (a,b),v in zip(p, vertices))) 274 if w > 0: 275 self.pos_points.append(real_p) 276 self.pos_weights.append(dim_factor*w) 277 else: 278 self.neg_points.append(real_p) 279 self.neg_weights.append(dim_factor*w) 280 281 self.points = self.pos_points + self.neg_points 282 self.weights = self.pos_weights + self.neg_weights 283 284 self.pos_info = zip(self.pos_points, self.pos_weights) 285 self.neg_info = zip(self.neg_points, self.neg_weights)
286 287
288 - def __call__(self, f):
289 return sum(f(x)*w for x, w in self.pos_info) \ 290 + sum(f(x)*w for x, w in self.neg_info)
291 292 293 294 295 if __name__ == "__main__": 296 cub = SimplexCubature(2, 3) 297 print cub.weights 298 for p in cub.points: 299 print p 300