How to Code Monte carlo simulation?
5 views (last 30 days)
Show older comments
In this code I am intending to get service life T(ser)
T(ser) = [ Pf(t) >= Pfmax ]
here , P is probability of an event.
Pfmax is generally between 7 to 10% and
Pf(t) = P[C >= Cr],
C = Cs[1-erf(Cd/(2*sqrt(D*T)))]
Here Cd, D, Cs, T, Cr are variables which are generated randomly, I have specified them in my code.
This is my code which is uncomplete. I am clueless on how to proceed further.
n = 20000;
%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
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]);
C = zeros(n,1);
K = 0;
%Create for loop
for i=1:n
C(i,1) = Cs(i,1)[1-erf(Cd(i,1)/(2*sqrt(D(i,1)*T(i,1))))]
if
C(i,1)>=Cr(i,1)
K = K + 1;
end
end
0 Comments
Answers (1)
Walter Roberson
on 30 Apr 2022
C(i,1) = Cs(i,1)[1-erf(Cd(i,1)/(2*sqrt(D(i,1)*T(i,1))))]
You forgot the multiplication; MATLAB does not have implied multiplication anywhere
Also, in MATLAB, [] always has to do with building lists or arrays. Although it is often possible to use [] as a form of prioritization brackets, it is more efficient to use () instead of [] for prioritization
You should probably be considering vectorization. Generate a vector of Cs values and the corresponding C values, and then if you
Pf = mean(C > Cs)
No need to loop counting the occurances one by one.
8 Comments
Walter Roberson
on 27 May 2022
Edited: Walter Roberson
on 27 May 2022
Cr = a.*randn(200,1) + b;
A column vector of length 200
C(i,1) = Cs(i,1).*(1-erf(Cd(i,1)./(2.*sqrt(con))));
C has not been initialized. In that context you would be creating a column vector.
C = a.*randn(20,1) + b;
That overwrites all of C with a column vector of length 20
Pf = mean(C > Cr);
column vectors of different lengths. Error.
If your actual code makes C the right length then the mean would be a scalar
Tser = mean(Pf >= Pfmax);
After the for loop, Pf is still a scalar. mean of a logical scalar is 0 or 1
mean of logical is a probability, never greater than 1
See Also
Categories
Find more on Digital Filter Analysis 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!