Monte Carlo integration of sin(x)

66 views (last 30 days)
Jose Morales
Jose Morales on 29 Oct 2019
Edited: James Tursa on 30 Oct 2019
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 Comments
Jose Morales
Jose Morales on 29 Oct 2019
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.

Sign in to comment.

Answers (2)

Fabio Freschi
Fabio Freschi on 29 Oct 2019
Edited: Fabio Freschi on 29 Oct 2019
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')

James Tursa
James Tursa on 30 Oct 2019
Edited: James Tursa on 30 Oct 2019
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 Comments
James Tursa
James Tursa on 30 Oct 2019
Edited: James Tursa on 30 Oct 2019
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.

Sign in to comment.

Categories

Find more on Operating on Diagonal Matrices in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by