Методы вычислений 20.10.22

#define _USE_MATH_DEFINES
#include <iostream>
#include <iomanip>
#include <vector>
#include <random>
#include <chrono>
#include <cmath>
#include <ctime>
using namespace std;
using namespace chrono;

typedef vector <double> vec;
typedef vector <vec> mat;

struct Polynomial
{
    double a, b, c, d;

    Polynomial() {}
    Polynomial(double a, double b, double c, double d) : a(a), b(b), c(c), d(d) {}
};

const double v = 10;

double P(double x)
{
    return (-x * x * x / 6.0) + (3 * x * x / 2.0) - (x / 3.0) + 10.0;
}

vec Gauss(mat& A, vec& b)
{
    int n = b.size();

    // прямой ход
    for (int i = 0; i < n; ++i)
    {
        if (fabs(A[i][i]) < 1e-5)
        {
            int m = i;

            for (int k = i + 1; k < n; ++k)
                if (fabs(A[k][i]) > fabs(A[m][i])) m = k;

            swap(b[i], b[m]);

            for (int j = i; j < n; ++j)
                swap(A[i][j], A[m][j]);
        }

        // нормируем строку
        for (int j = i + 1; j < n; ++j)
            A[i][j] /= A[i][i];

        // нормируем строку
        b[i] /= A[i][i]; A[i][i] = 1;

        // вычитаем из нижних строк текущую строку
        // получаем в i-ом столбце нули
        for (int k = i + 1; k < n; ++k)
        {
            for (int j = i + 1; j < n; ++j)
                A[k][j] -= A[i][j] * A[k][i];

            b[k] -= b[i] * A[k][i];
        }
    }

    vec x(n);

    // обратный ход
    for (int i = n - 1; i >= 0; --i)
    {
        x[i] = b[i];

        for (int j = i + 1; j < n; ++j)
            x[i] -= A[i][j] * x[j];
    }

    return x;
}

void print(mat& A)
{
    int n = A.size();

    cout << "Матрица:\n";

    for (int i = 0; i < n; ++i, cout << "\n")
        for (int j = 0; j < n; ++j)
            cout << left << setw(10) << A[i][j];
}

void print(vec& a)
{
    int n = a.size();

    cout << "Вектор: (";

    for (int i = 0; i < n - 1; ++i)
        cout << a[i] << ", ";

    cout << a[n - 1] << ")\n\n";
}

void print(mat& A, vec& b)
{
    int n = b.size();

    cout << "Матрица СЛУ имеет вид:\n";

    for (int i = 0; i < n; ++i, cout << "\n")
    {
        for (int j = 0; j < n; ++j)
            cout << left << setw(8) << A[i][j];

        cout << "| " << b[i];
    }

    cout << "\n";
}

double solve(int n, mat& A, vec& b, vec& x)
{
    vec checkx = Gauss(A, b);

    double res = 0.0;

    for (int i = 0; i < n; ++i)
        res = max(res, fabs(x[i] - checkx[i]));

    return res;
}

int main()
{
    srand(time(NULL));
    setlocale(LC_ALL, "Russian");
    cout << fixed << setprecision(2);

    int n = 3;
    double h = 1;

    vec x = { 0, 1, 2, 3 };
    vec y = { 10, 11, 14, 18 };

    mat A(4 * n, vec(4 * n, 0));
    vec b(4 * n, 0);

    // a[k] = x[4k]
    // b[k] = x[4k + 1]
    // c[k] = x[4k + 2]
    // d[k] = x[4k + 3]
    // k = 0 ... n - 1

    // Уравнения (1) : 1 * a[i] = y[i]
    for (int i = 0, j = 0; i < n; ++i, j += 4)
    {
        A[i][j] = 1;
        b[i] = y[i];
    }

    // Уравнения (2) : 1 * a[i - n] + h * b[i - n] + h^2 * c[i - n] + h^3 * d[i - n] = f[i - n + 1]
    for (int i = n, j = 0; i < 2 * n; ++i, j += 4)
    {
        A[i][j] = 1;
        A[i][j + 1] = h;
        A[i][j + 2] = h * h;
        A[i][j + 3] = h * h * h;

        b[i] = y[i - n + 1];
    }

    // Уравнения (3) : 1 * b[i - 2n] + 2h * c[i - 2n] + 3h^2 * d[i - 2n] - b[i - 2n + 1] = 0
    for (int i = 2 * n, j = 0; i < 3 * n - 1; ++i, j += 4)
    {
        A[i][j + 1] = 1;
        A[i][j + 2] = 2 * h;
        A[i][j + 3] = 3 * h * h;
        A[i][j + 5] = -1;
    }

    // Уравнения (4) : 1 * c[i - 3n + 1] + 3h * d[i - 3n + 1] - c[i - 3n + 2] = 0
    for (int i = 3 * n - 1, j = 0; i < 4 * n - 2; ++i, j += 4)
    {
        A[i][j + 2] = 1;
        A[i][j + 3] = 3 * h;
        A[i][j + 6] = -1;
    }

    // Уравнения (5)
    A[4 * n - 2][2] = 1;
    A[4 * n - 1][4 * n - 2] = 1;
    A[4 * n - 1][4 * n - 1] = 3 * h;

    print(A, b);

    vec sol = Gauss(A, b);

    vector <Polynomial> p(n);

    for (int i = 0; i < 4 * n; i += 4)
    {
        p[i / 4].a = sol[i];
        p[i / 4].b = sol[i + 1];
        p[i / 4].c = sol[i + 2];
        p[i / 4].d = sol[i + 3];
    }

    for (int i = 0; i < n; ++i)
    {
        cout << "S[" << i << "](x) : ";
        cout << left << setw(10) << p[i].a << setw(10) << p[i].b << setw(10) << p[i].c << setw(10) << p[i].d << "\n";
    }

    cout << '\n';

    x = { 0, 0.5, 1, 1.5, 2, 2.5, 3 };

    cout << left << "x: ";

    for (int i = 0; i < x.size(); ++i)
        cout << setw(10) << x[i];

    cout << "\n" << left << "y: ";

    for (int i = 0; i < p.size(); ++i)
    {
        Polynomial cur = p[i];

        for (int j = 0; j < 2; ++j)
            cout << setw(10) << cur.a + cur.b * (x[2 * i + j] - i) + cur.c * pow(x[2 * i + j] - i, 2) + cur.d * pow(x[2 * i + j] - i, 3);
    }

    Polynomial last = p[p.size() - 1];
    double last_x = x[x.size() - 1];

    cout << last.a + last.b * (last_x - 2) + last.c * pow(last_x - 2, 2) + last.d * pow(last_x - 2, 3);

    return 0;
}