integral_trapezium.c

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

/*自适应梯形公式求积分的程序
  author:黄翔 date：2008-5-11
*/

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

{
        double n   = 1;
        double h   = (b-a)/2;
        double Tn  = h*(f(a)+f(b));
        double T2n = 0;
        while(true)
        {
                double sumf = 0;
                for(int i=1;i<=n;++i)
                {
                        sumf += f(a+(2*i-1)*h);
                }
                T2n = Tn/2+h*sumf;
                if(fabs(T2n-Tn)<3*eps)break;
                Tn = T2n;
                n *= 2;
                h /= 2;
        }
        return T2n;
}

// 被积函数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;
}