C語言實現定積分求解方法

求定積分的方法有很多種，下面是我總結的幾種比較常用的方法。

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>

#define N		3

double fun(double x)
{
	double y;
	
	y = sqrt(4-(x)*(x));
	//y = sin(x);
	
	return y;
}

/*隨機點法求定積分*/
double Darts(int n)
{
	double x, y;
	time_t t;
	int i = 0;
	int count = 0;
	
	srand((unsigned)time(&t));
	
	for (i=0; i<n; i++)
	{
		x = rand()%100/100.0;
		y = rand()%100/100.0;
		
		if (y <= 1-pow(x,2))
		{
			count++;	
		}	
	}
	
	return (double)count/(double)n;
}

/*左矩形法求定積分*/
double LeftRect(double down, double up, int n)
{
	double h, s;
	int i;
	
	/*計算步長*/
	h = (up-down)/n;
	
	s = fun(down)*h;
	
	for (i=1; i<n; i++)
	{
		s = s + fun(down+i*h)*h;	
	} 
	
	return s;
}

/*梯形公式求定積分*/
double Trape(double down, double up, int n)
{
	double h, s;
	int i = 0;
	
	/*計算步長*/
	h = (up-down)/n;
	
	s = 0.5*(fun(down)+fun(down+h))*h;
	
	for (i=1; i<n; i++)
	{
		s = s + 0.5 * (fun(down+i*h) + fun(down+(i+1)*h))*h;
	}
	
	return s;	
} 

/*復合梯形公式*/
double T(double x, double y, int z)
{
	double h, Tn;
	int i = 0;
	
	h = (y-x)/z;
	Tn = (fun(x)+fun(y))/2;
	
	for (i=0; i<z; i++)
	{
		Tn = Tn+fun(x+i*h);	
	}	
	
	Tn = Tn*h;
	
	return Tn;
} 

/*辛普生公式求定積分，公式為：S[n]=(4*T[2*n]-T[n])/3，其中T[2n],T[n]為梯形公式計算結果*/
double Simposn(double down, double up, int n)
{
	double s;
	
	/*辛普生公式*/
	s = (4*T(down, up, 2*n) - T(down, up, n))/3;
	
	return s;	
}

/*高斯公式求定積分*/
double Gass(double (*func)(double x), double a, double b, int n)
{
	int i = 0;
	
	//高斯點及其求積系數列表
	float x1[1]={0.0};                                                  
	float A1[1]={2};
	float x2[2]={-0.5573503,0.5573503};                                
	float A2[2]={1,1};
	float x3[3]={-0.7745967,0.0,0.7745967};                             
	float A3[3]={0.555556,0.888889,0.555556};
	float x4[4]={0.3399810,-0.3399810,0.8611363,-0.8611363};            
	float A4[4]={0.6521452,0.6521452,0.3478548,0.3478548};
	float x5[5]={0.0,0.5384693,-0.5384693,0.9061799,-0.9061799};       
	float A5[5]={0.5688889,0.4786287,0.4786287,0.2369269,0.2369269};   


	float *p, *t;
	switch (n)
	{
		case 1:          
			p = x1;
			t = A1;
			break;
		case 2:          
			p = x2;
			t = A2;
			break;
		case 3:          
			p = x3;
			t = A3;
			break;
		case 4:          
			p = x4;
			t = A4;
			break;
		case 5:          
			p = x5;
			t = A5;
			break;
		default :   
			printf("intput wrong!");
	}
	
	float g = 0;
	for (i=0; i<n; i++)
	{
		g += (*func)((b-a)*p[i]/2+(a+b)/2)*t[i];
	} 
	g *= (b-a)/2;
	
	return g;
}

int main(int argc, char *argv[])
{
	printf("隨機點法積分值%f\n", Darts(10000));
	
	double down, up;
	int n;
	double sum = 0;
	
	printf("積分下限:\n");
	scanf("%lf", &down);
	printf("積分上限:\n");
	scanf("%lf", &up);
	printf("分隔數目:\n");
	scanf("%d", &n);
	
	sum = LeftRect(down, up, n);
	printf("左矩形法積分值為:%f\n", sum);
	
	sum = Trape(down, up, n);
	printf("梯形公式積分值為:%f\n", sum);
	
	sum = Simposn(down, up, n);
	printf("辛普生公式積分值為:%f\n", sum);
	
	sum = Gass(fun, down, up, N);
	printf("高斯公式積分值為:%f\n", sum);
	
	return 0;
}