[×ÊÔ´] ¡¾×ªÌû¡¿ÄâÃÉÌØ¿¨Âå·½·¨(Quasi-Monte Carlo)»ý·ÖʵÀý

%ʹÓÃMatlabÌṩµÄº¯ÊýÇó»ý·Ö,exp(-1/2*x^2)ÔÚ(0,1)¼ä»ý·Ö format long;

syms x

a = sym(1/2);

f = exp(-a*x^2);

ezplot(f)

disp(int(f,-1,1));

fprintf('integral result:%1.18f.\n',double(int(f,0,1)));

%disp(double(int(f,0,1)));
¸´ÖÆ´úÂë%ʹÓÃÄâÃÉÌØ¿¨Âå·½·¨»ý·Ö
%µÃµ½ÄâÃÉÌØ¿¨ÂåÐòÁÐ,¼´µÍÆ«²îÐòÁÐ,halton·¨
%Èç¹ûÓÐÏà¹ØµÄ¹¤¾ßÏäµÄ»°,¿ÉÒÔÓÃMatlabÀïÃæµÄhaltonset,faureset,sobolsetº¯ÊýʵÏÖ,

x=halton(10000,2,5577);

n=length(x);

mju=0;

for i=1:n

    mju=mju + exp(-0.5*x(i)^2);

end

mju=mju/n;

fprintf('Quasi-Monte Carlo result:%1.18f.\n',mju);

%disp(mju);

%ʹÓÃÃÉÌØ¿¨Âå·½·¨»ý·Ö

%µÃµ½UniformÐòÁÐ,

x=random('unif',0,1,10000,1);

n=length(x);

mju=0;

for i=1:n

    mju=mju + exp(-0.5*x(i)^2);

end

mju=mju/n;

fprintf('Monte Carlo result:%1.18f.\n',mju);



%=============Éú³ÉHALTONÐòÁÐ========================

function result = halton( m,base,seeder )

%Éú³ÉHALTONÐòÁÐ

% Check inputs

if nargin < 3

seeder = 0;

if nargin < 2

      error('MATLAB:Halton:NotEnoughInputs',...

             'Not enough input arguments. See Halton.');

end

end

res=0;

n=length(base);

for i=1:m

   

    for j=1:n

        element=0;

        temp=seeder+i;

        k=1;

        while temp>0

            element(k)=rem(temp,base(j));

            temp=fix(temp/base(j));

            k=k+1;

        end

        res(i,j)= 0;

        for k=1:length(element)

            res(i,j)=res(i,j)+element(k)/(base(j)^k);

        end

    end

   

end

result=res;