Monte Carlo integration of sin(x)

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;

[James Tursa]

Well, you could use this one-liner since sine is an odd function

sin_estimation = 0;

But I suspect that is not really what you want. Do you really intend to do the monte-carlo integration from -1 to 1 and then see how close to 0 you get?

[Jose Morales]

yes that's what i really inted to do

[Fabio Freschi]

Do you really wanted to have 0?

[Jose Morales]

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.

Answer 1

[James Tursa]

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.

[Jose Morales]

How woud the negitive side look like? Would my recent code be changed all to negative -1 rather than the positive 1 I'm using?

[James Tursa]

You just need to generate your points inside that lower left rectangle. E.g.,

    x = -rand;
    y = -rand*sin(1);

And then for the counting, you would increment (or decrement) m if y was greater than sin(x).

[Jose Morales]

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]

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.

Answer 2

[Fabio Freschi]

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')