Monte Carlo integration of sin(x)

46 次查看(过去 30 天)
Jose Morales
Jose Morales 2019-10-29
编辑: James Tursa 2019-10-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 个评论
Fabio Freschi
Fabio Freschi 2019-10-29
Do you really wanted to have 0?
Jose Morales
Jose Morales 2019-10-29
I am trying to do the integration from -1 to 1 for sinx. I want to see how close i can get to 0 since the answer is 0.

请先登录,再进行评论。

回答(2 个)

James Tursa
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 个评论
Jose Morales
Jose Morales 2019-10-30
clear
clc
n=100000;
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_estimation1 = m*sin(1)/n;
for i=1:1:n
x=-rand;
y=-rand*sin(1);
if y-sin(x)<=0
m=m+1;
end
end
format long
sin_estimation2 = m*sin(1)/n;
final_estimation = sin_estimation1+sin_estimation2;
Would the full code be like this then?
James Tursa
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
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')

类别

Help CenterFile Exchange 中查找有关 Preprocessing Data 的更多信息

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by