123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120 |
- #include <stdio.h>
- #include <stdlib.h>
- #include <math.h>
- #define A -3
- #define B 3
- double Func(double x, int *count){
- (*count)+=1;
- return x*x*x*x*x - 5.2*x*x*x + 5.5*x*x - 7*x;
- }
- double trapezoidalIntegral(int n, double previousResult, double (*func)(double x, int *count), int *count) {
- double result = 0.0;
- double h;
- int i;
- h = (B - A) / (double)n;
- if (n == 1){
- result = (func(A, count) + func(B, count)) * h / 2;
- for (i = 1; i < n; i++) {
- result += func(A + i * h, count) * h;
- }
- }else{
- result = previousResult/2.0;
- for (i = 1; i < n; i += 2) {
- result += func(A + i * h, count) *h;
- }
- }
- return result;
- }
- double Integral1(double eps, int *n, double (*func)(double x, int *count), int *count){
- (*n) = 1;
- double next= trapezoidalIntegral(*n, 0, func, count);
- double prev;
-
- do{
- (*n) = (*n)*2;
- prev = next;
- next = trapezoidalIntegral(*n, prev, func, count);
- } while (fabs(next - prev)/3 >= eps);
-
- return next;
- }
- double rado_2(double a, double b, double (*func)(double x, int *count), double fb, int *count){
- double x1 = (2*a+b)/3;
- return 3*(b-a)*func(x1, count)/4+(b-a)*fb/4;
- }
- double rado_3(double a, double b, double (*func)(double x, int *count), double fb, int *count){
- double x1 = ((6+sqrt(6))*a + (4-sqrt(6))*b)/10;
- double x2 = ((6-sqrt(6))*a + (4+sqrt(6))*b)/10;
- double A1 = (16-sqrt(6));
- double A2 = (16+sqrt(6));
- double A3 = 4;
- return (b-a)*(A1*func(x1, count) + A2*func(x2, count) + A3*fb)/36;
- }
- double Integral2(double eps, double a, double b, int *n, double (*func)(double x, int *count), int *count){
- int i;
- double left, right, step, fright;
- double R2, R3;
-
- (*n) = 1;
- fright = func(b, count);
-
- R2 = rado_2(a,b,func,fright, count);
- R3 = rado_3(a,b,func,fright, count);
-
- while (fabs(R2 - R3) >= eps){
- R2 = 0.0;
- R3 = 0.0;
- (*n) = (*n)*2;
- step = (b-a)/(*n);
- right = a;
-
- for (i=0;i<(*n);i++){
- left = right;
- right = left + step;
- fright = func(right, count);
-
- R2 += rado_2(left, right, func, fright, count);
- R3 += rado_3(left, right, func, fright, count);
- }
- }
- return R3;
- }
- int main(){
- FILE *f = fopen("data.txt", "w");
- if (f == NULL) {
- printf("Error opening file!\n");
- return 1;
- }
-
- int i, k, n, h;
- double Iexact = 99.0;
- double eps = 0.1;
- double result;
- for(i = 0; i<9; i++){
- k = 0;
- result = Integral1(eps, &n, Func, &k);
- fprintf(f, "%.15f %d\n", fabs(Iexact-result), k);
- eps = eps/10;
- }
-
- fprintf(f, "--------------------\n");
- eps = 0.1;
- for(i = 0; i<9; i++){
- h = 0;
- result = Integral2(eps, A, B, &n, Func, &h);
- fprintf(f, "%.15f %d\n", fabs(Iexact - result), h);
- eps = eps*0.1;
- }
-
- fclose(f);
- return 0;
- }
|