蒙特卡罗模拟法 —— matlab

目录

1.简介

2.实例分析

2.1 模拟求近似圆周率

2.2 估算定积分

2.3 求解整数规划

---

1.简介

        蒙特卡洛又称随机抽样或统计试验，就是产生随机变量，带入模型算的结果，寻优方面，只要模拟次数够多，最终是可以找到最优解或接近最优的解。

         基本思想：为了解决数学、物理、工程技术等方面的问题，首先建立一个概率模型或随机过程，使它的参数等于问题的解；然后通过对模型或过程的观察或抽样试验来计算所求参数的统计特征，最后给出所求解的近似值。

2.实例分析

2.1 模拟求近似圆周率

        绘制单位圆和外接正方形，正方形ABCD的面积为：2*2=4，圆的面积为：S=Π*1*1=Π，现在模拟产生在正方形ABCD中均匀分布的点n个，如果这n个点中有m个点在该圆内，则圆的面积与正方形ABCD的面积之比可近似为m/n

程序如下：

%%1.模拟求近似圆周率
clc;clear;close all
num=0:10:200000;
mypi=ones(1,length(num));
for j=1:length(num)
    n=num(j);m=0; 
    for i=1:n
        if (-1+2*rand)^2+(-1+2*rand)^2<=1
            m=m+1;
        end
    end
    mypi(j)=4*m/n;
end
plot(mypi)
hold on
line([0 ,length(num)*1.1],[pi,pi],'color','r')
text(0,pi,'\pi','color','r','fontsize',16)
legend('模拟π','实际π')
grid minor

返回：

2.2 估算定积分

程序如下：

%%用蒙特卡罗法估算定积分
clc;clear;close all;
num=0:500:10^6;
s=ones(1,length(num));
for j=1:length(num)
    n=num(j);
    a=0;b=1;
    d=max(a,b)+1;
    m=0;
    for i=1:n
        x=a+rand*(b-a);
        y=d*rand;
        if y<=x^2
            m=m+1;
        end
    end
    s(j)=m/n*d*(b-a);
end
plot(s)
hold on
line([0 ,length(num)*1.1],[1/3,1/3],'color','r')
text(0,pi,'1/3','color','r','fontsize',16)
legend('模拟','实际1/3')
grid minor

返回：

2.3 求解整数规划

要解的方程为：

条件如下：

程序如下：

%%用蒙特卡罗法求解整数规划
clc;clear;close all;
rand('state',sum(clock));
%设置该命令是因为每次产生随机数的时候，随机数生成器触发器的状态都会翻转一次。
%matlab生成的随机数是伪随机数，因此可生成时间相关的随机数，总之和当前时间相关。
%如果计算机运算太快的话，可能会生成相同随机数
p0=0;tic
%计时开始
for i=1:10^7%只要次数够高，最后肯定是实际的最优值    
    x=randi([0,99],1,3);%产生一行五列的区间在【0，99】上的随机整数    
    f=2*x(1)+3*x(1)^2+3*x(2)+x(2)^2+x(3);%主函数    
    g=[x(1)+2*x(1)^2+x(2)+2*x(2)^2+x(3)        
        x(1)+x(1)^2+x(2)+x(2)^2-x(3)        
        2*x(1)+x(1)^2+2*x(2)+x(3)        
        x(1)+2*x(2)];%条件    
    if g(1,1)<=100&&g(2,1)<=500&&g(3,1)<=400&&g(4,1)>=10        
        if p0<f%如果求最小值，则将if p0<f中的<改成>符号            
            x0=x;            
            p0=f;%记录当前较好的解        
        end
    end
end
x0 %x0为x1，x2，x3的值
p0 %p0为函数最大值toc%计时结束

返回：