2020-11-30

C语言–求定积分函数（龙贝格积分算法）

求定积分函数，经过计算，结果的精度还算可以，误差 不大

/********************************************  
    龙贝格积分算法
 ************************************************************************/    
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define N 20
//#define MAX 10     //数组存的最大行数
//#define a 0.0000001    //积分下限
//#define b 4.0     //积分上限
#define eps 1e-7 //精度
 
double f(double x)//所求积分公式
{
    return sin(x);
}
 
double computeT(double aa, double bb, long int n)//复化梯形公式，aa、bb区间两端的值，n为分成的份数
{
    int i;
    double sum=0, h = (bb - aa)/n;//h相当于宽度，把区间分成了n份
    for (i = 1; i < n; i++)
        sum += f(aa + i * h);//sum求的是高度
	
	    sum += (f(aa) + f(bb)) / 2;
    return (h * sum);//高度乘以宽度，得出面积
}
 
double lbg(double a, double b)
{
	int MAX=20;
	int i;
    long int n = N;//N=20,为分成的段数
	long int m=0;
	double M[1][1];
    double T[20 + 1][2];//定义一个T[][]值  double T[MAX + 1][2];
    T[0][1]=computeT(a,b,n);//梯形算法
      n*= 2;//为分成的段数
	  //大型for循环，MAX（数组存的最大行数）为20
    for (m = 1; m < MAX; m++)
    {
        for (i = 0; i < m; i++)
        {
            T[i][0] = T[i][1];
        }

	    T[0][1]=computeT(a,b,n);//梯形算法
        n *= 2;//为分成的段数
 
        for (i = 1; i <= m; i++) //T的m(h)
            T[i][1] = T[i - 1][1] + (T[i - 1][1] - T[i - 1][0]) / (pow(2.0, 2 * m) - 1);

        if ((T[m - 1][1] < T[m][1] + eps) && (T[m - 1][1] > T[m][1] - eps))
        {
	        M[0][0]=T[m][1];
            //printf("计算的数为:%lf\n", T[m][1]); //输出
        }
    }
	return M[0][0]; 

   
}
int main()
{
	double zh;
	zh=lbg(0.0001,6.2831852);
	printf("%f",zh);
	system("pause");

}

以下是收集的相关知识，帮助理解！
龙贝格算法
https://www.cnblogs.com/duye/p/9118957.html
复合梯形求积公式
https://wenku.baidu.com/view/425cd3ea76eeaeaad1f330df.html
龙贝格算法（matlab）
https://wenku.baidu.com/view/3929411750e2524de4187e01.html
龙贝格求定积分c语言
https://blog.csdn.net/landcruiser007/article/details/79314989
数值作业:龙贝格算法计算积分C语言实现
https://blog.csdn.net/Chen_dSir/article/details/71270316?utm_medium=distribute.pc_relevant.none-task-blog-BlogCommendFromBaidu-4.control&depth_1-utm_source=distribute.pc_relevant.none-task-blog-BlogCommendFromBaidu-4.control