How to Code Monte carlo simulation?
5 ビュー (過去 30 日間)
古いコメントを表示
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 件のコメント
回答 (1 件)
Walter Roberson
2022 年 4 月 30 日
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 件のコメント
Walter Roberson
2022 年 5 月 27 日
編集済み: Walter Roberson
2022 年 5 月 27 日
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
参考
カテゴリ
Help Center および File Exchange で Digital Filter Analysis についてさらに検索
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!