lab5.cpp

#define _USE_MATH_DEFINES
#include <math.h>
#include <iomanip>
#include <iostream>
using namespace std;

double Function(double x) {
    return log(x * x) * sin(x / 2) * exp(pow(x, 1.0/7.0));
}

double L(int n, double x) {
    double Ln_next, Ln = x, Ln_first = 1;

    if (n == 0) return Ln_first;
    if (n == 1) return Ln;

    int i = 1;
    while (i < n) {
        Ln_next = (1.0 * (2 * i + 1) / (i + 1)) * x * Ln - (1.0 * i / (i + 1)) * Ln_first;
        Ln_first = Ln;
        Ln = Ln_next;
        ++i;
    }
    return Ln;
}

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

double Simpson(int k1, int k2, double a, double b, int n, int N) {
    double sig1 = 0, sig2 = 0, y0 = func(N, a, k1, k2, a, b), yn = func(N, b, k1, k2, a, b);

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

    return h / 3 * (2 * sig2 + 4 * sig1 + y0 + yn);
}

double Integral(int k1, int k2, double a, double b, int N) {
    double eps = 1e-8;
    int n = ceil((b - a) / sqrt(sqrt(eps)));

    double In = Simpson(k1, k2, a, b, n, N);
    double I2n = Simpson(k1, k2, a, b, 2 * n, N);
    double r = fabs(In - I2n) / 15;

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

double* ChooseMainElement(double **matrix, int N) {
    double R, *M = new double[N];
    double *res = new double[N];
    double ** transformed = new double*[N + 1];
    for (int i = 0; i < N + 1; i++) {
        transformed[i] = new double[N + 1];
    }

    int line = 0, row = 0, _line = 0;
    for (int k = 0; k < N - 1; k++) {
        int a = matrix[0][0];
        for (int i = 0; i < N; i++) {
            for (int j = 0; j < N; j++) {
                if (fabs(matrix[i][j]) > a) {
                    a = fabs(matrix[i][j]);
                    line = i;
                    row = j;
                }
            }
        }

        for (int j = 0; j < N + 1; j++) {
            transformed[_line][j] = matrix[line][j];
        }
        _line++;

        for (int i = 0; i < N; i++) {
            M[i] = -(matrix[i][row] / a);
        }

        for (int i = 0; i < N; i++) {
            for (int j = 0; j < N + 1; j++) {
                if (i == line) continue;
                matrix[i][j] += matrix[line][j] * M[i];
            }
        }

        for (int j = 0; j < N + 1; j++) {
            matrix[line][j] = 0;
        }

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


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

    for (int j = 0; j < N + 1; j++) {
        transformed[_line][j] = matrix[line][j];
    }

    for (int j = 0; j < N; j++) {
        res[j] = 0;
    }

    for (int i = N - 1; i >= 0; i--) {
        R = transformed[i][N];
        for (int j = 0; j < N; j++) {
            if ((transformed[i][j] != 0) && (res[j] != 0)) {
                R -= transformed[i][j] * res[j];
            }
        }
        for (int j = 0; j < N; j++) {
            if ((transformed[i][j] != 0) && (res[j] == 0)) {
                res[j] = R / transformed[i][j];
            }
        }
    }

    return res;
}

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

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

    return ChooseMainElement(A, N);

}

double Pm(double x, double a, double b, double *A, int n) {
    double  t = (2 * x - b - a) / (b - a);
    double y = 0;
    for (int j = 0; j < n; ++j)
        y += A[j] * L(j, t);
    return y;
}

int main() {
    double a = 1, b = 35;
    double eps = 0.01;

    int k = (b - a) * 2;
    double step = (b - a) / k;
    int n = 2;

    while (true) {
        n *= 2;
        double deviation = 0;
        double x = a;
        double *A = MakeA(a, b, n);
        for (int i = 0; i <= k; ++i) {
            double y = Function(x) - Pm(x, a, b, A, n);
            deviation += y * y;
            x += step;
        }
        deviation = sqrt(deviation / (k + 1));
        if (deviation <= eps) break;
    }

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

    double *A = MakeA(a, b, n);
    cout << "|  X  | P(x)" << endl;
    for (double i = a; i <= b; i += step) {
        printf("%3.1f %.8f\n", i, Pm(i, a, b, A, n));
    }
}