lab5

#include <cstdio>
#include <cmath>

double a = 0, b = 1.2;
double eps = 1e-3;
int k = (b - a) * 50;
double step = (b - a) / (double) k;
double sqrError = eps * eps * (k + 1);

double Function(double x) {
    return 10 * x * x * cosh(x) * sin(13 * 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 Polinom(double x, double *ak, int N) {
    double t = (2 * x - b - a) / (b - a);
    double s = 0;
    for (int i = 0; i < N; i++) {
        s += ak[i] * Chebyshev(t, i);
    }
    return s;
}

double F(double x, int k1, int k2, int N) {
    double t = (2 * x - (a + b)) / (b - a);

    if (k2 == N) {
        return Chebyshev(t, k1) * Function(x);
    } else {
        return Chebyshev(t, k1) * Chebyshev(t, k2);
    }
}

double Simpson(int k1, int k2, int n, int N) {
    double s1 = 0;
    double s2 = 0;
    double h = (b - a) / (double) n;

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

        if (i % 2 == 0) {
            s2 += F(x, k1, k2, N);
        } else {
            s1 += F(x, k1, k2, N);
        }
    }

    return h / 3.0 * (2.0 * s2 + 4.0 * s1 + F(a, k1, k2, N) + F(b, k1, k2, N));
}

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

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

    loop:
        double In = I2n;
        I2n = Simpson(k1, k2, n *= 2, N);
        if (fabs(In - I2n) >= 15e-8) goto loop;

    return I2n;
}

void Ak(double *ak, int N) {
    double matrix[N][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 = 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] * ak[j];
        }
        ak[i] = matrix[i][N] - s;
    }
}

int main() {
    int N = 2;

    loop:
        double error = 0;
        double ak[N *= 2];
        Ak(ak, N);
        for (int i = 0; i <= k; i++) {
            double x = a + i * step;
            error += pow(Function(x) - Polinom(x, ak, N), 2);
        }

        if (error > sqrError) goto loop;

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

    for (int i = 0; i <= k; i++) {
        double x = a + i * step;

        printf("%-.2f\t%.8f\n", x, Polinom(x, ak, N));
    }
}