lab5.cpp

using namespace std;

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

double f(double x) {
    return 1 / x - 0.1 * x * x * sin(2 * x);
}

double *single_division(double *A, int N) {
    double *x = new double[N];

    for (int j = 0; j < N; j++) {
        double m = 1.0 / A[j * (N + 1) + j];
        for (int i = j; i < N + 1; i++) {
            A[j * (N + 1) + i] *= m;
        }
        for (int k = j + 1; k < N; k++) {
            double m = A[k * (N + 1) + j];
            for (int l = j; l < N + 1; l++) {
                A[k * (N + 1) + l] -= A[j * (N + 1) + l] * m;
            }
        }
    }
    for (int j = N - 1; j >= 0; j--) {
        double s = 0;
        for (int i = N - 1; i > j; i--) {
            s += A[j * (N + 1) + i] * x[i];
        }
        x[j] = A[j * (N + 1) + N] - s;
    }
    return x;
}

double chebyshev(double x, int n) {
    double Tn = x, Tn_minus_1 = 1;
    if (n == 0) return Tn_minus_1;

    for (int i = 1; i < n; i++) {
        double T = 2 * x * Tn - Tn_minus_1;

        Tn_minus_1 = Tn;
        Tn = T;
    }

    return Tn;
}

double func(double x, double a, double b, int k1, int k2, int N) {
    double t = (2 * x - (a + b)) / (b - a);
    return chebyshev(t, k1) * (k2 == N ? f(x): chebyshev(t, k2));
}

double simpson(double a, double b, int k1, int k2, int n, int N) {
    double s1 = 0, s2 = 0, h = (b - a) / n;
    for (int i = 1; i < n; i++) {
        if (i % 2 == 0) s2 += func(a + i * h, a, b, k1, k2, N);
        else s1 += func(a + i * h, a, b, k1, k2, N);
    }

    return h / 3 * (2 * s2 + 4 * s1 + func(a, a, b, k1, k2, N) + func(b, a, b, k1, k2, N));
}

double integral(double a, double b, int k1, int k2, int N) {
    int n = (int((b - a) / sqrt(1e-3)) >> 1) << 1;

    double I2n = simpson(a, b, k1, k2, n, N), In = 0;
    double error = fabs(In - I2n) / 3;

    while (error > 1e-3) {
        In = I2n;
        I2n = simpson(a, b, k1, k2, n *= 2, N);
        error = fabs(In - I2n) / 3;
    }
    return I2n;
}

double *Ak(double a, double b, int N) {
    double A[(N + 1) * (N + 1)];

    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N + 1; j++) {
            A[i * (N + 1) + j] = A[j * (N + 1) + i] = integral(a, b, i, j, N);
        }
        A[i * (N + 1) + N] = integral(a, b, i, N, N);
    }

    return single_division(A, N);
}

double P(double x, double a, double b, double *A, int N) {
    double t = (2 * x - b - a) / (b - a);
    double s = 0;
    for (int i = 0; i < N; i++) {
        s += A[i] * chebyshev(t, i);
    }
    return s;
}

int main() {
    double a = 2, b = 11;
    double eps = 1e-2;
    double step = 0.1;
    int k = k = (int)(b - a) / step;
    int N = 2;

    for(;;) {
        N *= 2;
        double error = 0;
        double x = a;
        for (int i = 0; i <= k; ++i) {
            auto A = Ak(a, b, N);
            double delta = f(x) - P(x, a, b, A, N);
            error += delta * delta;
            x += step;
            delete[] A;
        }
        if (sqrt(error / (k + 1)) <= eps) break;
    };

    cout << "N = " << N << endl;

    auto A = Ak(a, b, N);
    for (double x = a; x <= b; x += step) {
        printf("%-5.1f  %.8f\n", x, P(x, a, b, A, N));
    }
    delete[] A;
}