gauss_patterson.py 文件源码

def _get_weights(pts):
    '''Given a number of points in [-1, 1], according to

        On some Gauss and Lobatto based integration formulae,
        T. N. L. Patterson,
        Math. Comp. 22 (1968), 877-881,

    one can compute the corresponding weights. One reads there:

    > Thus the weights of an n-point integration formula [...] are given by
    >
    >     omega_i = int_{-1}^{1} L_i(x) dx,
    >
    > (where L_i is the Lagrange polynomial for the point x_i).
    > These weights can be evaluated exactly in a numerically stable fashion
    > using a Gauss formula with n/2 points when n is even and (n + 1)/2 points
    > when n is odd.
    '''
    n = len(pts)

    # Unnormalized Lagrange polynomial: Degree n, 0 at all x_j except x_i.
    def L(i, x):
        return numpy.prod([(x - pts[j]) for j in range(n) if j != i], axis=0)

    # Gauss-Legendre of order k integrates polynomials of degree 2*k-1 exactly.
    # L has degree n-1, so k needs to be n/2 if n is even, and (n+1)/2 if n is
    # odd.
    k = (n // 2) - 1 if n % 2 == 0 else (n+1) // 2
    return numpy.array([
        integrate(
            lambda x, i=i: L(i, x[0]),
            numpy.array([[-1.0], [1.0]]),
            GaussLegendre(k),
            sumfun=lambda a: numpy.array([math.fsum(a)])
            )[0]
        /
        numpy.prod([(pts[i] - pts[j]) for j in range(n) if j != i])
        for i in range(n)
        ])