Problem in Monte Carlo simulation.

Question

harsh Brar 2022-5-6

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1713215-problem-in-monte-carlo-simulation

评论： Walter Roberson 2022-5-27

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 