Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- using namespace std;
- #include <math.h>
- #include <iomanip>
- #include <iostream>
- double Function(double x) {
- return exp(pow(x - 15, 0.25)) * sin(1.2 * x);
- }
- double Legendre(int n, double x) {
- double Ln_next, Ln = x, Ln_first = 1;
- if (n == 0) return Ln_first;
- for (int i = 1; i < n; i++) {
- Ln_next = (1.0 * (2 * i + 1) / (i + 1)) * x * Ln - (1.0 * i / (i + 1)) * Ln_first;
- Ln_first = Ln;
- Ln = Ln_next;
- }
- return Ln;
- }
- double func(int N, double t, int k1, int k2, double a, double b) {
- double x = (2 * t - a - b) / (b - a);
- if (k2 == N) return Function(t) * Legendre(k1, x);
- else return Legendre(k1, x) * Legendre(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 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;
- }
- double *CompletePivoting(double **matrix, int n, int m) {
- int p, q, line_new = 0;
- double M[n], t[n][m];
- double *x = new double[n];
- for (int k = 0; k < n - 1; k++){
- double a = 0;
- for (int i = 0; i < n; i++) {
- for (int j = 0; j < m - 1; j++) {
- if (fabs(matrix[i][j]) > a) {
- a = fabs(matrix[i][j]);
- p = i;
- q = j;
- }
- }
- }
- for (int j = 0; j < m; j++) {
- t[line_new][j] = matrix[p][j];
- }
- line_new++;
- for (int i = 0; i < n; i++) {
- M[i] = -matrix[i][q] / a;
- }
- for (int i = 0; i < n; i++) {
- for (int j = 0; j < m; j++) {
- if (i != p) {
- matrix[i][j] += matrix[p][j] * M[i];
- }
- }
- }
- for (int j = 0; j < m; j++) {
- matrix[p][j] = 0;
- }
- for (int i = 0; i < n; i++) {
- matrix[i][q] = 0;
- }
- }
- for (int i = 0; i < n; i++) {
- if (matrix[i][m - 1] != 0) {
- p = i;
- }
- }
- for (int j = 0; j < m; j++) {
- t[line_new][j] = matrix[p][j];
- }
- for (int j = 0; j < m-1; j++) {
- x[j] = 0;
- }
- for (int i = n - 1; i >= 0; i--){
- double R = t[i][m - 1];
- for (int j = 0; j < m - 1; j++) {
- if ((t[i][j] != 0) && (x[j] != 0)) {
- R -= t[i][j] * x[j];
- }
- }
- for (int j = 0; j < m - 1; j++) {
- if ((t[i][j] != 0) && (x[j] == 0)) {
- x[j] = R / t[i][j];
- }
- }
- }
- return x;
- }
- double *MakeA(double a, double b, int N) {
- double ** A = new double*[N];
- for (int i = 0; i < N; i++)
- A[i] = new double[N + 1];
- for (int i = 0; i < N; i++) {
- for (int j = 0; j < N; j++) {
- A[i][j] = A[j][i] = Runge(i, j, a, b, N);
- }
- A[i][N] = Runge(i, N, a, b, N);
- }
- return CompletePivoting(A, N, N + 1);
- }
- double Pm(double x, double a, double b, double *A, int n) {
- double t = (2 * x - b - a) / (b - a), s = 0;
- for (int j = 0; j < n; j++) {
- s += A[j] * Legendre(j, t);
- }
- return s;
- }
- int main() {
- double a = 20, b = 30;
- 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));
- }
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement