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
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
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
61 """An abstract quadrature rule."""
63 self.weights = weights
64 self.points = points
65 self.data = zip(points, weights)
66
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
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 """
85
86 @staticmethod
88 """Return (nodes, weights) for an n-th order Gauss quadrature
89 with the Jacobi polynomials of type (alpha, beta).
90 """
91
92
93
94
95
96
97
98
99 from math import sqrt
100
101 apb = alpha+beta
102
103
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
157 """An M{N}th order Gauss quadrature associated with the Legendre polynomials.
158 """
161
181
186 """Return a tuple (p, a, b) such that p = aq + br,
187 where p is the greatest common divisor.
188 """
189
190
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):
212
217 gcd = _gcd(a,b)
218 return (a//gcd, b//gcd)
219
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 """
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
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