蒙特卡罗法求积分

文章目录

• 背景介绍
• 要求
• 题目
• 代码
• • 顺序代码
  • 并行代码
• 误差

背景介绍

蒙特卡罗法是一类随机算法的统称。随着二十世纪电子计算机的出现，蒙特卡洛法已经在诸多领域展现出了超强的能力。在机器学习和自然语言处理技术中，常常被用到的MCMC也是由此发展而来。蒙特卡洛法最为常见的一种应用就是使用概率方法近似求解定积分。

在本实验中，已知二维平面上一条连续曲线f(x)，需要求解一个区间内曲线的面积。需要使用蒙特卡洛法，随机采集区间内的N个点，按照公式面积 = ∫ x 2 x 1 f ( x ) d x = l i m N → ∞ 1 N ∑ r = 1 N f ( x r ) ( x 2 − x 1 ) =\int_{x_2}^{x_1}f(x)dx=lim_{N \to \infty}\frac{1}{N}\sum_{r=1}^{N}f(x_r)(x_2-x_1) =∫x2​x1​​f(x)dx=limN→∞​N1​∑r=1N​f(xr​)(x2​−x1​)，使用 1 N ∑ r = 1 N f ( x r ) ( x 2 − x 1 ) \frac{1}{N}\sum_{r=1}^{N}f(x_r)(x_2-x_1) N1​∑r=1N​f(xr​)(x2​−x1​)进行近似求解，并通过MPI并行化计算过程。

要求

• 自行设定二维平面上的曲线方程f(x)以及界定区域[x1,x2]；
• 计算所选定曲线在界定区域内的精确的面积大小；
• 给出蒙特卡罗法求解该曲线面积的顺序代码和并行代码；
• 统计不同的精度（N不同）下所求曲线面积和相应的误差；
• 在不同精度下，对并行代码进行性能分析（执行时间、加速比、计算/通信比）

题目

1. 选择曲线方程 f ( x ) = 9 x − x 2 f(x)=9x-x^2 f(x)=9x−x2，以及界定区域[0,9]。
2. 使用积分计算方程在[0,9]上的面积为 ∫ 0 9 ( 9 x − x 2 ) d x = 121.5 \int_0^9(9x-x^2)dx=121.5 ∫09​(9x−x2)dx=121.5
3. 使用蒙特卡罗法求解，近似公式为面积 = 1 N ∑ r = 1 N f ( x r ) ( x 2 − x 1 ) =\frac{1}{N}\sum_{r=1}^{N}f(x_r)(x_2-x_1) =N1​∑r=1N​f(xr​)(x2​−x1​)

代码

顺序代码

1.	#include <stdio.h>  
2.	#include <stdlib.h>  
3.	#include <time.h>  
4.	  
5.	#define N 100000000  
6.	  
7.	double rand_v(int x1, int x2);  
8.	int main(void)  
9.	{  
10.	    int i, x2=9, x1=0;  
11.	    double area, xr, acc, sum = 0;  
12.	    clock_t clk;  
13.	    srand(time(0));  
14.	  
15.	    clk = clock();  
16.	    for (i=0; i < N; i++) {  
17.	        xr = rand_v(x1, x2);  
18.	        sum += (9 - xr)*xr;  
19.	    }  
20.	    area = (sum / N) * (x2 - x1);  
21.	    clk = clock() - clk;  
22.	      
23.	    acc = (area - 121.5)/121.5;  
24.	    acc = acc>0 ? acc : -acc;  
25.	    printf("area = %f\tacc = %f\ttime = %f\n",   
26.	            area, acc, (double)clk/CLOCKS_PER_SEC);  
27.	    return 0;  
28.	}  
29.	double rand_v(int x1, int x2)  
30.	{  
31.	    return (double)rand()/RAND_MAX * (x2-x1)+x1;  
32.	}  
           

并行代码

1.	#include <stdio.h>  
2.	#include <stdlib.h>  
3.	#include <time.h>  
4.	#include "mpi.h"  
5.	#define N 100000000  
6.	  
7.	double rand_v(int x1, int x2);  
8.	int main(int argc, char **argv)  
9.	{  
10.	    MPI_Init(&argc, &argv);  
11.	    MPI_Status status;  
12.	    srand(time(0));  
13.	    int i, n=1000;//每个进程要计算的个数  
14.	    int rank, num_slaves;  
15.	    double sum=0, xr[n], area;  
16.	    enum tag {COMPUTE_TAG,  
17.	                REQ_TAG,  
18.	                STOP_TAG};  
19.	  
20.	    MPI_Comm_rank(MPI_COMM_WORLD, &rank);  
21.	    MPI_Comm_size(MPI_COMM_WORLD, &num_slaves);  
22.	  
23.	    if(!rank){//主进程  
24.	        int j, x2=9, x1=0;  
25.	        double acc, tm = MPI_Wtime();  
26.	  
27.	        for (i=0; i<N/n; i++) {//分成N/n份  
28.	            for (j=0; j<n; j++) {//准备n个随机数据给子进程  
29.	                xr[j] = rand_v(x1, x2);  
30.	            }  
31.	            MPI_Recv(NULL, 0, MPI_DOUBLE,/*等待子进程请求*/  
32.	                    MPI_ANY_SOURCE, REQ_TAG, MPI_COMM_WORLD, &status);  
33.	            MPI_Send(xr, n, MPI_DOUBLE,   
34.	                    status.MPI_SOURCE, COMPUTE_TAG, MPI_COMM_WORLD);  
35.	        }  
36.	  
37.	        for (i=1; i<num_slaves; i++) {  
38.	            MPI_Recv(NULL, 0, MPI_DOUBLE,  
39.	                    i, REQ_TAG, MPI_COMM_WORLD, &status);  
40.	            MPI_Send(NULL, 0, MPI_DOUBLE,  
41.	                    i, STOP_TAG, MPI_COMM_WORLD);//通知子进程停止  
42.	        }  
43.	  
44.	        MPI_Reduce(&sum, &area, 1,   
45.	                MPI_DOUBLE, MPI_SUM, 0,  MPI_COMM_WORLD);//规约求和  
46.	        //sum为输入，area为输出，操作为MPI_SUM，所有进程都要调用此函数  
47.	  
48.	        area = (area/N) * (x2 - x1);  
49.	        tm = MPI_Wtime()-tm;  
50.	  
51.	        acc = (area - 121.5)/121.5;  
52.	        acc = acc>0 ? acc : -acc;  
53.	        printf("area = %f\tacc = %f\ttime = %f\n",   
54.	                area, acc, tm);  
55.	  
56.	    }else{ //子进程  
57.	        MPI_Send(NULL, 0, MPI_DOUBLE,  
58.	                0, REQ_TAG, MPI_COMM_WORLD);//请求主进程  
59.	        MPI_Recv(xr, n, MPI_DOUBLE,   
60.	                0, MPI_ANY_TAG, MPI_COMM_WORLD, &status);  
61.	        while (status.MPI_TAG == COMPUTE_TAG) {  
62.	            for (i  = 0; i < n; i++)  
63.	            {  
64.	                sum = sum + xr[i]*(9-xr[i]);  
65.	            }  
66.	            MPI_Send(NULL, 0, MPI_DOUBLE,  
67.	                    0, REQ_TAG, MPI_COMM_WORLD);  
68.	            MPI_Recv(xr, n, MPI_DOUBLE,   
69.	                    0, MPI_ANY_TAG, MPI_COMM_WORLD, &status);  
70.	        }  
71.	  
72.	        MPI_Reduce(&sum, &area, 1,   
73.	                MPI_DOUBLE, MPI_SUM, 0,  MPI_COMM_WORLD);  
74.	    }  
75.	  
76.	    MPI_Finalize();  
77.	    return 0;  
78.	}  
79.	double rand_v(int x1, int x2)  
80.	{  
81.	    return (double)rand()/RAND_MAX * (x2-x1)+x1;  
82.	}  
           

误差

1. 不同的精度下所求曲线面积和相应的误差

误差	N=100	N=100000	N=100000000
顺序	0.025818	0.000426	0.000021
并行	0.025818	0.000426	0.000015

1. 在不同精度下，对并行代码进行性能分析（执行时间、加速比、计算/通信比）

• 当N=100时：

  

  执行时间：0.00196s
• 当N=100000时：

  

  执行时间：0.004633s
• 当N=100000000时：

  

  执行时间：1.43954s

  加速比：2.4/1.43954=1.667