这是我很多年前第一篇在scdn上写的博客,只是当时不会贴贴图,所以图都没有,是个不完整的博客。今天突然想起还一直留着这个问题没填,所以重新把图贴上了。了个事。
题目:
一个2×2的正方形,每次在其中随机取两点,这两点的距离小于√2的概率是多少?
分析:
从题目看两个点正方形中的分布应是平均分布。第一个点落入正方形中任意点的概率是F1=1/S正,第二个点与第一个点是独立的,它落入正方形任意点的概率同第一点,既F2=1/S正。 第二点与第一点距离小于√2的概率F1,2=以第一点为圆心,√2为半径的圆与正方形重合的部分的面积S重/S正。
总的概率为 F=∫∫F1×F1,2 S正。
解法:
如图将正方形分成8个三角形,由于是全对称的,只考虑第一点落在彩色三角中的情况。
现在要求出以这个点为圆心,半径为√2的圆与正方形重合的面积。不过由于以不同的点为圆心画圆时,与正方形各边的覆盖情况不同,所以需要区分对待。
通过分析,可以将这个三角再如图划分为5个部分,每个区域内的S重具有相同的形式。圆心在区域1时,这个圆会盖住左边两个顶点,与上下两边相交,不予右边相交。而在区域2时,会覆盖左上、左下、两个顶点并与右边有两个交点,区域3则不能覆盖左上顶点,但与上边有两个交点,区域4则与上边和右边都有两个交点,区域5则与上边和右边都没有交点。
下面写出各部分S重n的表达式:
其中R=√2
以后的工作就是对以上各式在其积分域内积分。这个工作可以利用Matlab来完成。
用Matlab编写.M文件来完成以上工作,M文件如下:
clear all;
r=sqrt(2);
sum=0;
%定义符号变量
syms x y;
%定义各部分面积计算式
s(1)=2*x+1/2*(2-y)*sqrt(2-(2-y)^2)+1/2*y*sqrt(2-y^2)+pi-acos((2-y)/r)-acos(y/r);
s(2)=2*x+1/2*(2-y)*sqrt(2-(2-y)^2)+1/2*y*sqrt(2-y^2)+pi-acos((2-y)/r)-acos(y/r)-2*acos((2-x)/r)+(2-x)*sqrt(2-(2-x)^2);
s(3)=x*y+1/2*x*sqrt(2-x^2)+1/2*y*sqrt(2-y^2)+pi-acos((2-y)/r)-acos(y/r)+asin(x/r)-acos((2-y)/r)+(2-y)*sqrt(2-(2-y)^2);
s(4)=x*y+1/2*x*sqrt(2-x^2)+1/2*y*sqrt(2-y^2)+pi-acos((2-y)/r)-acos(y/r)+asin(x/r)-acos((2-y)/r)-2*acos((2-x)/r) ...
+(2-y)*sqrt(2-(2-y)^2)+(2-x)*sqrt(2-(2-x)^2);
s(5)=x*y+1/2*x*sqrt(2-x^2)+1/2*y*sqrt(2-y^2)+pi+asin(x/r)-acos(y/r);
%定义各部分积分域
area_y=[2-sqrt(2-x^2),1;
2-sqrt(2-x^2),1;
2-r,2-sqrt(2-x^2);
x,2-sqrt(2-x^2);
x,2-r];
area_x=[0,2-r;
2-r,1;
0,2-r;
2-r,1;
0,2-r];
%计算积分
for n=1:5
temp(n)=int(s(n),y);
fx(n)=subs(temp(n),y,area_y(n,2))-subs(temp(n),y,area_y(n,1));
end
for n=1:5
temp(n)=int(fx(n),x);
f(n)=subs(temp(n),x,area_x(n,2))-subs(temp(n),x,area_x(n,1));
sum=sum+f(n);
end
然而很不幸的是,运行时提示有误,经检查发现是第二轮积分有式子得不到解析解,所以必须另想办法。
经过分析,我发现第一轮积分后,将y值按各积分域的值代入,使式子复杂化,所以产生了无法得到解析解的情况,但是实际上有些因子比如1/2× y√(2-y^2)在不止一个区域里出现,而将几个相邻的积分域合并的话,可以简化积分后的结果。
将所有出现的因子和出现的域列表
| 1 | 2 | 3 | 4 | 5 |
2x | Y | Y |
|
|
|
1/2*(2-y)*sqrt(2-(2-y)^2) | Y | Y |
|
|
|
1/2*y*sqrt(2-y^2) | Y | Y | Y | Y | Y |
pi | Y | Y | Y | Y | Y |
-acos(y/r) | Y | Y | Y | Y | Y |
-acos((2-y)/r) | Y | Y | Y | Y |
|
-2*acos((2-x)/r) |
| Y |
| Y |
|
(2-x)*sqrt(2-(2-x)^2) |
| Y |
| Y |
|
x*y |
|
| Y | Y | Y |
1/2*x*sqrt(2-x^2) |
|
| Y | Y | Y |
asin(x/r) |
|
| Y | Y | Y |
-acos((2-y)/r) |
|
| Y | Y |
|
(2-y)*sqrt(2-(2-y)^2) |
|
| Y | Y |
|
| 积分结果 |
|
2x | -5/3+4/3*2^(1/2) |
|
1/2*(2-y)*sqrt(2-(2-y)^2) | 1/8 |
|
1/2*y*sqrt(2-y^2) | pi/16 |
|
pi | pi/2 |
|
-acos(y/r) | -1/4*pi+1/4 |
|
-acos((2-y)/r) | 1/2*pi-7/4 |
|
-2*acos((2-x)/r) | -1/2*pi+3/2 |
|
(2-x)*sqrt(2-(2-x)^2) | -1/3+1/8*pi |
|
x*y | 23/12-4/3*2^(1/2) |
|
1/2*x*sqrt(2-x^2) | -17/24-1/16*pi+2/3*2^(1/2) |
|
asin(x/r) | 0.04124498 | 得不到解析解 |
-acos((2-y)/r) | -0.01347756 | 得不到解析解 |
(2-y)*sqrt(2-(2-y)^2) | 5/12-1/8*pi |
|
和= | 1.50597463 |
|
总概率=和/2= | 0.75298731 |
|
精确值 | 2/3√2+pi/4-1/4+∫∫…dΩ |
|
为便于计算,写了两个M文件来进行积分,这两个文件分别是先积x方向和先积y方向。程序如下:
file name intxy:
fxy=int(f,x);
f1=subs(fxy,x,area_xy(2))-subs(fxy,x,area_xy(1));
f2=simple(f1);
fs=int(f2,y,area_xy(3),area_xy(4));
file name intyx:
fxy=int(f,y);
f1=subs(fxy,y,area_yx(2))-subs(fxy,y,area_yx(1));
f2=simple(f1);
fs=int(f2,x,area_yx(3),area_yx(4));
在实践过程中,发现有时在积分前用 simple 函数对被积函数化简是很重要的,否则会得不到正确解。
对于两个不能得到解析解的因子,只好采用数值积分的方法得到近似解了。所用的函数是trapz(x,y)。x是取值在积分域上的一个向量,所取得间隔是0.0001,y是因子注变量取x时得到的值向量。
也可以用quad(@fx,a,b)指令,fx是一个函数M文件。两个函数文件分别是:
function y=fx(x)
y=(sqrt(2-(2-x).^2)-x).*acos((2-x)./sqrt(2));
function y=fx1(x)
y=(2-(2-x.^2).^(1/2)-x).*asin(x./sqrt(2));
注:这个因子是已经对一个变量积分一次,将该变量积分域代入后得到一个一元函数。
通过模拟,证实了以上结果的正确。模拟程序是随机取了1G对(2^30)的点进行统计,得到的结果非常吻合。模拟程序如下:
clear all;
sum=[0,0];
tic
for n=1:1024
a=rand(2,1024*1024);
b=rand(2,1024*1024);
a=b-a;
b=a(1,:).*a(1,:)+a(2,:).*a(2,:);
bound=[0,1];
y=hist(b,bound);
sum=sum+y;
end
prob=sum./2^30
t=toc
2003-8-8 补充
今天使用Mathematica5.0软件计算用Matlab无法计算出的两个积分,发现可以算出。
约等于0.7529872852