程序代码：

alpha=4

L=100;
K=50;
R=1;

r1=R*sqrt(rand(1,L));
seta1=2*pi*rand(1,L);
x1=r1.*cos(seta1);
y1=r1.*sin(seta1);
% plot(x1,y1,'rd')


r2=R*sqrt(rand(1,K));
seta2=2*pi*rand(1,K);
x2=r2.*cos(seta2);
y2=r2.*sin(seta2);
% plot(x2,y2,'b*')


gamma1=zeros(K,L);
d=zeros(1,K);
for u=1:K
    for b=1:L
        gamma1(u,b) = sqrt((x2(u)-x1(b))^2+(y2(u)-y1(b))^2);
%         gamma(u,b) = gamma1(u,b)^(-alpha/2);
    end
    d(u) = min(gamma1(u,:));
end

f=zeros(1,K);
syms R  t;
R=@(t) exp(-t)/t;
for k = 1:K
    f(k)=integral(R,d(k),Inf);
end
plot(d,f,'r*');