Hello everyone i have this code, it generates 1 plot , now i want to have 10000 samples of it of different values by running a loop and then i need to find the pdf of it, how can i do that do that? please help me, thank you

1 view (last 30 days)
Arpit Sola
Arpit Sola on 4 Sep 2020
Answered: Ayush Gupta on 7 Sep 2020
clear; clc;
n_s = 1e4;
Y2S = @ (x) x * 365 * 24 * 60 * 60; % Transform year into seconds
T_year = 40; % in years
% Generate onset and formulation time
pd_t_form = makedist ('Normal', 'mu', 9.3, 'sigma', 3.162);
t_form_year = random (pd_t_form); % Generate crack formulation time: in year
t_form = Y2S (t_form_year); % in seconds
% Tube parameters
pd_t_p = makedist ('Normal', 'mu', 8.3, 'sigma', 1/3); % Pressure difference: Mpa
P_diff = random (pd_t_p);
pd_t_d = makedist ('Normal', 'mu', 22.23, 'sigma',. 5/3); % Outer diameter: mm
d = random (pd_t_d);
pd_t_t = makedist ('Normal', 'mu', 1.27, 'sigma',. 125 * 1.27 / 3); % inner diameter mm
thickness = random (pd_t_t);
% S_u = 713; % Ultimate tensile strength, MPa
% S_y = 362; % Yeild strength, MPa
T = Y2S (T_year);
% Scotte model parameters: the velocity unit is mm / s, for non-cold workded
% Alloy 600
alpha = 2.8e-9;
K_th = 9;
beta = 1.16;
F = .93;
a_0 = .1; % Initial crack length: mm
a_rup = 62.5; % threshold value mm
t_current = t_form;
sigma_stress = P_diff * d / 2 / thickness; % Stress factor
% Initial values ​​for recording the crack evolution process
a_k = a_0; % a_k: current crack length
t = t_current;
a = a_k;
time_step = 1;
index = 1; % Counter for saving the result
if t <T% If t_current does not exceed $ T $
while t
K = F * sigma_stress * sqrt (pi * a_k / 2); % Stress factor
a_kp = a_k + alpha * (K-K_th) ^ beta * time_step; % Calculate the crack length at the next time point
t_current = t_current + time_step; % Update the time point
a_k = a_kp; % Update current crack length
index = index + 1; % Update counter
time_step = time_step + 1;
if a_kp> a_rup% If the rupture happens
t (index) = t_current;
a (index) = a_kp;
break;
else if t_current> T
t (index) = T;
a (index) = a_kp;
break;
else
t (index) = t_current;
a (index) = a_kp;
end
end
end
else% If t_current exceeds $ T $
t = T; % Censored
end
end
plot (t, a)
xlabel ('Time'), ylabel ('Crack length')

Answers (1)

Ayush Gupta
Ayush Gupta on 7 Sep 2020
The 1000 instances data of the variable can be stored in a global variable which can be later used for calculating a much better and approximate pdf. C is the global array and a is the variable for which we want to calculate the pdf of so after every iteration of the for loop we can add
C = [C a];
After getting C we can use fitdist function and fit a distribution accordingly. After this pdf function in MATLAB can used to find the pdf. Examples on how to use the pdf function can be accessed here.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!