How to Code Monte carlo simulation?
10 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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 commentaires
Réponses (1)
Walter Roberson
le 30 Avr 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 commentaires
Walter Roberson
le 27 Mai 2022
Modifié(e) : Walter Roberson
le 27 Mai 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
Voir également
Catégories
En savoir plus sur Digital Filter Analysis dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!