Problem for doing Monte Carlo Integration

12 views (last 30 days)
Deck Zhan Sim
Deck Zhan Sim on 23 Nov 2021
Edited: Deck Zhan Sim on 23 Nov 2021
Hi all, I have a problem for doing the monte carlo integration.
Below here are all the codes, the %% Monte Carlo is the section for monte carlo integration:
%% Monte Carlo
monte_carlo=0; % Initialize Input
for i=1:h
monte_carlo=c*(sum(y)/h(i));
end
toc % Time Execution Stop
hold on
The outputs is looks something like this
I would be grateful that someone can help me to solve this matter. Thanks :)

Accepted Answer

KSSV
KSSV on 23 Nov 2021
Edited: KSSV on 23 Nov 2021
Your y in the ocde is a anonymous function:
function_handle with value:
@(x)sin(x)
You need to input some value into it....Like
y(1), y(2), y(pi) etc...
You cannot use sum(y) straight away.
May be you have to repalce this line:
monte_carlo=c*(sum(y)/h(i));
with
monte_carlo=c*(sum(y(x(1:i)))/h);
  2 Comments
KSSV
KSSV on 23 Nov 2021
But you need to improve the code. If h = 1, it will throw error.
Thanks is accepting/ voting the answer. :)

Sign in to comment.

More Answers (1)

VBBV
VBBV on 23 Nov 2021
clc
clear
close all
%% Initialize Input
tic % Time Execution Start
syms x
a = 0; % Lower Limit
b = pi; % Upper Limit
c = b-a; % Subtraction of Upper Limit and Lower Limit
prompt = 'Enter subinterval,h >>'; % User Input of Interval
h = 2;%input(prompt); % Scan User Input
delta_h = c/h; % Calculating the width of segment
x = [a:delta_h:b]; % Vector
y = @(x) sin(x); % Function
%% Midpoint Rule
midpoint = 0; % Initialize Input
for i = 1:h
midpoint = midpoint+y((2*i+1)/2)*delta_h;
end
toc % Time Execution Stop
Elapsed time is 1.070949 seconds.
%% Trapezoidal Rule
trapezoidal = 0; % Initialize Input
for i = 1:h
trapezoidal = trapezoidal+(y(i)+y(i+1))/2*delta_h;
end
toc % Time Execution Stop
Elapsed time is 1.086139 seconds.
%% Simpson's 1/3 Rule
if h > 1 % The vector is full number form
N = h-1;
simpson = delta_h/3*(y(a)+4*sum(y(2:2:end))+2*sum(y(3:2:end-2))+y(end))
else % The vector is not full number form
simpson = delta_h/3*(y(x(a))+2*sum(y(x(3:2:end-2)))+4*sum(y(x(2:2:end)))+y(x(end)))
end
simpson = 0.4406
toc % Time Execution Stop
Elapsed time is 1.273722 seconds.
%% Monte Carlo
monte_carlo=0; % Initialize Input
for i=1:h
monte_carlo=c*(sum(y(1:i))/h);
end
toc % Time Execution Stop
Elapsed time is 1.280316 seconds.
%% Outputs
fprintf('Midpoint Rule approximation = %f.\n',midpoint);
Midpoint Rule approximation = 2.506939.
fprintf('Trapezoidal Rule approximation = %f.\n',trapezoidal);
Trapezoidal Rule approximation = 2.200046.
fprintf("Simpson's 1/3 Rule approximation = %f.\n",simpson);
Simpson's 1/3 Rule approximation = 0.440593.
fprintf('Monte Carlo Rule approximation = %f.\n',monte_carlo);
Monte Carlo Rule approximation = 2.750101.

Products


Release

R2021a

Community Treasure Hunt

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

Start Hunting!