Monte Carlo integration of sin(x)
37 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
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 Commenti
Risposte (2)
James Tursa
il 30 Ott 2019
Modificato: James Tursa
il 30 Ott 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 Commenti
James Tursa
il 30 Ott 2019
Modificato: James Tursa
il 30 Ott 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.
Fabio Freschi
il 29 Ott 2019
Modificato: Fabio Freschi
il 29 Ott 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')
0 Commenti
Vedere anche
Categorie
Scopri di più su Preprocessing Data in Help Center e File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!