task6

import numpy as np
from scipy import integrate
from scipy.special import jacobi
from tabulate import tabulate

a = -1.0
b = 1.0

Q1 = np.poly1d([1, 1.5, 0, -14.5])
Q2 = np.poly1d([1, 2, -11])


def f(x):
    return 1. - x / 2.


def points(n):
    return np.linspace(a, b, n)


def colloc_points(n):
    return [np.cos(((2 * i - 1) * np.pi) / (2 * n)) for i in range(1, n + 1)]


def w(n):
    if n == 1:
        return [Q1]
    elif n == 2:
        return [Q1, Q2]
    else:
        jacobi_tail = [np.poly1d([1, 0, -2, 0, 1]) *
                       jacobi(i - 3, 2, 2) for i in range(3, n + 1)]
        return [Q1, Q2] + jacobi_tail


def L_diff(w, x):
    w_2 = np.polyder(w, 2)
    w_1 = np.polyder(w, 1)
    ans = -1 / (3 + np.sin(x)) * w_2(x) - x * w_1(x) + (1 + x) * w(x)
    return ans


def L_colloc(n):
    L = np.zeros((n, n))
    g = np.zeros(n)
    W = w(n)
    X = colloc_points(n)

    for i in range(0, n):
        for j in range(0, n):
            L[i, j] = L_diff(W[j], X[i])
        g[i] = f(X[i])
    return L, g


def L_gal(n):
    L = np.zeros((n, n))
    g = np.zeros(n)
    W = w(n)

    for j in range(0, n):
        for i in range(0, n):
            L[j, i] = integrate.quad(lambda x: W[j](x) * L_diff(W[i], x), a, b)[0]
        g[j] = integrate.quad(lambda x: W[j](x) * f(x), a, b)[0]
    return L, g


def coeff(l, n):
    L, g = l(n)
    return np.linalg.solve(L, g)


def solution(W, C, n):
    return lambda x: np.dot(np.array([W[i](x) for i in range(0, n)]), C)


def print_table(n, a, b, l):
    print((' ' * 6) + "|" + (' ' * 14) + 'y_n(x)' + (' ' * 16))
    table = []
    for i in range(3, n + 1):
        W = w(i)
        C = coeff(l, i)
        sol = solution(W, C, i)
        table.append([i, sol(a), sol((a + b) / 2), sol(b)])

    head = ['x = {0}'.format(a), 'x = {0}'.format((a + b) / 2), 'x = {0}'.format(b)]
    print(tabulate(table, headers=['n'] + head * 2, tablefmt="grid", floatfmt="f"))

print('Galerkin method')
print('')
print_table(7, a, b, L_gal)
print('')

print('Collocation method')
print('')
print_table(7, a, b, L_colloc)