#include <iostream>
#include <algorithm>
#include <cmath>
#include <functional>

using namespace std;

function <double(double)> y = [](double x) {return 3.0 / 2.0 * exp(-x) + x / 2.0 - 1.0 / 2.0;};
function <double(double)> dy = [](double x) {return 1.0 / 2.0 - 3.0 / 2.0 * exp(-x);};
function <double(double)> ddy = [](double x) {return 3.0 / 2.0 * exp(-x);};
function <double(double)> dddy = [](double x) {return -3.0 / 2.0 * exp(-x);};
function <double(double)> ddddy = ddy;
function <double(double)> dddddy = dddy;
function <double(double)> ddddddy = ddddy;
function <double(double, double)> f = [](double x, double y) {return -y + x / 2.0;};

struct Point {
    double node;
    double value;
};

struct DistanceFromPointToX
{
    double distance;
    int pointIndex;
};

bool cmp(DistanceFromPointToX a, DistanceFromPointToX b) {
    return a.distance < b.distance;
}

int factorial (int n) {
    int result = 1;
    for (int i = 1; i <= n; ++i)
        result *= i;
    return result;
}

vector <Point> orderPoints(const double x,
                           vector <Point> points) {
    vector <DistanceFromPointToX> distances;
    for (int i = 0; i < points.size(); ++i)
        distances.push_back({fabs(points[i].node - x), i});
    sort(distances.begin(), distances.end(), cmp);
    vector <Point> sortedPoints;
    for (int i = 0; i < distances.size(); ++i)
        sortedPoints.push_back(points[distances[i].pointIndex]);
    return sortedPoints;
}

void getFiniteDifferenceTable(vector <Point> points,
                              const int n,
                              vector <double> *coefficients,
                              double h) {
    for (int i = 0; i <= n; ++i)
        coefficients[0].push_back(points[i].value * h);
    for (int i = 1; i <= n; ++i)
        for (int j = 0; j <= n - i; ++j)
            coefficients[i].push_back(coefficients[i - 1][j + 1] - coefficients[i - 1][j]);
}

double AdamsRule(const double x,
                 vector <Point> points,
                 const int n,
                 const double h,
                 vector <double> ys) {
    vector <double> coefficients[n + 1];
    getFiniteDifferenceTable(points, n, coefficients, h);
    vector <double> etas;
    for (int i = 0; i <= 4; ++i)
        etas.push_back(coefficients[i][n - i]);
    double result = ys[ys.size() - 1] + etas[0] + etas[1] * 1.0/2.0 + etas[2] * 5.0/12.0 + etas[3] * 3.0/8.0 + etas[4] * 251.0/720.0;
    return result;
}

void FirstMeth(const int n, const double x_0, const double h, const double y_0) {
    vector <double> ys;
    ys.clear();
    ys.push_back(y_0);
    for (int i = 1; i <= n; ++i) {
        double x = x_0 + (i - 1) * h;
        double currentY = y(x) + h / 2.0 * f(x, y(x));
        ys.push_back(y(x) + h * f(x + h / 2.0, currentY));
    }

    for (int i = 0; i < ys.size(); ++i) {
        double x = x_0 + i * h;
        cout << "y(" << x << ") = " << ys[i] << ", абсолютная погрешность " << fabs(y(x) - ys[i]) << endl;
    }
    cout << endl << endl << endl;
}

void SecondMeth(const int n, const double x_0, const double h, const double y_0) {
    vector <double> ys, cauchyY;
    ys.clear(), cauchyY.clear();
    ys.push_back(y_0), cauchyY.push_back(y_0);
    for (int i = 1; i <= n; ++i) {
        double x = x_0 + (i - 1) * h;
        ys.push_back(ys[i - 1] + h * f(x - h, ys[i - 1]));
    }
    for (int i = 1; i <= n; ++i) {
        double x = x_0 + (i - 1) * h;
        cauchyY.push_back(cauchyY[i - 1] + h / 2.0 * (f(x - h, cauchyY[i - 1]) + f(x, ys[i])));
    }

    for (int i = 0; i < cauchyY.size(); ++i) {
        double x = x_0 + i * h;
        cout << "y(" << x << ") = " << cauchyY[i] << ", абсолютная погрешность " << fabs(y(x) - cauchyY[i]) << endl;
    }
    cout << endl << endl << endl;
}


int main() {

    setlocale(LC_ALL, "Russian");
    cout << "Численное решение Задачи Коши для обыкновенного дифференциального уравнения первого порядка" << endl;

    cout << "Дифференциальное уравнение: y'(x) = -y(x) + x / 2" << endl;
    cout << "Начальное условие: y(0) = 1" << endl;
    cout << "Шаг: h = 0.1" << endl;
    cout << "Количество узлов: N = 10" << endl << endl;

    cout << "Решение дифференциального уравнения в общем виде: y(x) = 3 / 2 * exp(-x) + x / 2 - 1 / 2." << endl;

    double x_0 = 0.0;
    double y_0 = 1.0;
    double h = 0.1;
    const int n = 10;

    cout << "Точное решение: " << endl;
    vector <double> ys;
    for (int i = -2; i <= n; ++i) {
        if (i <= 2)
            ys.push_back(y(x_0 + i * h));
        cout << "y(" << x_0 + (double)i * h << ") = " << y(x_0 + (double)i * h) << endl;
    }
    cout << endl << endl << endl;

    cout << "Метод разложения в ряд Тейлора: " << endl << endl;

    for (int i = -2; i <= 2; ++i){
        double x = x_0 + i * h;
        double curResult = y(x_0) * pow(x - x_0, 0) / (double)factorial(0) +
                dy(x_0) * pow(x - x_0, 1) / (double)factorial(1) +
                ddy(x_0) * pow(x - x_0, 2) / (double)factorial(2) +
                dddy(x_0) * pow(x - x_0, 3) / (double)factorial(3) +
                ddddy(x_0) * pow(x - x_0, 4) / (double)factorial(4) +
                dddddy(x_0) * pow(x - x_0, 5) / (double)factorial(5) +
                ddddddy(x_0) * pow(x - x_0, 6) / (double)factorial(6);
        cout << "y(" << x << ") = " << curResult << ", абсолютная погрешность " << fabs(y(x) - curResult) << endl;
    }
    cout << endl << endl << endl;


    cout << "Экстраполяционный метод Адамса 4-го порядка: " << endl;
    vector <Point> points;
    for (int i = -2; i <= 2; ++i) {
        double x = x_0 + i * h;
        points.push_back({x, -y(x) + x / 2.0});
    }
    for (int i = 3; i <= n; ++i) {
        double x = x_0 + i * h;
        double result = AdamsRule(x, points, i + 1, h, ys);
        ys.push_back(result);
        points.push_back({x, -result + x / 2.0});
        cout << "y(" << x << ") = " << result << ", абсолютная погрешность " << fabs(y(x) - result) << endl;
    }
    cout << endl << endl << endl;


    cout << "Метод Рунге-Кутта 4-го порядка: " << endl;
    ys.clear();
    ys.push_back(y_0);
    for (int i = 1; i <= n; ++i) {
        double x = x_0 + (i - 1) * h;
        double k1 = h * f(x, y(x));
        double k2 = h * f(x + h / 2.0, y(x) + k1 / 2.0);
        double k3 = h * f(x + h / 2.0, y(x) + k2 / 2.0);
        double k4 = h * f(x + h, y(x) + k3);
        ys.push_back(ys[ys.size() - 1] + (k1 + 2.0 * k2 + 2.0 * k3 + k4) / 6.0);
    }
    for (int i = 0; i < ys.size(); ++i) {
        double x = x_0 + i * h;
        cout << "y(" << x << ") = " << ys[i] << ", абсолютная погрешность " << fabs(y(x) - ys[i]) << endl;
    }
    cout << endl << endl << endl;


    cout << "Методы Эйлера: " << endl;
    cout << "Нулевой метод: " << endl;
    for (int i = 1; i <= n; ++i) {
        double x = x_0 + (i - 1) * h;
        double result = y(x) + h * f(x, y(x));
        cout << "y(" << x + h << ") = " << result << ", абсолютная погрешность " << fabs(y(x) - result) << endl;
    }
    cout << endl << endl << endl;

    cout << "Первый метод: " << endl;
    FirstMeth(n, x_0, h, y_0);


    vector <double> ys1;

    ys1.clear();
    cout << "Второй метод: " << endl;
    SecondMeth(n, x_0, h, y_0);

    return 0;
}