Basis_Function.py

# -*- coding: utf-8 -*-
"""
@author: Laura C. Kühle

"""
import numpy as np
from sympy import Symbol, integrate

x = Symbol('x')
xi = Symbol('z')


class Vector(object):
    def __init__(self, polynomial_degree):
        self._polynomial_degree = polynomial_degree
        self._basis = self._build_basis_vector(x)
        self._wavelet = self._build_wavelet_vector(xi)

    def get_basis_vector(self):
        return self._basis

    def _build_basis_vector(self, eval_point):
        return []

    def get_wavelet_vector(self):
        return self._wavelet

    def _build_wavelet_vector(self, eval_point):
        return []

    def get_basis_projections(self):
        pass

    def get_multiwavelet_projections(self):
        pass


class Legendre(Vector):
    def _build_basis_vector(self, eval_point):
        return self._calculate_legendre_vector(eval_point)

    def _calculate_legendre_vector(self, eval_point):
        vector = []
        for degree in range(self._polynomial_degree+1):
            if degree == 0:
                vector.append(1.0 + 0*eval_point)
            else:
                if degree == 1:
                    vector.append(eval_point)
                else:
                    poly = (2.0*degree - 1)/degree * eval_point * vector[-1] - (degree-1)/degree * vector[-2]
                    vector.append(poly)
        return vector


class OrthonormalLegendre(Legendre):
    def _build_basis_vector(self, eval_point):
        leg_vector = self._calculate_legendre_vector(eval_point)
        return [leg_vector[degree] * np.sqrt(degree+0.5) for degree in range(self._polynomial_degree+1)]

    def _build_wavelet_vector(self, eval_point):
        degree = self._polynomial_degree

        if degree == 0:
            return [np.sqrt(0.5) + eval_point*0]
        if degree == 1:
            return [np.sqrt(1.5) * (-1 + 2*eval_point), np.sqrt(0.5) * (-2 + 3*eval_point)]
        if degree == 2:
            return [1/3 * np.sqrt(0.5) * (1 - 24*eval_point + 30*(eval_point**2)),
                    1/2 * np.sqrt(1.5) * (3 - 16*eval_point + 15*(eval_point**2)),
                    1/3 * np.sqrt(2.5) * (4 - 15*eval_point + 12*(eval_point**2))]
        if degree == 3:
            return [np.sqrt(15/34) * (1 + 4*eval_point - 30*(eval_point**2) + 28*(eval_point**3)),
                    np.sqrt(1/42) * (-4 + 105*eval_point - 300*(eval_point**2) + 210*(eval_point**3)),
                    1/2 * np.sqrt(35/34) * (-5 + 48*eval_point - 105*(eval_point**2) + 64*(eval_point**3)),
                    1/2 * np.sqrt(5/34) * (-16 + 105*eval_point - 192*(eval_point**2) + 105*(eval_point**3))]
        if degree == 4:
            return [np.sqrt(1/186) * (1 + 30*eval_point + 210*(eval_point**2)
                                      - 840*(eval_point**3) + 630*(eval_point**4)),
                    0.5 * np.sqrt(1/38) * (-5 - 144*eval_point + 1155*(eval_point**2)
                                           - 2240*(eval_point**3) + 1260*(eval_point**4)),
                    np.sqrt(35/14694) * (22 - 735*eval_point + 3504*(eval_point**2)
                                         - 5460*(eval_point**3) + 2700*(eval_point**4)),
                    1/8 * np.sqrt(21/38) * (35 - 512*eval_point + 1890*(eval_point**2)
                                            - 2560*(eval_point**3) + 1155*(eval_point**4)),
                    0.5 * np.sqrt(7/158) * (32 - 315*eval_point + 960*(eval_point**2)
                                            - 1155*(eval_point**3) + 480*(eval_point**4))]

        raise ValueError('Invalid value: Alpert\'s wavelet is only available \
                         up to degree 4 for this application')

    def get_basis_projections(self):
        basis_projection_left = self._build_basis_matrix(xi, 0.5 * (xi - 1))
        basis_projection_right = self._build_basis_matrix(xi, 0.5 * (xi + 1))
        return basis_projection_left, basis_projection_right

    def _build_basis_matrix(self, first_param, second_param):
        matrix = []
        for i in range(self._polynomial_degree + 1):
            row = []
            for j in range(self._polynomial_degree + 1):
                entry = integrate(self._basis[i].subs(x, first_param) *
                                  self._basis[j].subs(x, second_param),
                                  (xi, -1, 1))
                row.append(np.float64(entry))
            matrix.append(row)
        return matrix

    def get_multiwavelet_projections(self):
        wavelet_projection_left = self._build_multiwavelet_matrix(xi, -0.5*(xi-1), True)
        wavelet_projection_right = self._build_multiwavelet_matrix(xi, 0.5*(xi+1), False)
        return wavelet_projection_left, wavelet_projection_right

    def _build_multiwavelet_matrix(self, first_param, second_param, is_left_matrix):
        matrix = []
        for i in range(self._polynomial_degree+1):
            row = []
            for j in range(self._polynomial_degree+1):
                entry = integrate(self._basis[i].subs(x, first_param) * self._wavelet[j].subs(xi, second_param),
                                  (xi, -1, 1))
                if is_left_matrix:
                    entry = entry * (-1)**(j + self._polynomial_degree + 1)
                row.append(np.float64(entry))
            matrix.append(row)
        return matrix