montecarlo_area.m

来自「monte carlo方法在各个方面中的小应用」· M 代码 · 共 52 行

%求两曲线间的面积第一种方法
n=input('n=?');%总数
m=input('m=?');%步长
k=1;
j=1;
h=0;%求两曲线间的面积
for i=1:n
    x=3*rand();y=30*rand();
    f1=0.3*x.^4+0.5*x.^3-x.^2-x+2;
    f2=-0.8*x.^2+0.5*x+3;
    if ((y<=f1&&y>=f2)|(y>=f1&&y<=f2))
         h=h+1;
    end;
    if(i==k*m)%记录单位步长的pi值
        aps(j)=m*k;
        k=k+1;
        est(j)=h/aps(j);%面积估计值
        std(j)=90*sqrt(est(j)*(1.-est(j))/aps(j));%标准差
        j=j+1;    
    end    
end
plot(aps,est,'r')
hist(est,800)


%求两曲线间的面积第2种方法
n=input('n=?');%总数
m=input('m=?');%步长
k=1;
j=1;
sum=0;
h=0;%求两曲线间的面积
for i=1:n
    x=3*rand();y=30*rand();
    f1=0.3*x.^4+0.5*x.^3-x.^2-x+2;
    f2=-0.8*x.^2+0.5*x+3;
    if ((y<=f1&&y>=f2)|(y>=f1&&y<=f2))
         h=1;
    else
        h=0;
    end
    sum=sum+h;
    if(i==k*m)%记录单位步长的pi值
        aps(j)=m*k;
        k=k+1;
        est(j)=sum/aps(j);%面积估计值
        std(j)=90*sqrt(est(j)*(1.-est(j))/aps(j));
        j=j+1;    
    end    
end