自适应积分的函数的编写

不念不忘少年蓝@ 2022-10-30 01:21 152阅读 0赞

最近在写软件，需要用到积分。我们采用多个矩形累加的方法，很容易实现积分的功能，但是这个方法有一个问题：我们不知道把a~b之间分成多少份，也就是不知道n取多少。如果n太小，计算结果精度较低，n太大，则计算量太大，且过高的精度在工程上也没有意义。

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dhbmdlaWwwMDc_size_16_color_FFFFFF_t_70]

这里我们采用自适应的方式，让程序自己选择n的大小。我们首先定义两个要积分的函数f1和f2。

double f1(double x)
    {
        
    	return 1/x;
    }
    
    double f2(double x)
    {
        
    	return 1/(x * x);
    }

下面是积分函数：  
4个参数的意义：  
a：积分的起点；  
b：积分的终点；  
n：从a到b，一共分成n份。这个n是另一个函数指定的。  
myFun: 函数指针，可以用来指定对f1还是f2进行积分。

//从a到b的积分，但不一定精确 
    //第三个参数是一个函数指针，申明的时候一定要注意函数名后面带上参数列表
    //这时，函数名就可以被当做参数传递到函数中
    double jifen(double a, double b, int n, double (*myFun)(double))
    {
        	
    	//把被积分的区域a~b等分成n块。dx就是每一块的宽度。
    	//在数学上，并没有要求每一块的宽度相等。这里为了简单起见，采取等分 
    	double dx = (b - a) /n;  
    	
    	double sum = 0;
    	double x; //具体的x的位置 
    	
    	for(int i = 0; i < n ; i++)
    	{
        
    		x = a + dx * i;	//动态更新x的值 
    		
    		//myFun相当于C#中的delegation (委托)
    		//直接使用形参中的函数指针调用不同的函数，例如f1和f2 
    		//矩形法计算一小块矩形的面积，并进行累加 
    		sum += myFun(x) * dx; 
    	}
    
    	return sum;
    }

在编写完基本的积分函数之后，我们开始编写能决定jifen()中n的值的函数jifen2()。这是这段程序的关键所在。  
在积分的定义上，并没有规定把a和b之间的距离进行等分，只规定每一份之间的距离都是无穷小 (dx->0)。为了简便起见，我们把a和b之间的距离进行等分。  
我们采用的思路是：如果把a和b之间的距离等分成n份进行计算积分值，得到一个结果r1；如果把a和b之间的距离等分成10 \* n份进行计算积分值，得到一个结果r2。随着n的增大，积分值距离真正的精确值越来越近，提高的幅度也越来越小。如果r1和r2之间的差距太大，则需要继续增大n的值。如果r1和r2之间的差距 < 0.1%，则我们认为r1和r2都是比较准确的值。r1和r2之间的差距很小时，r1和r2都与真实值很接近，它们与真实值之间的差距也很小。据此，我们写下如下程序：我们要求两个结果之间的差距 < 0.1%。

//最大误差为0.1%的积分 
    //第三个参数是一个函数指针，申明的时候一定要注意函数名后面带上参数列表 
    //这时，函数名就可以被当做参数传递到函数中 
    double jifen2(double a, double b, double (*myFun)(double))
    {
        
    	int n = 1000; //从1000开始计算。10、100等数据时，误差很大，一般需要到5000或者10000才能得到比较精确的结果 
    	double r1, r2; //两次积分得到的结果 
    	double error = 100; //初始的error 
    	
    	while(!(-0.001 < error && error < 0.001)) //误差太大就继续提高n的值继续计算 
    	{
        
    		r1 = jifen(a, b, n, myFun); //等分成n份的时候的结果 
    		r2 = jifen(a, b, n * 10, myFun);//等分成10n份的时候的结果
    		error = (r2 - r1) / r2; //两次计算之间的误差 
    		
    		printf("分成n = %d份, 结果r1 = %f, r2 = %f, 两个结果的误差：%f,\r\n\r\n", n, r1, r2, error);
    		n *= 5; //n扩大的幅度是n = n * 5 
    	}
    	
    	printf("积分的结果：sum = %f\r\n", r1);
    	return r1;
    }

完成了函数的定义之后就可以进行简单的测试了：

int main()
    {
        
    	//函数指针是指向函数的指针变量，即本质是一个指针变量
    	//double (*myFun) (double a); /*声明一个函数指针 */
    	jifen2(1.0, 10.0, f1);
    	
    	printf("Hello\r\n");
    	return 0;
    }

运行结果如下：

![在这里插入图片描述][2021021616541918.png]  
我们用牛顿-莱布尼茨公式计算一下y = 1/x在1~10之间的积分值：

![在这里插入图片描述][20210216165853930.png]

对比这两个数据，我们的程序的计算值与真实值之间的误差为0.035%，符合我们的要求。这时n的值是5000。

我们也可以把f2作为参数代入jifen2()进行计算：

jifen2(1.0, 10.0, f2);

运行结果为：  
![在这里插入图片描述][20210216170338975.png]

我们同样用数学方法计算一下：

![在这里插入图片描述][20210216170813547.png]  
计算值与真实值之间的误差为0.099%，也符合我们的要求。这时n的值也是5000。上面两个实例说明，我们的计算方法是有效且准确的：用这个方法可以得到与真实值误差约为0.1%的结果 (实测对于不同的待积分函数，误差有可能超过0.1%，但不超过0.2%)。当然，上述程序也并非没有问题，我们在计算两个n时，差别为10倍，最终n=5000时，其实对n=50000也进行了计算，所以可以降低两者的倍数，例如降低为2倍；而且前一次的计算结果也可以用于下一轮计算，这样也可以提高计算效率。

总体代码如下：

/*
    本文首先定义了f1和f2这两个函数。
    接下来在主函数中调用jifen2，同时把f1这个函数名作为参数传递给jifen2。jifen2通过精度上的判断，
    决定给jifen函数一个合适的n，然后jifen函数对f1函数，在a~b进行积分，使得最终误差<0.1% 
    
    这里使用了函数指针，
    函数指针可以作为参数传递，这是要主要在函数名后面带上参数列表。
    也可以直接用函数指针加上参数，调用不同的函数。这就相当于C#中的委托 (delegation)
    使用函数指针之后，就可以把jifen2和jifen这两个函数完全封装起来了，不需要为不同的被积分函数写不同的jifen函数 
    */
    
    #include <stdio.h>
    #include <stdlib.h>
    #include <math.h>
    
    
    double f1(double x)
    {
        
    	return 1/x;
    }
    
    double f2(double x)
    {
        
    	return 1/(x * x);
    }
    
    //从a到b的积分，但不一定精确 
    //第三个参数是一个函数指针，申明的时候一定要注意函数名后面带上参数列表
    //这时，函数名就可以被当做参数传递到函数中
    double jifen(double a, double b, int n, double (*myFun)(double))
    {
        	
    	//把被积分的区域a~b等分成n块。dx就是每一块的宽度。
    	//在数学上，并没有要求每一块的宽度相等。这里为了简单起见，采取等分 
    	double dx = (b - a) /n;  
    	
    	double sum = 0;
    	double x; //具体的x的位置 
    	
    	for(int i = 0; i < n ; i++)
    	{
        
    		x = a + dx * i;	//动态更新x的值 
    		
    		//myFun相当于C#中的delegation (委托)
    		//直接使用形参中的函数指针调用不同的函数，例如f1和f2 
    		//矩形法计算一小块矩形的面积，并进行累加 
    		sum += myFun(x) * dx; 
    	}
    
    	return sum;
    }
    
    //最大误差为0.1%的积分 
    //第三个参数是一个函数指针，申明的时候一定要注意函数名后面带上参数列表 
    //这时，函数名就可以被当做参数传递到函数中 
    double jifen2(double a, double b, double (*myFun)(double))
    {
        
    	int n = 1000; //从1000开始计算。10、100等数据时，误差很大，一般需要到5000或者10000才能得到比较精确的结果 
    	double r1, r2; //两次积分得到的结果 
    	double error = 100; //初始的error 
    	
    	while(!(-0.001 < error && error < 0.001)) //误差太大就继续提高n的值继续计算 
    	{
        
    		r1 = jifen(a, b, n, myFun); //等分成n份的时候的结果 
    		r2 = jifen(a, b, n * 10, myFun);//等分成10n份的时候的结果
    		error = (r2 - r1) / r2; //两次计算之间的误差 
    		
    		printf("分成n = %d份, 结果r1 = %f, r2 = %f, 两个结果的误差：%f,\r\n\r\n", n, r1, r2, error);
    		n *= 5; //n扩大的幅度是n = n * 5 
    	}
    	
    	printf("积分的结果：sum = %f\r\n", r1);
    	return r1;
    }
    
    int main()
    {
        
    	//函数指针是指向函数的指针变量，即本质是一个指针变量
    	//double (*myFun) (double a); /*声明一个函数指针 */
    	jifen2(1.0, 10.0, f2);
    	
    	printf("Hello\r\n");
    	return 0;
    }

以下代码是按照上述缺点对jifen2的改进：

//最大误差为0.1%的积分 
    //第三个参数是一个函数指针，申明的时候一定要注意函数名后面带上参数列表 
    //这时，函数名就可以被当做参数传递到函数中 
    //改进点：通过r1=r2的方式，避免重复计算r1 
    //简化while循环的判断条件 
    double jifen2(double a, double b, double (*myFun)(double))
    {
        
    	int n = 1000; //从1000开始计算。一般需要到5000或者10000才能得到比较精确的结果 
    	double r1, r2; //两次积分得到的结果 
    	int times = 3; //数据增大的倍数 
    	r1 = jifen(a, b, n, myFun);//等分成n份的时候的结果
    	double error = 1; //初始的error 
    	
    	//误差太大就提高n的值继续计算
    	while(error < -0.001 || error > 0.001) 
    	{
        
    		n *= times; //n扩大的幅度是n = n * times 
    		
    		//等分成n * times份的时候的结果。r2的结果一定比r1更精确 
    		r2 = jifen(a, b, n, myFun); 
    		error = (r2 - r1) / r2; //两个result之间的error 
    		
    		printf("分成n = %d份, 结果r1 = %f, r2 = %f, 两个结果的误差：%f,\r\n\r\n", n, r1, r2, error);
    		
    		//r2的精度一定比r1更高，所以用r1 = r2的方式更新r1 
    		//这样可以避免重复计算r1。但是不要把r1和r2存放到数组中，因为不知道数组的容量是多少 
    		r1 = r2;  
    	}
    	
    	printf("积分的结果：sum = %f\r\n", r2);
    	return r2;
    }

[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dhbmdlaWwwMDc_size_16_color_FFFFFF_t_70]: /images/20221024/1c1d21f1b0af43d385e8a47e0227b638.png
[2021021616541918.png]: /images/20221024/0b6c44a9154a41d3aec01e2f330bf8ab.png
[20210216165853930.png]: /images/20221024/8bf556851d594078abefd3a3252378bd.png
[20210216170338975.png]: /images/20221024/2495066d1c03405db3f9f586fb09dace.png
[20210216170813547.png]: /images/20221024/aa5e157cfb134779b0b9d0bf73774b07.png