-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy path_stroud_1966.py
110 lines (83 loc) · 3.38 KB
/
_stroud_1966.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
import math
import numpy
import sympy
from ..helpers import article, combine, fsd, pm, untangle, z
from ._helpers import SnScheme
source = article(
authors=["A.H. Stroud"],
title="Some Fifth Degree Integration Formulas for Symmetric Regions",
journal="Mathematics of Computation",
volume="20",
number="93",
month="jan",
year="1966",
pages="90-97",
url="https://doi.org/10.1090/S0025-5718-1966-0191094-8",
)
def stroud_1966_a(n, symbolic=False):
sqrt = sympy.sqrt if symbolic else math.sqrt
a = sqrt(2 * (n + 4))
r2 = (n + 4 - a) / (n + 4)
s2 = (n * (n + 4) + 2 * a) / (n ** 2 + 2 * n - 4) / (n + 4)
B1 = 1 / r2 ** 2 / (n + 2) / (n + 4)
B2 = 1 / s2 ** 2 / 2 ** n / (n + 2) / (n + 4)
r = sqrt(r2)
s = sqrt(s2)
data = [(B1, fsd(n, (r, 1))), (B2, pm(n * [s]))]
points, weights = untangle(data)
points = numpy.ascontiguousarray(points.T)
return SnScheme("Stroud 1966a", n, weights, points, 5, source)
def stroud_1966_b(n, symbolic=False):
frac = sympy.Rational if symbolic else lambda a, b: a / b
sqrt = sympy.sqrt if symbolic else math.sqrt
alpha = 0
s = sqrt(frac(n + alpha + 2, (n + 2) * (n + alpha + 4)))
data = []
B0 = 1
for k in range(1, n + 1):
B = frac(
2 ** (k - n) * (n + 2) * (n + alpha) * (n + alpha + 4),
n * (k + 1) * (k + 2) * (n + alpha + 2) ** 2,
)
B0 -= 2 ** (n - k + 1) * B
r = sqrt(frac((k + 2) * (n + alpha + 2), (n + 2) * (n + alpha + 4)))
v = numpy.concatenate(
[
numpy.zeros((2 ** (n - k + 1), k - 1), dtype=int),
pm(numpy.array([r] + (n - k) * [s])),
],
axis=-1,
)
data.append((B, v))
data.append((B0, z(n)))
points, weights = untangle(data)
points = numpy.ascontiguousarray(points.T)
return SnScheme("Stroud 1966b", n, weights, points, 5, source)
def stroud_1966_c(n, symbolic=False):
frac = sympy.Rational if symbolic else lambda a, b: a / b
sqrt = sympy.sqrt if symbolic else math.sqrt
a = sqrt(2 * (n + 2))
r = sqrt((n + 2 + (n - 1) * a) / (n * (n + 4)))
s = sqrt((n + 2 - a) / (n * (n + 4)))
B0 = frac(4, (n + 2) ** 2)
B1 = frac(n + 4, 2 ** n * (n + 2) ** 2)
data = [(B0, z(n)), (B1, combine(((+r, -r), 1), ((+s, -s), (n - 1))))]
points, weights = untangle(data)
points = numpy.ascontiguousarray(points.T)
return SnScheme("Stroud 1966c", n, weights, points, 5, source)
def stroud_1966_d(n, symbolic=False):
frac = sympy.Rational if symbolic else lambda a, b: a / b
sqrt = sympy.sqrt if symbolic else math.sqrt
a = 2 * sqrt(n + 4)
b = sqrt(2 * (n + 1) * (n + 2) * (n + 4))
r = sqrt((n * (n + 4) + a + (n - 1) * b) / (n * (n + 2) * (n + 4)))
s = sqrt((n * (n + 4) + a - b) / (n * (n + 2) * (n + 4)))
t = sqrt((n + 4 - a) / ((n + 2) * (n + 4)))
B = frac(1, 2 ** n * (n + 1))
# The data is given symbolically, and for large n, those are thousands of points and
# weights. Converting them to float takes a long time. A better approach would be to
# be convert r, s, t first, and assemble the data afterwards.
data = [(B, combine(((+r, -r), 1), ((+s, -s), (n - 1)))), (B, pm(n * [t]))]
points, weights = untangle(data)
points = numpy.ascontiguousarray(points.T)
return SnScheme("Stroud 1966d", n, weights, points, 5, source, 1.019e-14)