[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
6
8  This program is free software: you can redistribute it and/or modify
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(
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
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
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)
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
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
157      """An M{N}th order Gauss quadrature associated with the Legendre polynomials.
158      """
159 -    def __init__(self, N):
161
162
163
164
166      """A quadrature rule on an arbitrary interval M{(a,b)}. """
167
168 -    def __init__(self, quad, left, right):
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
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
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
<!--
expandto(location.href);
// -->

```

 Generated by Epydoc 3.0.1 on Sat Aug 29 14:33:47 2009 http://epydoc.sourceforge.net