用蒙特卡罗方法计算定积分（随机投点法）matlab实现

蒙特·卡罗方法（Monte Carlo method），也称统计模拟方法，是二十世纪四十年代中期由于科学技术的发展和电子计算机的发明，而被提出的一种以概率统计理论为指导的一类非常重要的数值计算方法。

算法解释：

设     求在区间    上的积分值

                                                       

设二维随机变量服从正方形上的均匀分布，可知服从均匀分布，且相互独立，又记事件

                                                    

根据边际分布函数推导发生概率等于积分值（积分顺序交换）：

                               

 根据如上推导，定积分的值即为事件发生概率(见图1)。而事件是一个二点分布，既每一次实验都只有1,0两种情况，且每次实验独立，根据伯努利大数定律我们可以用重复试验出现的频率作为的估计值，这种求积分的方法也叫随机投点法。既：

                                             

 

 

                                                            图1  随机投点法

 

 当被积函数值域、积分区间不在上时，可通过对值域、积分区间做线性变换，即可转化为蒙卡特罗方法的领域之内，代码中包含线性转换，可通用于定积分计算。

代码如下：

%%蒙特卡罗方法计算定积分（随机投点法）
%%k累计数、x随机点、y随机点代表一组随机点(x(i)，y(i)),z为转化后的被积函数

function [jifen]=kj(a,b,g,q)     %%输入(a,b)积分区域，g为t的匿名函数,q实验次数
k=0;                              %%累计y小于积分函数的次数
x=rand(1,q);                      %%生成随机数x，y组成(x(i)，y(i))
y=rand(1,q);
t=x*(b-a)+a;                      %%积分区域与（0,1）区间线性转换

for i=1:1:length(x)
    g1(i)=g(t(i));
end

s=[t;g1];                         %% 得g(t)
c=min(s(2,:));                    %%得出被积函数在积分区域的上下限
d=max(s(2,:));

for i=1:1:q                      %%转化函数,线性转化将g（x）函数转化为指定定义域，值域的函数
    m=find(s(1,:)==t(i));        %%解决问题一
    if length(m)>2               %%解决问题二
        m=m(1);
        y1(i)=(s(2,m)-c)*(1/(d-c));
    else
        y1(i)=(s(2,m)-c)*(1/(d-c));
    end
end

for i=1:1:q                        %%转化后的函数进行卡蒙特罗随机试验
    if y(i)<y1(i)
        k=k+1;
    end
end
y2=k/length(x);                     %%卡蒙特罗实验结果
jifen=(b-a)*(d-c)*y2+c*(b-a)        %%g(X)逆向推理，根据积分性质得到积分值
end



%%三个问题1、x=y时，2、两个x相等时，3、给出确定的最大值最小值,4、形参为输入函数

 使用方法：输入积分上限b，积分下限a，投点次数q，变量为t的匿名函数g

举例：求积分             

>> a=2;
>> b=3;
>> g=@(t) t.^3+1/2*t.^2+5*t;
>> q=10000;
>> [jifen]=kj(a,b,g,q);

jifen =

   32.0070

 求解完毕！

 

如对线性转换具体操作方面感兴趣，欢迎留言。