123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163 |
- #define _CRT_SECURE_NO_WARNINGS
- #include <stdlib.h>
- #include <stdio.h>
- #include <math.h>
- #include "header.h"
- 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) {
- 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]);
- }
- s = ((b - a) / 2) * ans * A;
- return s;
- }
- double adaptive_div(double(*f)(double),double** xx_tab,double a,double b,int n,int power, double S_n,double S_n_plus_1,int *N,double* len) {
- if (fabs(S_n_plus_1 - S_n) < pow(10,-power)/(*N)) {
- return S_n;
- }
- else {
- double S_n_plus_1_r = chebyshev_method(f, a, (a + b) / 2, n+1, xx_tab);
- double S_n_r = chebyshev_method(f, a, (a + b) / 2, n, xx_tab);
- double S_n_plus_1_l = chebyshev_method(f, (a + b) / 2, b, n+1, xx_tab);
- double S_n_l = chebyshev_method(f, (a + b) / 2, b,n, xx_tab);
- *N += 1;
-
- if (*len > (b - a) / (*N)) {
- *len = (b-a)/(*N);
- }
- return (adaptive_div(f, xx_tab, a, (a + b) / 2, n, power, S_n_l, S_n_plus_1_l,N,len) + adaptive_div(f, xx_tab, (a + b) / 2, b, n, power, S_n_r, S_n_plus_1_r,N,len));
- }
-
- }
- double adaptive_div_Richardson(double** xx_tab, double a, double b, double(*f)(double), int n, int power, double S_prev, double S_new, int* N, double* len) {
- if (fabs(S_new - S_prev)/(pow(2,n)-1) < pow(10, -power) / ( (*N))) {
- return S_new;
- }
- else {
- double S_newl = chebyshev_method(f, a, (a + b) / 2, n + 1, xx_tab);
- double S_prevl = chebyshev_method(f, a, (a + b) / 2, n, xx_tab);
- double S_newr = chebyshev_method(f, (a + b) / 2, b, n + 1, xx_tab);
- double S_prevr = chebyshev_method(f, (a + b) / 2, b, n, xx_tab);
- *N += 1;
- if (*len > (b - a) / (*N)) {
- *len = (b - a) / (*N);
- }
- return (adaptive_div_Richardson(xx_tab, a, (a + b) / 2, f, n, power, S_prevl, S_newl, N, len) + adaptive_div_Richardson(xx_tab, (a + b) / 2, b, f, n, power, S_prevr, S_newr, N, len));
- }
- }
- double find_approx_integral(double(*f)(double), double a, double b, int n, int power, int* counter,double* length,int richardson) {
- int N = 1;
- double** xx_table = table_xx(a, b);
-
- double S2;
- double len = 10000;
- if (richardson) {
- S2 = adaptive_div_Richardson(xx_table, a, b, f, n, power, 1000000, 100, &N, &len);
- }
- else {
- S2 = adaptive_div(f,xx_table, a, b, n, power, 1000000, 100, &N, &len);
- }
- printf("N = %ld, n = %d, S = %.15lf, eps = %e, len=%.16lf\n", N,n, S2, pow(10, -power),len);
- *counter = N;
- *length = len;
- return S2;
- }
|