Методы вычислений 10.11.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;

mt19937 gen(time(NULL));
uniform_real_distribution<double> interval_double(-10, 10);

const double v = 10;

vec operator * (const mat& A, const vec& x)
{
    int n = x.size();

    vec b(n);

    for (int i = 0; i < n; ++i)
        for (int j = 0; j < n; ++j)
            b[i] += A[i][j] * x[j];

    return b;
}

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

    vec p, q;

    // добавляем коэффициенты P2 и Q2
    // a_i = A[i - 1][i - 2]
    // b_i = -A[i - 1][i - 1]
    // c_i = A[i - 1][i]
    // d_i = b[i - 1]
    p.push_back(-A[0][1] / A[0][0]);
    q.push_back(b[0] / A[0][0]);

    for (int i = 2; i < n; ++i)
    {
        // добавляем коэффициенты P_(i + 1) и Q_(i + 1)
        q.push_back((A[i - 1][i - 2] * q.back() - b[i - 1]) / (-A[i - 1][i - 1] - A[i - 1][i - 2] * p.back()));
        p.push_back(A[i - 1][i] / (-A[i - 1][i - 1] - A[i - 1][i - 2] * p.back()));
    }

    vec res;

    // добавляем коэффициент Q_(n + 1) и сразу записываем его в решение
    res.push_back((A[n - 1][n - 2] * q.back() - b[n - 1]) / (-A[n - 1][n - 1] - A[n - 1][n - 2] * p.back()));

    for (int i = n - 2; i >= 0; --i)
        res.push_back(p[i] * res.back() + q[i]);

    reverse(res.begin(), res.end());

    return res;
}

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

    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";
}

mat task_matrix(int n)
{
    mat res(n, vec(n));

    for (int i = 0; i < n; ++i)
        for (int j = 0; j < n; ++j)
        {
            if (i == j) res[i][j] = v + i;
            else if (0 < abs(i - j) && abs(i - j) < 2) res[i][j] = (v + i) / 100.0;
        }

    return res;
}

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

    int n = 10;
    mat A = task_matrix(n);

    print(A);

    vec x(n);
    for (int i = 0; i < n; ++i)
        x[i] = interval_double(gen);

    vec b = A * x;

    vec checkx = Tridiagonal(A, b);

    print(x);
    print(checkx);

    return 0;
}