integral_simpson.c

平时完成作业做的几个小程序。有数值计算方面的梯形公式求积、simpson算法求积、Jacobi算法求解线性方程组、Gauss-Seidel法求解线性方程组

/*自适应辛普森公式求积分的程序
  author:黄翔 date：2008-5-12
*/

#include <stdio.h>
#include <math.h>
// 辛普森公式的逐次分半法求积函数
double integral(
                double a,//积分区间起始值
                double b,//积分区间终止值
                double(*f)(double),//被积函数
                double eps//积分允许误差
                )

{
        double n   = 1,k =1 ;
        double h   = (b-a);
		double s0  = f(a) + f(b) , s1 = 0 , s2 = 0;
        double In  = 0;
        double I2n = 0;
        while(true)
        {  
                long i,m=n/2;		
                for( i=1;i<m;++i){ 			
						s2 += f(a+(2*i-1)*h);
                }		
                for(i=1;i<m;++i){ 			
						s1 += f(a+(2*i)*h);
                }						
                I2n = (s0 + 2*s1 + 4*s2)*(h/6);
                if(fabs(I2n-In)<15*eps) break;
                In = I2n;
                n *= 2;
				k += 1;
                h /= 2;
        }
        return I2n;
}

// 被积函数f1(x)=log(x)
double func1(double x)
{
        return log(x);
}

// 被积函数f2(x)=1/(1+x*x)
double func2(double x)
{
        return(1/(1+x*x));
}

// 被积函数f3(x)=1/x
double func3(double x)
{
        return(1/x);
}


int main()
{        
		double eps;
		printf("请给出求解的精度，用小数表示：\n");
		scanf("%lf",&eps);
        printf("%s%.8f\n","f1(x)=log(x)在[1,10]的积分值为:",integral(1,10,func1,eps));
		printf("%s%.8f\n","f2(x)=1/(1+x*x)在[-4,4]的积分值为:",integral(-4,4,func2,eps));
		printf("%s%.8f\n","f3(x)=1/x在[1,2]的积分值为:",integral(1,2,func3,eps));
        return 0;
}