NumericalMethods/Курсовая работа/kp.c

#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;    
}