MATLAB: I’m working on bearing capacity formula described in the model below. Two variables have been choose as a random variable (i.e c and Unit_weight). I want to calculate probability of failure vs Qall to plot a histogram. Please kindly do assist me in

bearing capacity for strip footing

I'm working on bearing capacity formula described in the model below. Two variables have been choose as a random variable (i.e c and Unit_weight). I want to calculate probability of failure vs Qall to plot a histogram. Please kindly do assist me in that regard.
The only condition is if Qallowable is greater than Qultimate, if fails. I have attached the code. Please, kindly do help. If you need any clarification from my side. Thank you
% Soil properties
% Soil type: Cohesive Soil
% c is a random variable
% Unit_weight is a random variable
B = 1;
D = 1;
m1 = 100; 
s1 = 50;
c = m1 + s1.*randn(1000,1);
m2 = 10;
s2 = 5;
Unit_weight = m2 + s2.*randn(1000,1);
Phi = 0;
FOS = 3; 
% From Table
Nc = 5.71;
Nq = 1.0;
Ngamma = 0.0;
Qult = c*Nc + Unit_weight*D*Nq + 0.5*B*Ngamma;
Qall = Qult/3;
% FOS = Qult(:,1)/Qall(:,1);
m1_normal = 0.47;
s1_normal = 6.13;
m2_nomal = 0.20;
s2_normal = 3.23;
% Monte Carlo Reliability Assessment
failure = 0;
for i=1:1000
    if Qult(i)<Qall(i)
        failure = failure +1;
    end
    beta(i) = norm(1-failure/1000);
end
figure(1),
hist(beta);
xlabel("probability of failure")

Best Answer

Don't need loop to compute the failure fraction--
nFail=sum(Qult<Qall);
But, since you've written
Qult = c*Nc + Unit_weight*D*Nq + 0.5*B*Ngamma;
Qall = Qult/3;
the ratio is fixed at 1/3 and Qult<Qall can never be anything but false.
There should be a fixed failure value that is dependent upon the materials and whatever the reference is; not enough info to tell just what you're modeling here for certain. Then there should be an estimated value compared to that, the 3 looks more like a safety factor used for loadings in many cases. Either way, comparing the computed value to 1/3 the same value on a case-by-case basis isn't what you're intending nor needing to do here.
As you've discovered, beta below is constant for all i, so
beta = norm(1-nFail/1000);
ADDENDUM
Follow up to conversation below -- if run your code above stopping a computing the estimate Qult values, then the traditional way of determining percentiles from the sampled distribution is by way of order statistics--
>> P=length(Qult)*0.995;  % index to 99.5% estimate
>> QFail=sort(Qult); QFail=QFail(P)  % select 995th-largest
QFail =
 1.2532e+03
Alternatively, fit a normal distribution to the values and use it to estimate--
>> [mu,sd]=normfit(Qult)  % estimate parameters of normal
mu =
585.2466
sd =
279.7926
>> z=norminv(0.995,mu,sd)  % compute the tail z 
z =
 1.3059e+03
>> ypdf=normpdf(-500:1500,mu,sd);  % calculate a distribution
>> histogram(Qult,20)        % histogram the simulation
>> hold all
>>hpdf=plot(-500:1500,max(n)*ypdf/max(ypdf),'b-','linewidth',2);   % add the fitted line normalized to peak
>> hL995=line([QFail QFail],[0 100],'color','m');  % the two
>> hLZ995=line([z z],[0 100],'color','r');  % point estimates
>>
While the order statistics gives approximately the same value as the z-value of the fitted distribution they're not the same. There's only so close you can get with the sample statistics based simply on the number of trials -- 1-1/1000 is 0.999; anything after that is 1.000; often a 1/2 is used as a continuum correction.
The above results in the following illustration:
Related Question