lab5.cpp

using namespace std;

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

double f(double x) {
    return 10.0 * asin(x) * cos(11.0 * x) * cos(x);
}

vector<double> single_division(double **A, int N) {
    int i, j, k, l;
    double tmp1, tmp2, tmp3;
    vector<double> res(N);

    for (int j = 0; j < N; j++) {
        tmp1 = A[j][j];
        for (int i = j; i < N + 1; i++) {
            A[j][i] /= tmp1;
        }
        for (int k = j + 1; k < N; k++) {
            tmp2 = A[k][j];
            for (l = j; l < N + 1; l++) {
                A[k][l] -= A[j][l] * tmp2;
            }
        }
    }
    for (int j = N - 1; j >= 0; j--) {
        tmp3 = 0;
        for (int i = N - 1; i > j; i--) {
            tmp3 += A[j][i] * res[i];
        }
        res[j] = A[j][N] - tmp3;
    }
    return res;
}

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, int k1, int k2, double a, double b, int N){
    double t = (2 * x - (a + b)) / (b - a);
    return chebyshev(t, k1) * (k2 == N ? f(x): chebyshev(t, k2));
}

double trapezoid(double a, double b, int k1, int k2, int n, int N){
    double h = (b - a) / n, s = (func(a, k1, k2, a, b, N) + func(b, k1, k2, a, b, N)) * 0.5;
    for (double x = a + h; x < b; x += h) {
        s += func(x, k1, k2, a, b, N);
    }
    return s * h;
}

double integral(double a, double b, int k1, int k2, int N){
    double error, In, I2n, eps = 1e-3;
    int n = 2 * ceil((b - a) / sqrt(eps) * 0.5);
    In = trapezoid(a, b, k1, k2, n, N);
    I2n = trapezoid(a, b, k1, k2, 2 * n, N);
    error = fabs(In - I2n) / 3;

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

// double Runge(int k1, int k2, double a, double b, int N) {
//  double error = 1e-8;

//  int n = 2 * ceil((b - a) / pow(error, 0.25) / 2.0);
//  double h = (b - a) / n;

//  double R = 0, In, I2n = Simpson(k1, k2, a, b, n, N);
//  do {
//      In = I2n;
//      I2n = Simpson(k1, k2, a, b, n *= 2, N);
//      R = fabs(In - I2n) / 15;
//  } while (R > error);

//  return I2n;
// }

vector<double> Ak(double a, double b, int N){
    double ** A = new double*[N];
    for (int i = 0; i < N + 1; i++) {
        A[i] = new double[N + 1];
    }

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

    return single_division(A, N);
}

double P(double x, double a, double b, int n, vector<double> A) {
    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 N_P(double eps, double a, double b){
    double y, x, error, step = 0.2;
    int n = 2, k = (int)(b - a) / step;

    do {
        n *= 2;
        error = 0;
        x = a;
        for (int i = 0; i <= k; ++i) {
            y = f(x) - P(x, a, b, n, Ak(a, b, n));
            error += y * y;
            x += step;
        }
        error = sqrt(error / (k + 1));
    } while (error > eps);

    return n;
}


int main(){
    double a = -1, b = 1;
    double eps = 1e-2;

    int N = N_P(eps, a, b) * 2;
    cout << "N = " << N << endl;

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

    }
}