Source code for recursivenodes.quadrature

'''Quadrature nodes used in ``recursivenodes``.'''

import numpy as np
from math import lgamma
from numpy.linalg import eigh
from itertools import product
from recursivenodes.polynomials import _jacobi_recurrence


def _jacobi_matrix(n, a=0., b=0.):
    '''return as numpy arrays, the diagonals for the Jacobi recurrence
    relationship'''
    A, B, C = _jacobi_recurrence(list(range(1,n+2)), a, b)
    alpha = -B[0:n] / A[0:n]
    beta = C[1:(n+1)] / (A[0:n] * A[1:(n+1)])
    return (alpha, beta)


def _monicjacobi(n, x, alpha, beta):
    '''evaluate the monic Jacobi polynomial:

    - for the interval [-1,1]
    - with degree n
    - with the diagonals generated by _jacobi_matrix(n,a,b) for the Jacobi
      weight (1+x)^a (1-x)^b
    - at the points x (a numpy array)

    returns a numpy array with the same shape as x
    '''
    pm2 = np.ones(np.shape(x))
    pm1 = (x - alpha[0]) * pm2
    for i in range(1, n):
        pm2 *= -beta[i-1]
        pm2 += (x - alpha[i]) * pm1
        p = pm2
        pm2 = pm1
        pm1 = p
    return pm1


[docs]def gaussjacobi(n, a=0., b=0.): '''The `n`-point Gauss-Jacobi quadrature points and weights for weight function `(1+x)^a (1-x)^b` for the interval [-1,1], which exactly integrates polynomials with degree up to `2n-1`. Args: n (int): The number of points. a (float, optional): the left exponent of the weight function. b (float, optional): the right exponent of the weight function. Returns: (ndarray, ndarray): Each of shape (`n`,), the points (in ascending order) and weights of the quadrature rule. Reference: :cite:`GoWe69` ''' alpha, beta = _jacobi_matrix(n, a, b) sqbeta = (beta[:-1])**0.5 A = np.zeros((n, n)) A[range(n), range(n)] = alpha A[range(n-1), range(1, n)] = sqbeta A[range(1, n), range(n-1)] = sqbeta x, V = eigh(A) if (a == b): x = (x - x[::-1]) / 2. if (n % 2): x[n // 2] = 0. mu0 = 2**(a+b+1)*np.exp(lgamma(a+1) + lgamma(b+1) - lgamma(a+b+2)) w = V[0, :]**2. * mu0 if (a == b): w = (w + w[::-1]) / 2. return (x, w)
[docs]def lobattogaussjacobi(n, a=0., b=0.): '''The `n`-point Lobatto-Gauss-Jacobi quadrature points and weights for weight function `(1+x)^a (1-x)^b` for the interval [-1,1], which includes the endpoints `-1` and `1` and exactly integrates polynomials with degree up to `2n-3`. Args: n (int): The number of points. a (float, optional): the left exponent of the weight function. b (float, optional): the right exponent of the weight function. Returns: (ndarray, ndarray): Each of shape (`n`,), the points (in ascending order) and weights of the quadrature rule. Reference: :cite:`Golu73` ''' mu0 = 2**(a+b+1)*np.exp(lgamma(a+1) + lgamma(b+1) - lgamma(a+b+2)) if (n == 2): return (np.array([-1., 1.]), np.array([mu0/2., mu0/2.])) alpha, beta = _jacobi_matrix(n, a, b) ends = np.array([-1., 1.]) p_nm2 = _monicjacobi(n-2, ends, alpha, beta) p_nm1 = _monicjacobi(n-1, ends, alpha, beta) rhs = np.array([-p_nm1[0], p_nm1[1]]) Ainv = np.array([[p_nm2[1], -p_nm2[0]], [-p_nm1[1], p_nm1[0]]]) Ainv /= (p_nm1[0] * p_nm2[1] - p_nm2[0] * p_nm1[1]) sol = Ainv.dot(rhs) alpha[-1] = sol[0] beta[-2] = sol[1] sqbeta = (beta[:-1])**0.5 A = np.zeros((n, n)) A[range(n), range(n)] = alpha A[range(n-1), range(1, n)] = sqbeta A[range(1, n), range(n-1)] = sqbeta x, V = eigh(A) if (a == b): x = (x - x[::-1]) / 2. if (n % 2): x[n // 2] = 0. w = V[0, :]**2. * mu0 if (a == b): w = (w + w[::-1]) / 2. x[0] = -1. x[-1] = 1. return (x, w)
[docs]def simplexgausslegendre(d, n): '''An `n^d`-point Gaussian quadrature rule for the biunit simplex that exactly integrates polynomials with degree up to `2n-1`. Args: d (int): The spatial dimension. n (int): The points per coordinate direction Returns: (ndarray, ndarray): The points (shape (`n, \\dots, n, d`)) and weights (shape (`n, \\dots, n,`)) of the quadrature rule. ''' if (d == 1): x, w = gaussjacobi(n, 0., 0.) return x.reshape((n, 1)), w xt, wt = simplexgausslegendre(d-1, n) xz, wz = gaussjacobi(n, d-1, 0.) x = np.ndarray((n,)*d+(d,)) w = np.ndarray((n,)*d) for i in product(range(n), repeat=d-1): for j in range(n): x[i+(j,)][0:d-1] = (xt[i][0:d-1] + 1.) * ((1. - xz[j])/2.) - 1. x[i+(j,)][d-1] = xz[j] w[i+(j,)] = wt[i] * wz[j] / 2**(d-1) return x, w