diff options
Diffstat (limited to 'vquad/tables.py')
| -rw-r--r-- | vquad/tables.py | 183 |
1 files changed, 183 insertions, 0 deletions
diff --git a/vquad/tables.py b/vquad/tables.py new file mode 100644 index 0000000..267ac51 --- /dev/null +++ b/vquad/tables.py @@ -0,0 +1,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() |