〖数学算法〗积分算法（一）

当我上小学的时，就学习了球的体积公式V=(4/3)πR³，当时觉得它实在太神奇了，是不是求得这个公式的人把一个铁球熔成铁水，放在一个矩形容器中求的？直到大上学才知道是利用积分算得的，当然微积分这个东西对于包括我在内的广大同学们来说可能是恨大于爱，但不可否认是积分在几乎所有理工学科都有着无可替代的作用，所以博主就写一写求积分的算法，由于算法过多，为了避免篇幅过长，给读者造成疲劳感，我决定分4篇写积分算法。

唯一方便统一，本篇各算法均以这个最基本式子的来作例子

1.随机投点法（蒙特卡洛算法）

在求圆周率的文章中已经提及过一次此方法

链接：http://blog.csdn.net/nash_/article/details/8199357

在a到b和函数组成的矩形的范围内，随机投N个点，落到绿色阴影点的个数为M个，对于此图来说便可以容易得知积分的值（绿色阴影）为(M/N)*矩形面积。

代码清单：

[java] view plaincopy

1. public class JiFen {  
2.       
3.     public static void main(String[] args){  
4.       
5.         int N = 1000000;  
6.         int count = 0;  
7.         for(int i = 0; i < N; i++){  
8.             double x = Math.random();  
9.             double y = Math.random();  
10.             if(f(x) >= y)  
11.                 count++;  
12.         }  
13.         System.out.println((double)count/N);  
14.     }  
15.   
16.     private static double f(double x) {  
17.           
18.         return x*x;  
19.     }  
20.       
21. }  

输出：0.333302

2.另一种蒙特卡洛法

第一种方法视乎情况比较特殊，如果是积分形式如下图：

对于普通情况，如果用第一种算法，就要判断随机点是投在了x轴的上方还是下方，对于矩形的选取还要分f(a),f(b)是否同号，以及f(a),f(b)的绝对值的大小，比较麻烦。于是产生了另一种投点法：在a到b的范围内随机生成N个点x1~xn，则积分的值为(f(x1) + f(x2) + ...+ f(xn)) * (b-a)/m。

代码清单：

[java] view plaincopy

1. public class JiFen_toudian2 {  
2.       
3.     public static void main(String[] args){  
4.       
5.         double a = 0;  
6.         double b = 1;  
7.           
8.         double area = getArea(a,b);  
9.           
10.         System.out.println(area);  
11.     }  
12.     public static double getArea(double a, double b){  
13.         int sign = 1;//正负号标志位  
14.         if(a > b){  
15.             double t = a;  
16.             a = b;  
17.             b = t;  
18.             sign = -1;  
19.         }  
20.         int N = 10000;  
21.         double sum = 0;  
22.         for(int i = 0; i < N; i++){  
23.             double x = a + (b-a)*Math.random();//随机生成1个在（a,b)范围内的点  
24.             sum += f(x);  
25.         }  
26.         return (sum * (b-a)/N) * sign;   
27.     }  
28.     private static double f(double x) {  
29.           
30.         return x*x;  
31.     }  
32.       
33. }  

输出：0.33966325465503505

3.定义求积分法

回忆一下最初是怎么求积分的，是把一块面积分解成N个小矩形，然后求面积和，最后求极限。我们就利用这种方式来实现它，但我们的N毕竟有限，为了结果更准确，把求矩形面积改成求梯形面积（当然矩形也是可以的），如下图：

把(a,b)分成N等分，积分值等于S1+S2+....+Sn，其中Si = (f(xi) + f(x(i+1))) * (b-a)/n / 2(矩形面积公式）

有了之前的基础，就可以比较容易的写程序了。

代码清单：

[java] view plaincopy

1. public class JiFen3 {  
2.   
3.     /** 
4.      * @param args 
5.      */  
6.     public static void main(String[] args) {  
7.           
8.         double a;  
9.         double b;  
10.         double y;  
11.         a = 0;  
12.         b = 1;  
13.         y = getArea(a, b);  
14.         System.out.println(y);  
15.     }  
16.     static double f(double x){  
17.         return x*x;  
18.     }  
19.     static double getArea(double a, double b){  
20.         int sign = 1;  
21.         if(a > b){  
22.             double t = a;  
23.             a = b;  
24.             b = t;  
25.             sign = -1;  
26.         }  
27.         double h;  
28.         double x;  
29.         double sum = 0;  
30.         int n = 10000;  
31.         h = Math.abs(a - b)/n;  
32.         x = a;  
33.         for(int i = 0; i < n; i++){  
34.               
35.             sum += f(x) + f(x + h);  
36.             x = x + h;  
37.         }  
38.         return sum * h / 2 * sign;  
39.     }  
40. }  

输出：0.3333333349999429

4.变步长梯形求积分法

用定义求积分法，要确定积分的误差是很难的，我们无法找到一个正好的N值符合我们要求的误差值，所以就需对定义法进行改进，改进的方式就是把原来固定的步长改为变化的步长，利用二分法，如下图：

                                                      图4-1

                                                    图4-2

                                               图4-3

。

。

。

我们要分到什么时候呢？ 分到 后一个面积和 - 前一个面积和 < 规定误差 时。这样我们就达到了精确的目的。

代码清单：

[java] view plaincopy

1. public class Jifen_bianchang {  
2.   
3.     static double e = 0.00001;// 误差  
4.   
5.     public static void main(String[] args) {  
6.         double a = 0;// 积分下限  
7.         double b = 1;// 积分上限  
8.         double area = getArea(a, b);  
9.         System.out.println(area);  
10.     }  
11.   
12.     public static double getArea(double a, double b) {  
13.   
14.         int sign = 1;// 正负标志位  
15.         if (a > b) {  
16.             double t = a;  
17.             a = b;  
18.             b = t;  
19.             sign = -1;  
20.         }  
21.         double s1 = 0;// 前一个面积和  
22.         double s2 = 0;// 后一个面积和  
23.         double h = b - a;  
24.         s2 = getOneArea(a, b, h);  
25.   
26.         for (int i = 2; Math.abs(s1 - s2) > e; i *= 2) {  
27.             double t = h / i;// 每个梯形高  
28.             double sum = 0;  
29.             double x = a;  
30.             for (int j = 0; j < i; j++) {// 求梯形和  
31.                 sum += getOneArea(x, x + t, t);  
32.                 x = x + t;  
33.             }  
34.             s1 = s2;  
35.             s2 = sum;  
36.         }  
37.         return sign * s2;  
38.     }  
39.   
40.     public static double getOneArea(double a, double b, double h) {  
41.   
42.         return (f(a) + f(b)) * h / 2;  
43.     }  
44.   
45.     public static double f(double x) {  
46.         return x * x;  
47.     }  
48.   
49. }  

输出：0.33333587646484375

==================================================================================================

  作者：nash_  欢迎转载，与人分享是进步的源泉！

  转载请保留原文地址：http://blog.csdn.net/nash_/article/details/8200349

===================================================================================================