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
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
|
# Copyright 2017 Christoph Groth (CEA).
#
# This file is part of Vquad. It is subject to the license terms in the file
# LICENSE.rst found in the top-level directory of this distribution.
from fractions import Fraction as Frac
from collections import defaultdict
import numpy as np
from scipy.linalg import norm, inv
__all__ = ['sizes', 'nodes', 'newton_coeffs', 'inv_Vs', 'V_cond_nums', 'Ts',
'alpha', 'gamma']
def legendre(n):
"""Return the first n Legendre polynomials.
The polynomials have *standard* normalization, i.e.
int_{-1}^1 dx L_n(x) L_m(x) = delta(m, n) * 2 / (2 * n + 1).
The return value is a list of list of fraction.Fraction instances.
"""
result = [[Frac(1)], [Frac(0), Frac(1)]]
if n <= 2:
return result[:n]
for i in range(2, n):
# Use Bonnet's recursion formula.
new = (i + 1) * [Frac(0)]
new[1:] = (r * (2*i - 1) for r in result[-1])
new[:-2] = (n - r * (i - 1) for n, r in zip(new[:-2], result[-2]))
new[:] = (n / i for n in new)
result.append(new)
return result
def newton(n):
"""Compute the monomial coefficients of the Newton polynomial over the
nodes of the n-point Clenshaw-Curtis quadrature rule.
"""
# The nodes of the Clenshaw-Curtis rule are x_i = -cos(i * Pi / (n-1)).
# Here, we calculate the coefficients c_i such that sum_i c_i * x^i
# = prod_i (x - x_i). The coefficients are thus sums of products of
# cosines.
#
# This routine uses the relation
# cos(a) cos(b) = (cos(a + b) + cos(a - b)) / 2
# to efficiently calculate the coefficients.
#
# The dictionary 'terms' descibes the terms that make up the
# monomial coefficients. Each item ((d, a), m) corresponds to a
# term m * cos(a * Pi / n) to be added to prefactor of the
# monomial x^(n-d).
mod = 2 * (n-1)
terms = defaultdict(int)
terms[0, 0] += 1
for i in range(n):
newterms = []
for (d, a), m in terms.items():
for b in [i, -i]:
# In order to reduce the number of terms, cosine
# arguments are mapped back to the inteval [0, pi/2).
arg = (a + b) % mod
if arg > n-1:
arg = mod - arg
if arg >= n // 2:
if n % 2 and arg == n // 2:
# Zero term: ignore
continue
newterms.append((d + 1, n - 1 - arg, -m))
else:
newterms.append((d + 1, arg, m))
for d, s, m in newterms:
terms[d, s] += m
c = (n + 1) * [0]
for (d, a), m in terms.items():
if m and a != 0:
raise ValueError("Newton polynomial cannot be represented exactly.")
c[n - d] += m
# The check could be removed and the above line replaced by
# the following, but then the result would be no longer exact.
# c[n - d] += m * np.cos(a * np.pi / (n - 1))
cf = np.array(c, float)
assert all(int(cfe) == ce for cfe, ce in zip(cf, c)), 'Precision loss'
cf /= 2.**np.arange(n, -1, -1)
return cf
def scalar_product(a, b):
"""Compute the polynomial scalar product int_-1^1 dx a(x) b(x).
The args must be sequences of polynomial coefficients. This
function is careful to use the input data type for calculations.
"""
la = len(a)
lc = len(b) + la + 1
# Compute the even coefficients of the product of a and b.
c = lc * [a[0].__class__()]
for i, bi in enumerate(b):
if bi == 0:
continue
for j in range(i % 2, la, 2):
c[i + j] += a[j] * bi
# Calculate the definite integral from -1 to 1.
return 2 * sum(c[i] / (i + 1) for i in range(0, lc, 2))
def newton_legendre(sizes):
"""Calculate the decompositions of Newton polynomials (over the nodes
of the n-point Clenshaw-Curtis quadrature rule) in terms of
Legandre polynomials.
The parameter 'sizes' is a sequence of numers of points of the
quadrature rule. The return value is a corresponding sequence of
normalized Legendre polynomial coefficients.
"""
legs = legendre(max(sizes) + 1)
result = []
for n in sizes:
poly = []
a = list(map(Frac, newton(n)))
for b in legs[:n + 1]:
igral = scalar_product(a, b)
# Normalize & store. (The polynomials returned by
# legendre() have standard normalization that is not
# orthonormal.)
poly.append(np.sqrt((2*len(b) - 1) / 2) * igral)
result.append(np.array(poly))
return result
def vandermonde(nodes):
V = [np.ones(nodes.shape), nodes.copy()]
for i in range(2, len(nodes)):
V.append((2*i-1) / i * nodes * V[-1] - (i-1) / i * V[-2])
for i in np.arange(len(nodes)):
V[i] *= np.sqrt(i + 0.5)
return np.array(V).T
def precalculate(n_levels=5):
"""Precalculate tables for adaptive quadrature based on the
Clenshaw-Curtis rule and Legendre polynomials.
The Clenshaw-Curtis rule with three points is the lowest rule that
contains the center of the interval, hence we define it as level 0.
"""
global sizes, nodes, newton_coeffs, inv_Vs, V_cond_nums, Ts, alpha, gamma
# Points of the Clenshaw-Curtis rule.
sizes = [2**(level + 1) + 1 for level in range(n_levels)]
nodes = [-np.cos(np.pi / (n - 1) * np.arange(n)) for n in sizes]
# Set central rule points precisely to zero. This does not really
# matter in practice, but is useful for tests.
for l in nodes:
l[len(l) // 2] = 0.0
# Vandermonde-like matrices and their condition numbers
V = list(map(vandermonde, nodes))
inv_Vs = inv_Vs = list(map(inv, V))
V_cond_nums = [norm(a, 2) * norm(b, 2) for a, b in zip(V, inv_Vs)]
# Shift matrices
Ts = [inv_Vs[-1] @ vandermonde((nodes[-1] + a) / 2) for a in [-1, 1]]
# Newton polynomials
newton_coeffs = newton_legendre(sizes)
# Other downdate matrices
k = np.arange(sizes[-1])
alpha = np.sqrt((k+1)**2 / (2*k+1) / (2*k+3))
gamma = np.concatenate([[0, 0],
np.sqrt(k[2:]**2 / (4*k[2:]**2-1))])
precalculate()
|
