lab5

#include <cstdio>
#include <cmath>
#include <vector>

double a = 1;
double b = 6.5;

double F(double x) {
    return exp(pow(x + 5.0, 0.25)) * sin(3.0 * x) * cos(x);
}

double Chebyshev(double x, int n) {
    if (n == 0) return 1;

    double tn = x;
    double tnm1 = 1;

    for (int i = 1; i < n; i++) {
        double tnp1 = 2 * x * tn - tnm1;

        tnm1 = tn;
        tn = tnp1;
    }

    return tn;
}

double F(double x, 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(int k1, int k2, int n, int N) {
    double h = (b - a) / (double)n;
    double s = F(a, k1, k2, N) + F(b, k1, k2, N);

    for (int i = 1; i < n; i++) {
        double x = a + i * h;

        s += F(x, k1, k2, N) * (i % 2 == 0 ? 2.0 : 4.0);
    }

    return h * s / 3.0;
}

double Runge(int k1, int k2, int N) {
    int n = (b - a) * 100;

    double I2n = Simpson(k1, k2, n, N);

    for (;;) {
        double In = I2n;
        I2n = Simpson(k1, k2, n *= 2, N);
        if (fabs(In - I2n) < 15e-8) break;
    }

    return I2n;
}

double Polinom(double x, const std::vector<double> &p) {
    double t = (2 * x - b - a) / (b - a);
    double s = 0;

    for (int i = 0, N = p.size(); i < N; i++) {
        s += p[i] * Chebyshev(t, i);
    }
    return s;
}

int main() {
    std::vector<double> p(4);

    while (true) {
        int N = p.size();

        double **matrix = new double*[N];
        for (int i = 0; i < N; i++) {
            matrix[i] = new double[N + 1];
        }

        for (int i = 0; i < N; i++) {
            for (int j = 0; j <= N; j++) {
                matrix[i][j] = Runge(i, j, N);
            }
        }

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

        for (int i = N - 1; i >= 0; i--) {
            double s = 0;
            for (int j = i + 1; j < N; j++) {
                s += matrix[i][j] * p[j];
            }
            p[i] = matrix[i][N] - s;
        }

        for (int i = 0; i < N; i++) {
            delete[] matrix[i];
        }
        delete[] matrix;

        double err = 0;
        for (int i = 0; i <= 60; i++) {
            double x = a + i * 0.02;
            double d = F(x) - Polinom(x, p);
            err += d * d;
        }

        if (err < 61e-6) break;
        p.resize(N << 1);
    }

    printf("N = %i\n", p.size());

    for (int i = 0; i <= 60; i++) {
        double x = a + i * 0.02;

        printf("%-.2f\t%.8f\n", x, Polinom(x, p));
    }
    return 0;
}