Template for the composite high-order Legendre discretization toolkit for HW4 of the 2013 integral equations class at UIUC

legendre_discr.py

from __future__ import division
import numpy as np
import numpy.linalg as la
import scipy.special as sp


class CompositeLegendreDiscretization:
    """A discrete function space on a 1D domain consisting of multiple
    subintervals, each of which is discretized as a polynomial space of
    maximum degree *order*. (see constructor arguments)

    There are :attr:`nintervals` * :attr:`npoints` degrees of freedom
    representing a function on the domain. On each subinterval, the
    function is represented by its function values at Gauss-Legendre
    nodes (see :func:`scipy.speical.legendre`) mapped into that
    subinterval.

    .. note::

        While the external interface of this discretization
        is exclusively in terms of vectors, it may be practical
        to internally reshape (see :func:`numpy.reshape`) these
        vectors into 2D arrays of shape *(nintervals, npoints)*.

    The object has the following attributes:

    .. attribute:: intervals

        The *intervals* constructor argument.

    .. attribute:: nintervals

        Number of subintervals in the discretization.

    .. attribute:: npoints

        Number of points on each interval. Equals *order+1*.

    .. attributes:: nodes

        A vector of ``(nintervals*npoints)`` node locations, consisting of
        Gauss-Legendre nodes that are linearly (or, to be technically correct,
        affinely) mapped into each subinterval.

    """

    def __init__(self, intervals, order):
        """
        :arg intervals: determines the boundaries of the subintervals to
            be used. If this is ``[a,b,c]``, then there are two subintervals
            :math:`(a,b)` and :math:`(b,c)`.
            (and the overall domain is :math:`(a,c)`)

        :arg order: highest polynomial degree being used
        """
        # add code here

    def integral(self, f):
        r"""Use Gauss-Legendre quadrature on each subinterval to approximate
        and return the value of

        .. math::

            \int_a^b f(x) dx

        where :math:`a` and :math:`b` are the left-hand and right-hand edges of
        the computational domain.

        :arg f: the function to be integrated, given as function values
            at :attr:`nodes`
        """

        # add code here

    def left_indefinite_integral(self, f):
        r"""Use a spectral integration matrix on each subinterval to
        approximate and return the value of

        .. math::

            g(x) = \int_a^x f(x) dx

        at :attr:`nodes`, where :math:`a` is the left-hand edge of the
        computational domain.

        The computational cost of this routine is linear in
        the number of degrees of freedom.

        :arg f: the function to be integrated, given as function values
            at :attr:`nodes`
        """

        # add code here

    def right_indefinite_integral(self, f):
        r"""Use a spectral integration matrix on each subinterval to
        approximate and return the value of

        .. math::

            g(x) = \int_x^b f(x) dx

        at :attr:`nodes`, where :math:`b` is the left-hand edge of the
        computational domain.

        The computational cost of this routine is linear in
        the number of degrees of freedom.

        :arg f: the function to be integrated, given as function values
            at :attr:`nodes`
        """

        # add code here