Problem in Monte Carlo simulation.

12 views (last 30 days)
harsh Brar
harsh Brar on 6 May 2022
Commented: Walter Roberson on 27 May 2022
My code is not giving me the desired result. It always shows 1. For your ease let me explain the base of code again
The failure occur when "C >= Cr". so I intend to calculate the probabilty of that failure, lets say its P_fail. Here, C will be Calculated using
C(i,1) = Cs(i,1).*(1-erf(Cd(i,1)./(2.*sqrt(con))));
and Cr = 1.2
and ultimately final output is calculated when P_fail becomes greater than P_fail_Max (which equals to 0.10).
clear all
clc
n = 200;
%cover depth
mean_Cd = 16.1;
sigma_Cd = 0.23;
%diffusion coefficient at 28 days
mean_D = 3.87*10^(-12);
sigma_D = 0.52;
%surface chloride concn
mean_Cs = 13.1;
sigma_Cs = 0.103;
%Time exponent
mean_T = 0.2;
sigma_T = 0.2;
%critical chloride content
mean_Cr = 1.2;
sigma_Cr = 0.2;
%Now use normrnd function and lognrnd to generate values randomly
Cd = normrnd(mean_Cd, sigma_Cd, [n,1]);
D = lognrnd(mean_D, sigma_D, [n,1]);
Cs = normrnd(mean_Cs, sigma_Cs, [n,1]);
T = normrnd(mean_T, sigma_T, [n,1]);
Cr = normrnd(mean_Cr, sigma_Cr, [n,1]);
%Create for loop
for i=1:n
rng(0,'twister'); %generating a vector of 200 random values for Cr values
a = 0.52;
b = 3.87*10^(-12);
D = a.*randn(200,1) + b;
%display(D)
rng(0,'twister'); %generating a vector of 200 random values for Cr values
a = 0.2;
b = 0.2;
T = a.*randn(200,1) + b;
% display(T)
if (D(i,1).*T(i,1) > 0 ) % Check to see if this number is real and positive before calculating sqrt
con = D(i,1).*T(i,1);
else
con = 1;
%break;
end
rng(0,'twister'); %generating a vector of 200 random values for Cr values
a = 0.2;
b = 1.2;
Cr = a.*randn(200,1) + b;
%display(Cr)
rng(0,'twister'); %generating a vector of 200 random values for Cr values
a = 0.23;
b = 16.1;
Cd = a.*randn(200,1) + b;
% display(Cd)
rng(0,'twister'); %generating a vector of 1000 random values for Cr values
a = 0.103;
b = 13.1;
Cs = a.*randn(200,1) + b;
%display(Cs)
% MAIN EQUATION TO CALCULATE C
C(i,1) = Cs(i,1).*(1-erf(Cd(i,1)./(2.*sqrt(con))));
a = 1;
b = i;
C = a.*randn(20,1) + b;
display(C)
Pf = mean(C > Cr);
display(Pf)
end
Pfmax = 0.10;
Tser = mean(Pf >= Pfmax);
display(Tser) % final output
  2 Comments
Walter Roberson
Walter Roberson on 27 May 2022
https://www.mathworks.com/matlabcentral/answers/1707460-how-to-code-monte-carlo-simulation#comment_2181870

Sign in to comment.

Answers (0)

Products


Release

R2020a

Community Treasure Hunt

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

Start Hunting!

Translated by