123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166 |
- #define _CRT_SECURE_NO_WARNINGS
- #include <stdlib.h>
- #include <stdio.h>
- #include <math.h>
- #include "header.h"
- double integrate_left_rects(double(*f)(double), double prev_S, long int N, double a, double b, int* f_counter) {
- double S;
- long double H = (long double)(b - a) / N;
- S = 0;
- for (int i = 2;i <= N;i += 2) {
- S += f(a + (i - 1) * H, f_counter);
- }
- return H * S + (0.5) * prev_S;
- }
- void runge_left_rects(double eps, double(*f)(double), double a, double b,int* f_counter) {
- double S, S2 = 0;
- long int N = 1;
- S2 = (b - a) * f(0,f_counter);
- do {
- S = S2;
- N *= 2;
- S2 = integrate_left_rects(f, S, N, a, b, f_counter);
- } while (fabs(S2 - S) > eps);
- return S2;
- }
- double polyval(double* coeffs, double x, int n) {
- double ans = 0;
- for (int i = 0;i < n + 1;i++) {
- ans += coeffs[i] * pow(x, n - i);
- }
- return ans;
- }
- double* polyder(double* coeffs, int n) {
- double* der = malloc(sizeof(double) * (n));
- for (int i = 0;i < n;i++) {
- der[i] = coeffs[i] * (n - i);
- }
- return der;
- }
- double newton_solve(double* coeffs, int n, double a, double b) {
- double x_next = b, x_prev = b;
- double* der = polyder(coeffs, n);
- do {
- x_prev = x_next;
- x_next = x_prev - polyval(coeffs, x_prev, n) / polyval(der, x_prev, n - 1);
- } while (fabs(x_next - x_prev) >= 10e-10);
- free(der);
- return x_next;
- }
- double* solve_poly(double* coeffs, int n, double a, double b) {
- double* xx = calloc(n, sizeof(double));
- int c = 0;
- double x1 = a, x2 = a;
- while (x1 < b) {
- while (x2 < x1) {
- if (polyval(coeffs, x1, n) * polyval(coeffs, x2, n) < 0) {
- xx[c] = newton_solve(coeffs, n, x2, x1);
- c++;
- x1 = x2;
- x2 += (b - a) / 10000;
- break;
- }
- x2 += (b - a) / 10000;
- }
- x1 += (b - a) / 10000;
- }
- return xx;
- }
- double** table_xx() {
- double** ret = malloc(sizeof(double*) * 4);
- for (int n = 2;n <= 5;n++) {
- double* coeffs = malloc(sizeof(double) * (n + 1));
- double* sk = malloc(sizeof(double) * (n + 1));
- if (coeffs == NULL || sk == NULL) {
- exit(16);
- }
- coeffs[0] = 1;
- double A = 2.0 / n;
- for (int i = 1;i < (n + 1);i++) {
- sk[i] = (1.0 / (i + 1) - pow(-1, i + 1) / (i + 1)) / A;
- coeffs[i] = sk[i];
- for (int j = 1;j < i;j++) {
- coeffs[i] += coeffs[j] * sk[i - j];
- }
- coeffs[i] *= -1.0 / i;
- }
- ret[n - 2] = solve_poly(coeffs, n, -1, 1);
- free(sk);
- free(coeffs);
- }
- return ret;
- }
- double chebyshev_method(double(*f)(double), double a1, double b1, int n, double** xx_tab, int* f_counter) {
- double s = 0, a, b, ans;
- a = a1;
- b = b1;
- double* xx = xx_tab[n - 2];
- double A = 2.0 / n;
- ans = 0;
- for (int k = 0;k < n;k++) {
- ans += f((a + b) / 2 + ((b - a) / 2) * xx[k], f_counter);
- }
- s = ((b - a) / 2) * ans * A;
- return s;
- }
- double adaptive_div(double(*f)(double), double** xx_tab, double a, double b, int n, double eps, double S_n, double S_n_plus_1, int* N, int* f_counter) {
- if (fabs(S_n_plus_1 - S_n) < eps / (*N)) {
- return S_n;
- }
- else {
- double S_n_plus_1_r = chebyshev_method(f, a, (a + b) / 2, n + 1, xx_tab, f_counter);
- double S_n_r = chebyshev_method(f, a, (a + b) / 2, n, xx_tab, f_counter);
- double S_n_plus_1_l = chebyshev_method(f, (a + b) / 2, b, n + 1, xx_tab, f_counter);
- double S_n_l = chebyshev_method(f, (a + b) / 2, b, n, xx_tab, f_counter);
- *N += 1;
- return (adaptive_div(f, xx_tab, a, (a + b) / 2, n, eps, S_n_l, S_n_plus_1_l,N, f_counter) + adaptive_div(f, xx_tab, (a + b) / 2, b, n, eps, S_n_r, S_n_plus_1_r, N, f_counter));
- }
- }
- double cheb_approx(double(*f)(double), double a, double b, int n, double eps,int* f_counter) {
-
- double** xx_table = table_xx(a, b);
- int N = 1;
- return adaptive_div(f, xx_table, a, b, n, eps, 100, 1,&N, f_counter);
- }
- double* counters(double(*f)(double), double a, double b, double eps) {
- int counter_recs = 0,counter_cheb = 0;
- cheb_approx(f, a, b, 3, eps, &counter_cheb);
- runge_left_rects(eps, f, a, b, &counter_recs);
- printf("eps = %e, counter_recs = %d, counter_cheb = %d\n",eps, counter_recs, counter_cheb);
- int* counter = malloc(sizeof(int) * 2);
- counter[0] = counter_recs;
- counter[1] = counter_cheb;
-
- return counter;
- }
|