lab5.cpp

using namespace std;

#include <math.h>
#include <iomanip>
#include <iostream>

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 *SingleDivision(double **A, int N) {
	double *x = new double[N];

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

double *MakeA(double a, double b, int N) {
	double ** A = new double*[N + 1];
	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(i, j, a, b, N);
		}
		A[i][N] = Integral(i, N, a, b, N);
	}

	return SingleDivision(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 = 1e-3;

	int k = (b - a) * 2;
	double step = (b - a) / (double)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;
		}

		if (sqrt(deviation / (k + 1)) <= eps) break;
	}

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

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