Monte Carlo integration of sin(x)
46 次查看(过去 30 天)
显示 更早的评论
This is a code for the integration of sin(x) from 0 to 1. How do I change it so it can be from -1 to 1?
clear
clc
n=10000;
m=0;
for i=1:1:n
x=rand;
y=rand*sin(1);
if y-sin(x)<=0;
m=m+1;
end
end
format long
sin_estimation = m*sin(1)/n;
4 个评论
回答(2 个)
James Tursa
2019-10-30
编辑:James Tursa
2019-10-30
This is what you are currently doing with the "random counting in an area" method:
To make it go from -1 to +1 instead, you could define two rectangles, one on the positive side and one on the negative side.
You could use your current code for the positive side, and then do something similar for the negative side (generate your x's and y's slightly differently and make your test slightly different) and then add (or subtract depending on how you do the counting) the two values to get the final result. The lower left value will be negative, so either subtract 1 from m at each counting point, or add 1 at each counting point and then negate the total at the end.
4 个评论
James Tursa
2019-10-30
编辑:James Tursa
2019-10-30
No. I wrote the following for the lower left rectangle counting: "... if y was greater than sin(x) ..."
Remember, you are trying to count the points above the sin(x) curve when you are in that lower left rectangle. But here is your coded test:
if y-sin(x)<=0
So your test is backwards.
And then, I also wrote the following: "... The lower left value will be negative ..."
You are using the sin_estimation2 value as positive, which again is backwards.
Fix these two errors.
Fabio Freschi
2019-10-29
编辑:Fabio Freschi
2019-10-29
This code evaluates the integral using the Monte Carlo method with increasing number of random samples, compare the result with exact integration and plots the relative error
% function to integrate
f = @(x)sin(x);
% interval
a = 0;
b = 1;
% numebr of samples
N = logspace(1,8,8);
% preallocation
intMC = zeros(size(N));
% integral value
for i = 1:length(N)
% samples
xS = a+(b-a)*rand(N(i),1);
% Monte Carlo integration
intMC(i) = (b-a)*sum(f(xS))/N(i);
end
% exact integration
fsym = sym(f);
syms x
intEx = vpaintegral(sym(f),x,[a b]);
% relative error
relErr = abs((intMC-intEx)/intEx);
% plot results
figure
loglog(N,relErr)
xlabel('number of samples')
ylabel('relative error')
0 个评论
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Preprocessing Data 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!