How to terminate the loop

I want to terminate the loop. Condition is satisfying at i = 96, but loop is still running till the end.

after the termination, i want to know the value of i.

point_load(1) = 0;
for i = 1 : 100
    point_load(i + 1) = point_load(i) + 1;
    PL_BM(i) = point_load(i) * beam_length / 4;
    PL_BM = repelem(PL_BM(i),n);
    failure_domain_BM = ( RV_BM - PL_BM(i) ) < 0;
    p_f_BM = sum(failure_domain_BM (:) == 1) / n; 
    if p_f_BM == 1
        break
        disp(point_load(i+1))
    end
end

Answer 1

[the cyclist]

Due to floating-point error, it may be that the condition

if p_f_BM == 1

is not exactly met. With floating-point numbers, it's better to check for equality within some tolerance. Try something like

if abs(p_f_BM-1) < tol

for some appropriately small value of tol.

[Jaydeep Kansara]

Thank you for the reply. but, still its running till the end.

[Jaydeep Kansara]

This is the entire code, if it helps.

b18 = 3+0+0+7+6+0+4+9; % b18 = 29

b58 = 6+0+4+9; % b58 = 19

n = 100000000; % No of samples = 10,000,000

rng default

beam_length = 5; % Beam length = 5 m

% Bending Moment Capacity follows a log-normal distribution

mean_BM = 2 * b18; % mean of BM = 58 kNm

COV_BM = 0.01 * b58; % COV of BM = 0.19

variance_BM = (mean_BM * COV_BM)^2; % Variance of BM = 121.4404

mu_BM = log(mean_BM^2 / sqrt(mean_BM^2 + variance_BM)); % Location parameter of BM = 4.0427

sigma_BM = sqrt(log(1 + (variance_BM / mean_BM^2))); % Scale paramaeter of BM = 0.1883

RV_BM = lognrnd(mu_BM, sigma_BM, 1, n); % Random variables of BM

% Load Test

% Apply 50 kN load gradually in 1 kN load increment

point_load(1) = 0;

breaking_point = 10^-10;

for i = 1 : 100

point_load(i + 1) = point_load(i) + 1;

PL_BM(i) = point_load(i) * beam_length / 4;

PL_BM = repelem(PL_BM(i),n);

failure_domain_BM = ( RV_BM - PL_BM(i) ) < 0;

p_f_BM = sum(failure_domain_BM (:) == 1) / n;

if abs(p_f_BM - 1) < breaking_point

break

end

end

[Ridwan Alam]

n is too big

you need a larger breaking point, maybe 1e-5?

[Jaydeep Kansara]

Hello Mr. Cyclist,

Thanks for your guidance. I got my answer based on your solution.

Answer 2

[Ridwan Alam]

b18 = 3+0+0+7+6+0+4+9; % b18 = 29
b58 = 6+0+4+9; % b58 = 19
n = 10000000; % No of samples = 10,000,000
rng default
beam_length = 5; % Beam length = 5 m
% Bending Moment Capacity follows a log-normal distribution
mean_BM = 2 * b18; % mean of BM = 58 kNm
COV_BM = 0.01 * b58; % COV of BM = 0.19
variance_BM = (mean_BM * COV_BM)^2; % Variance of BM = 121.4404
mu_BM = log(mean_BM^2 / sqrt(mean_BM^2 + variance_BM)); % Location parameter of BM =  4.0427
sigma_BM = sqrt(log(1 + (variance_BM / mean_BM^2))); % Scale paramaeter of BM = 0.1883
RV_BM = lognrnd(mu_BM, sigma_BM, 1, n); 

point_load(1) = 0;
for i = 1 : 100
    point_load(i + 1) = point_load(i) + 1;
    PL_BM(i) = point_load(i) * beam_length / 4;
    %PL_BM = repelem(PL_BM(i),n);
    failure_domain_BM = sum( RV_BM < PL_BM(i) );
    p_f_BM = failure_domain_BM / n 
    if  p_f_BM >= 0.9999
        disp(['i=' num2str(i)]) 
        break
    end
end

this gave me i=93

[Jaydeep Kansara]

Thanks for the reply.

actually i want to break it at 1.00.

it is breaking before that.

[Ridwan Alam]

If you increase the end limit of i to 120 (i=1:120), you will see that the break point is reached at i = 118. Because the display is rounding up the floating points, you thought the condition is met earlier, whereas it is not. Hope you understand.

b18 = 3+0+0+7+6+0+4+9; % b18 = 29
b58 = 6+0+4+9; % b58 = 19
n = 10000000; % No of samples = 10,000,000
rng default
beam_length = 5; % Beam length = 5 m
% Bending Moment Capacity follows a log-normal distribution
mean_BM = 2 * b18; % mean of BM = 58 kNm
COV_BM = 0.01 * b58; % COV of BM = 0.19
variance_BM = (mean_BM * COV_BM)^2; % Variance of BM = 121.4404
mu_BM = log(mean_BM^2 / sqrt(mean_BM^2 + variance_BM)); % Location parameter of BM =  4.0427
sigma_BM = sqrt(log(1 + (variance_BM / mean_BM^2))); % Scale paramaeter of BM = 0.1883
RV_BM = lognrnd(mu_BM, sigma_BM, 1, n); 
point_load(1) = 0;
for i = 1 : 120
    point_load(i + 1) = point_load(i) + 1;
    PL_BM(i) = point_load(i) * beam_length / 4;
    %PL_BM = repelem(PL_BM(i),n);
    failure_domain_BM = sum( RV_BM < PL_BM(i) );
    p_f_BM = failure_domain_BM / n 
    if  p_f_BM >= 1
        disp(['i=' num2str(i)]) 
        break
    end
end

[Ridwan Alam]

Jaydeep, did you get your answer?

[Image Analyst]

He did mark the cyclist's answer as "Accepted".

[Jaydeep Kansara]

Hello Ridwan,

Thabks for the reply. I got the answer. It was a floating point error. When I kept the difference too small around 5e-4, i got the answer. Thanks again.