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

const double v = 10, EPS = 1e-9;

double p(double x)
{
    return x * x;
}

double q(double x)
{
    return x;
}

double f(double x)
{
    return 4 * pow(x, 4) - 3 * v * pow(x, 3) + 6 * x - 2 * v;
}

double y(double x)
{
    return x * x * (x - v);
}

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

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

    double h = 0.1, x = h; int n = v / h + 1;

    mat A(n, vec(n)); vec b(n), expected_y(n);

    A[0][0] = 1; A[n - 1][n - 1] = 1;

    for (int i = 1; i < n - 1; ++i, x += h)
    {
        double pi = p(x), qi = q(x), fi = f(x), yi = y(x);

        A[i][i - 1] = (1 / h / h - pi / 2 / h);
        A[i][i] = (-2 / h / h + qi);
        A[i][i + 1] = (1 / h / h + pi / 2 / h);

        b[i] = fi;
        expected_y[i] = yi;
    }

    // print(A, b);
    // print(expected_y);

    vec observed_y = Tridiagonal(A, b);
    // print(observed_y);

    cout << "+------------+------------+------------+\n";
    cout << "|  expected  |  observed  |   delta    |\n";
    cout << "+------------+------------+------------+\n";

    for (int i = 0; i < n; ++i)
        cout << left << "| " << setw(10) << expected_y[i] << " | " << setw(10) << observed_y[i] << " | " << setw(10) << fabs(expected_y[i] - observed_y[i]) << " |\n";

    cout << "+------------+------------+------------+";

    return 0;
}