|
马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?注册
x
本帖转载于Matlab技术论坛,原帖参见http://www.matlabsky.com/thread-673-1-1.html
%使用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; |
|