midpoint_value.py

'''
This script investigates midpoint value of nodal basis functions of Lagrange
finite element on [0, 1].
'''

from numpy import zeros, linspace
from numpy.linalg import inv
import matplotlib.pyplot as plt


def get_value(q, midpoint=1./2):
    # Points where polynomials are evaluated. The degrees of freedom are point
    # evaluations in points
    points = [i/float(q) for i in range(q+1)]

    # Monomial basis functions of polynomials of order q
    m_basis = [lambda x, i=i: x**i for i in range(q + 1)]

    # Matrix = point evaluation at point for all basis
    B = zeros((q+1, q+1))
    for i, point in enumerate(points):
        for j, base_f in enumerate(m_basis):
            B[i, j] = base_f(point)
    print 'B'
    print B
    print

    # Compute the coefficient matrix to get nodal basis
    A = inv(B).T
    print 'A'
    print A
    print

    # Evaluate the basis functions
    midpoint = 1./2
    for i in range(q+1):
        coefs = A[i, :]
        value = sum(coef*basis(midpoint) for coef, basis in zip(coefs, m_basis))
        print '%d-th nodal basis value at x=%g is %g' % (i, midpoint, value)

    # Plot
    fig = plt.figure()
    ax = plt.subplot(111)
    x = linspace(0, 1, 20*q)
    for i in range(q+1):
        coefs = A[i, :]
        y = zeros(len(x))
        for coef, basis in zip(coefs, m_basis):
            y += coef*basis(x)
        plt.plot(x, y, label=r'$\phi_{%d}$' % i)
    ax.legend(loc='center left', bbox_to_anchor=(1, 0.5), fancybox=True)
    plt.xlabel(r'$x$')
    plt.ylabel(r'$y$')
    plt.title('Nodal basis functions of $CG_{%d}$' % i)
    plt.show()


if __name__ == '__main__':
    import sys
    q = int(sys.argv[1])
    get_value(q)