MATLAB蒙特卡洛（Monte Carlo）方法求椭圆面积

在某个规定的范围内随机打点，找到满足条件的点，并数一下这些点的数量与总的随机点数量的比，就OK了。关键是设置条件。

代码

clear;clc;
n=10000; %随机数的个数

a=51/2; %长半轴
b=29/2; %短半轴
%生成随机数
x=rand(1,n);
y=rand(1,n);
%改变随机数的范围
x=2*a.*x-a;
y=2*b.*y-b;

r=(1/(a*a)).*(x.*x)+(1/(b*b)).*(y.*y); %椭圆边界
m1=find(r<=1);%满足条件（椭圆内）的点

X=x(m1);Y=y(m1);
mm=length(m1); 
area=a*b*mm/n;%椭圆面积

%把椭圆点出来
figure(1);hold on;
plot(X,Y,'.b');
xlim([-30 30]);ylim([-15 15]);

ax=gca;
fs=18;
ax.XColor = 'red';
ax.YColor = 'red';
ax.XAxisLocation='Origin';
ax.YAxisLocation='Origin';
ax.FontSize = fs;

%把椭圆边缘画出来
fplot(@(x)sqrt((1-x.^2./650.25).*210.25),[-25.5 25.5],...
    'Color','red');

fplot(@(x)-sqrt((1-x.^2./650.25).*210.25),[-25.5 25.5],...
    'Color','red');

如果用一条直线把椭圆切割一下。。。

代码

%Monte Carlo 方法计算椭圆被一条直线切割后其中一部分的面积

clear;clc;
n=50000;

a=51/2; %长半轴
b=29/2; %短半轴
x=rand(1,n);
y=rand(1,n);
x=2*a.*x-a;
y=2*b.*y-b;

r1=(1/(a*a)).*(x.*x)+(1/(b*b)).*(y.*y);
r2=2.*x+5-y;
m1=find(r1<=1);%第一个条件
m2=find(r2>=0);%第二个条件
m3=find(ismember(m2,m1));
M=m2(m3);
X=x(M);
Y=y(M);
mm=length(M); 
area=a*b*mm/n;

figure(1);hold on;
plot(X,Y,'.b');
xlim([-30 30]);
ylim([-15 15]);

ax=gca;
fs=18;
ax.XColor = 'red';
ax.YColor = 'red';
ax.XAxisLocation='Origin';
ax.YAxisLocation='Origin';
ax.FontSize = fs;

% 把椭圆的边缘画出来
fplot(@(x)sqrt((1-x.^2./650.25).*210.25),[-25.5 25.5],...
    'Color','red');

fplot(@(x)-sqrt((1-x.^2./650.25).*210.25),[-25.5 25.5],...
    'Color','red');

%把直线画出来
fplot(@(x)2.*x+5,[-30 30],...
    'Color','blue');

结果如下图：