#include <iostream>
#include <cmath>
#include <vector>

using namespace std;

// коэффициенты для прогоночной матрицы
// tma (tridiagonal matrix algorithm) - метод прогонки
struct tma{
    double a;
    double b;
    double c;
    double f;
};

// вспомогательные прогоночные коэффициенты
// forward sweep - прямая прогонка
struct sweep{
    double X;
    double Z;
};

// приближённое решение в узлах сетки
struct dots{
    double x;
    double y;
};

// граничные условия
struct boundary{
    double alpha_0;
    double alpha_1;
    double beta_0;
    double beta_1;
    double gamma_0;
    double gamma_1;
};

// вспомогательные коэффициенты
struct additional{
    double a;
    double d;
    double phi;
};

// функция при первой производной
double k(double x)
{
    return 1.0;
}

// функция при "нулевой" производной
double q(double x)
{
    return -4.0;
}

// неоднородность в правой части уравнения
double f(double x)
{
    return -cos(2.0 * x);
}

// инициализация сетки
void initGrid(vector<dots> &grid, double a, double h, size_t N)
{
    dots temp;

    for (size_t i = 0; i <= N; i++)
    {
        temp.x = a + i * h;
        temp.y = 0.0;
        grid.push_back(temp);
    }
}

// вычисление вспомогательных коэффициентов для метода прямоугольников
void fillRectangle(vector<additional> &rectangle, const vector<dots> &grid, double h, size_t N)
{
    additional temp;
    for (size_t i = 0; i < N; i++)
    {
        temp.a = k(grid[i + 1].x - 0.5 * h);
        temp.d = q(grid[i + 1].x);
        temp.phi = f(grid[i + 1].x);

        rectangle.push_back(temp);
    }
}

// вычисление вспомогательных коэффициентов для метода трапеций
void fillTrapezium(vector<additional> &trapezium, const vector<dots> &grid, double h, size_t N)
{
    additional temp;
    for (size_t i = 0; i < N; i++)
    {
        temp.a = 2.0 * k(grid[i + 1].x) * k(grid[i].x) / (k(grid[i + 1].x) + k(grid[i].x));
        temp.d = (q(grid[i + 1].x - h / 2.0) + q(grid[i + 1].x + h / 2.0)) / 2.0;
        temp.phi = (f(grid[i + 1].x - h / 2.0) + f(grid[i + 1].x + h / 2.0));

        trapezium.push_back(temp);
    }
}

// заполняется матрица коэффициентов
void fillMatrix(vector<tma> &matrix, const vector<additional> &method, double h, size_t N)
{
    tma temp;

    for (size_t i = 0; i < N - 1; i++)                                  // Из N (количества точек разбиения) вычитается 1,
    {                                                                   // так как число строк в прогоночной матрице на 1 меньше N
        temp.a = method[i + 1].a / pow(h, 2.0);                         // и коэффициенты нумеруются с цифры 1.
        temp.b = -(method[i + 1].a + method[i].a) / pow(h, 2.0)         // Коэффициенты с индексами 0 и N
                 - method[i].d;                                         // используются при задании граничных условий.
        temp.c = method[i].a / pow(h, 2.0);                             //
        temp.f = method[i].phi;                                         //
        matrix.push_back(temp);
    }
}

// заполняется вектор прогоночных коэффициентов
void fillSweepCoefficients(vector<sweep> &coefficients, const vector<tma> &matrix, boundary conditions)
{
    size_t N = matrix.size();
    sweep temp;
    temp.X = -conditions.beta_0 / conditions.alpha_0;
    temp.Z = conditions.gamma_0 / conditions.alpha_0;
    coefficients.push_back(temp);

    for (size_t i = 0; i < N; i++)
    {
        temp.X = -matrix[i].a / (matrix[i].b + matrix[i].c * coefficients[i].X);
        temp.Z = (matrix[i].f - matrix[i].c * coefficients[i].Z) / (matrix[i].b + matrix[i].c * coefficients[i].X);
        coefficients.push_back(temp);
    }
}

// обратная прогонка, заполнение вектора значений искомой функции
void fillSolution(const vector<sweep> &coefficients, vector<dots> &grid, boundary conditions, double a, double h)
{
    size_t N = coefficients.size();

    double y = (conditions.gamma_1 - conditions.alpha_1 * coefficients.back().Z) / (conditions.beta_1 + conditions.alpha_1 * coefficients.back().X);

    grid.back().y = y;

    for (size_t i = N - 1; i < N; i--)                          // Начинается отсчёт с N-1, так как
    {                                                           // точка с индексом N уже посчитана.
        y = coefficients[i].X * grid[i + 1].y + coefficients[i].Z;
        grid[i].y = y;
    }
}

int main()
{
    double pi = 3.1415926;
    size_t N = 2;                                   // число точек разбиения
    double a = .0, b = pi / 4.0;                    // начало и конец отрезка
    double h = (b - a) / N;                         // шаг разбиения
    double mu_1 = pi / 8.0, mu_2 = 0.0, beta = 1.0;

    vector<dots> grid_rectangle, grid_trapezium;
    vector<additional> rectangle, trapezium;
    vector<tma> matrix_rectangle, matrix_trapezium;
    vector<sweep> coefficients_rectangle, coefficients_trapezium;

    // вычисление коэффициентов с использованием метода прямоугольников
    initGrid(grid_rectangle, a, h, N);
    fillRectangle(rectangle, grid_rectangle, h, N);
    fillMatrix(matrix_rectangle, rectangle, h, N);

    boundary conditions_rectangle;
    conditions_rectangle.alpha_0 = beta + q(a) * h / 2.0 + rectangle[0].a / h;
    conditions_rectangle.alpha_1 = 0.0;
    conditions_rectangle.beta_0 = -rectangle[0].a / h;
    conditions_rectangle.beta_1 = 1.0;
    conditions_rectangle.gamma_0 = mu_1 + f(a) * h / 2.0;
    conditions_rectangle.gamma_1 = mu_2;

    fillSweepCoefficients(coefficients_rectangle, matrix_rectangle, conditions_rectangle);
    fillSolution(coefficients_rectangle, grid_rectangle, conditions_rectangle, a, h);

    // вычисление коэффициентов с использованием метода трапеций
    initGrid(grid_trapezium, a, h, N);
    fillTrapezium(trapezium, grid_trapezium, h, N);
    fillMatrix(matrix_trapezium, trapezium, h, N);

    boundary conditions_trapezium;
    conditions_trapezium.alpha_0 = beta + q(a) * h / 2.0 + trapezium[0].a / h;
    conditions_trapezium.alpha_1 = 0.0;
    conditions_trapezium.beta_0 = -trapezium[0].a / h;
    conditions_trapezium.beta_1 = 1.0;
    conditions_trapezium.gamma_0 = mu_1 + f(a) * h / 2.0;
    conditions_trapezium.gamma_1 = mu_2;

    fillSweepCoefficients(coefficients_trapezium, matrix_trapezium, conditions_trapezium);
    fillSolution(coefficients_trapezium, grid_trapezium, conditions_trapezium, a, h);

    // вывод результатов
    cout.precision(10);
    std::cout.setf(std::ios::fixed, std::ios::floatfield);
    cout << "\t" << "x" << "\t" << "Rectangle method" << "\t" << "Trapezium method" << "\t" << "Analitic" << endl;
    for (size_t i = 0; i <= N; i++)
        cout << grid_rectangle[i].x << "\t" << grid_rectangle[i].y << "\t \t" << grid_trapezium[i].y << "\t \t"
             << sin(2.0 * grid_trapezium[i].x) * (grid_trapezium[i].x / 4.0 - pi / 16.0) << endl;

    system("PAUSE");
    return 0;
}