fork download
copy 
  1. #include <stdio.h>
  2. #include <stdlib.h>
  3. #include <math.h>
  4. #include <omp.h> // Для параллельных вычислений
  5.  
  6. // Целевая функция
  7. double target_function(double x) {
  8.     return (x * x - sin(x)) * cos(2.5 * x);
  9. }
  10.  
  11. // Вычисление полиномов Лежандра
  12. void compute_legendre_polynomials(double* polynomials, int degree, double x) {
  13.     polynomials[0] = 1; // P0(x)
  14.     if (degree > 0) {
  15.         polynomials[1] = x; // P1(x)
  16.     }
  17.     for (int n = 2; n <= degree; n++) {
  18.         polynomials[n] = ((2 * n - 1) * x * polynomials[n - 1] - (n - 1) * polynomials[n - 2]) / n;
  19.     }
  20. }
  21.  
  22. // Интегрирование методом трапеций
  23. double trapezoidal_integration(double (*func)(double, int, int), int i, int j, double a, double b, int steps) {
  24.     double h = (b - a) / steps;
  25.     double sum = 0.0;
  26.  
  27.     for (int k = 0; k <= steps; k++) {
  28.         double x = a + k * h;
  29.         double weight = (k == 0 || k == steps) ? 0.5 : 1.0;
  30.         sum += weight * func(x, i, j);
  31.     }
  32.  
  33.     return sum * h;
  34. }
  35.  
  36. // Вычисление произведения полиномов Лежандра
  37. double legendre_product(double x, int i, int j) {
  38.     double polynomials[50]; // Ограничение на степень
  39.     compute_legendre_polynomials(polynomials, fmax(i, j), x);
  40.     return polynomials[i] * polynomials[j];
  41. }
  42.  
  43. // Вычисление произведения целевой функции и полинома Лежандра
  44. double target_product(double x, int i, int j) {
  45.     double polynomials[50]; // Ограничение на степень
  46.     compute_legendre_polynomials(polynomials, i, x);
  47.     return target_function(x) * polynomials[i];
  48. }
  49.  
  50. // Построение нормальной матрицы и вектора правой части
  51. void build_normal_matrix(double* matrix, double* rhs, int degree, double a, double b, int steps) {
  52.     #pragma omp parallel for
  53.     for (int i = 0; i <= degree; i++) {
  54.         for (int j = i; j <= degree; j++) { // Используем симметрию
  55.             double value = trapezoidal_integration(legendre_product, i, j, a, b, steps);
  56.             matrix[i * (degree + 1) + j] = value;
  57.             if (i != j) {
  58.                 matrix[j * (degree + 1) + i] = value;
  59.             }
  60.         }
  61.         rhs[i] = trapezoidal_integration(target_product, i, 0, a, b, steps);
  62.     }
  63. }
  64.  
  65. // Метод Гаусса с выбором главного элемента
  66. void gauss_elimination(double* matrix, double* rhs, double* solution, int size) {
  67.     for (int i = 0; i < size; i++) {
  68.         int max_row = i;
  69.         for (int k = i + 1; k < size; k++) {
  70.             if (fabs(matrix[k * size + i]) > fabs(matrix[max_row * size + i])) {
  71.                 max_row = k;
  72.             }
  73.         }
  74.         for (int k = i; k < size; k++) {
  75.             double temp = matrix[i * size + k];
  76.             matrix[i * size + k] = matrix[max_row * size + k];
  77.             matrix[max_row * size + k] = temp;
  78.         }
  79.         double temp = rhs[i];
  80.         rhs[i] = rhs[max_row];
  81.         rhs[max_row] = temp;
  82.  
  83.         for (int k = i + 1; k < size; k++) {
  84.             double factor = matrix[k * size + i] / matrix[i * size + i];
  85.             rhs[k] -= factor * rhs[i];
  86.             for (int j = i; j < size; j++) {
  87.                 matrix[k * size + j] -= factor * matrix[i * size + j];
  88.             }
  89.         }
  90.     }
  91.     for (int i = size - 1; i >= 0; i--) {
  92.         solution[i] = rhs[i] / matrix[i * size + i];
  93.         for (int k = i - 1; k >= 0; k--) {
  94.             rhs[k] -= matrix[k * size + i] * solution[i];
  95.         }
  96.     }
  97. }
  98.  
  99. // Оценка среднеквадратичной ошибки
  100. double calculate_error(double* coeffs, int degree, double a, double b, int steps) {
  101.     double error = 0;
  102.     double step = (b - a) / steps;
  103.     for (int i = 0; i <= steps; i++) {
  104.         double x = a + i * step;
  105.         double approx = 0;
  106.         double polynomials[50];
  107.         compute_legendre_polynomials(polynomials, degree, x);
  108.         for (int j = 0; j <= degree; j++) {
  109.             approx += coeffs[j] * polynomials[j];
  110.         }
  111.         error += pow(target_function(x) - approx, 2);
  112.     }
  113.     return sqrt(error / (steps + 1));
  114. }
  115.  
  116. // Сохранение коэффициентов в файл
  117. void save_coefficients_to_file(const double* coeffs, int degree) {
  118.     FILE* file = fopen("coefficients.txt", "w");
  119.     if (!file) {
  120.         printf("Не удалось открыть файл для записи коэффициентов.\n");
  121.         return;
  122.     }
  123.     fprintf(file, "Коэффициенты многочлена Лежандра (степень %d):\n", degree);
  124.     for (int i = 0; i <= degree; i++) {
  125.         fprintf(file, "a[%d] = %lf\n", i, coeffs[i]);
  126.     }
  127.     fclose(file);
  128.     printf("Коэффициенты успешно сохранены в файл coefficients.txt\n");
  129. }
  130.  
  131. int main() {
  132.     int max_degree = 40;
  133.     double a = 0, b = 7;
  134.     int steps = 1000;
  135.     double tolerance = 0.01;
  136.  
  137.     for (int degree = 1; degree <= max_degree; degree++) {
  138.         double* matrix = (double*)malloc((degree + 1) * (degree + 1) * sizeof(double));
  139.         double* rhs = (double*)malloc((degree + 1) * sizeof(double));
  140.         double* coeffs = (double*)malloc((degree + 1) * sizeof(double));
  141.  
  142.         build_normal_matrix(matrix, rhs, degree, a, b, steps);
  143.         gauss_elimination(matrix, rhs, coeffs, degree + 1);
  144.  
  145.         double error = calculate_error(coeffs, degree, a, b, steps);
  146.         printf("Degree: %d, Error: %lf\n", degree, error);
  147.  
  148.         if (error <= tolerance) {
  149.             printf("Optimal degree: %d\n", degree);
  150.             save_coefficients_to_file(coeffs, degree);
  151.             free(matrix);
  152.             free(rhs);
  153.             free(coeffs);
  154.             break;
  155.         }
  156.         free(matrix);
  157.         free(rhs);
  158.         free(coeffs);
  159.     }
  160.  
  161.     return 0;
  162. }
  163.  
Success #stdin #stdout 2.06s 5276KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
Degree: 1, Error: 14.201725
Degree: 2, Error: 13.406704
Degree: 3, Error: 12.933561
Degree: 4, Error: 12.913922
Degree: 5, Error: 12.066789
Degree: 6, Error: 8.047765
Degree: 7, Error: 5.736353
Degree: 8, Error: 5.127578
Degree: 9, Error: 1.782684
Degree: 10, Error: 1.674198
Degree: 11, Error: 0.563190
Degree: 12, Error: 0.363394
Degree: 13, Error: 0.349523
Degree: 14, Error: 0.245960
Degree: 15, Error: 0.896356
Degree: 16, Error: 0.097016
Degree: 17, Error: 0.115056
Degree: 18, Error: 0.126031
Degree: 19, Error: 0.171374
Degree: 20, Error: 0.068058
Degree: 21, Error: 0.064201
Degree: 22, Error: 0.058142
Degree: 23, Error: 0.069285
Degree: 24, Error: 0.120590
Degree: 25, Error: 0.078761
Degree: 26, Error: 0.136907
Degree: 27, Error: 0.042303
Degree: 28, Error: 0.132902
Degree: 29, Error: 0.082553
Degree: 30, Error: 0.055851
Degree: 31, Error: 0.153969
Degree: 32, Error: 0.045659
Degree: 33, Error: 0.104272
Degree: 34, Error: 0.104347
Degree: 35, Error: 0.050368
Degree: 36, Error: 0.043915
Degree: 37, Error: 0.026945
Degree: 38, Error: 0.048973
Degree: 39, Error: 0.181343
Degree: 40, Error: 0.108479
https://ideone.com/Aegy9L
language:
C (gcc 8.3)
created:
22 days ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!