I'm simulating the dynamic response of a linear system with uncertainty in one of the parameters of the stiffness matrix. In the code it is the ds(17) parameter of the data set array, which is substituted in the KK stiffness matrix of the dynamic system. All other parameters are constant.
First I calculate a random Gaussian normal distribution of the parameter and plot the result just to see the probability of the value falling into a certain interval. Even if I calculate 1e6 random samples it doesn't take much to compute.
Then I move into calculating the natural frequencies with the mean value of k__p stiffness (as denoted in the KK matrix) or ds(17) for my base system as a reference.
After defining my number of needed random samples num = 1e3; and designate mean value of Km = 1e8;, I move into the first for loop that substitutes the new random ds(17) value into the dynamic matrix KK (which is 18 degrees of freedom) and gives me the new natural frequencies, which I then plot into the next for loop.
With a num = 1e2 my computer takes 30 seconds to give me the results, but 100 samples is way too few for what I need and for representing a Gaussian distribution, I would like to try with at least a thousand but my cpu can't handle, it's currently been running for more than 20 hours for 1e3 but I think I'm gonna stop it.
I have two questions regarding general programming practices: 1. What is the correct way to preallocate acell array of matrices (Kcell and Fcell). 2. Is there a better way of creating an array or list or variables for a data set to input it into a function than what I did (the ds array)?
I tried to leave my code to the minimum for this example. I would appreciate any comments on how to make my code faster, thanks.
% clc;
% clear;
% close all;
%% ~~~~~~~~~~~~~~~~~~~~~~ Simulation Data Array ~~~~~~~~~~~~~~~~~~~~~~~~~ %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
ds(1) = degtorad(24.6);
ds(2) = 0.40;
ds(3) = 0.39;
ds(4) = 0.0774/2;
ds(5) = 2.35;
ds(6) = 3.00;
ds(7) = 0.2750/2;
ds(8) = 5.43;
ds(9) = 6.29;
ds(10) = (0.1762/2)/cos(ds(1));
ds(11) = 0.66;
ds(12) = 0.61;
ds(13) = 0.1003/2;
ds(14) = 1e8;
ds(15) = 1e8;
ds(16) = 1e8;
ds(17) = 1e8;
ds(18) = 0;
ds(19) = 1e9;
ds(20) = 0;
ds(21) = 5e8;
ds(22) = (2*pi/60)*[0];
%% ~~~~~~~~~~~~~~~~~~~~~~~~~ Calculate Matrices ~~~~~~~~~~~~~~~~~~~~~~~~~ %
[q_std,Kcell_std,MM_std,radius_std] = gear_eigen (1,ds);
%% Simulate Model Properties
disp('Standard Eigenproblem Solution w^2*M*phi=K*phi');
% Natural Frequencies(Eigenvalues) & Mode Shapes(Eigenvectors)
[frqs_std] = eigenprob_KM (q_std,cell2mat(Kcell_std),MM_std,ds(22));
%% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ K Uncertainty ~~~~~~~~~~~~~~~~~~~~~~~~~~~ %
num = 1e2;
Kl=6e7; Ku=14e7; % lower and upper limits
Km = 1e8; % Planet bearing stiffness (mean value)
width = 1e5; % Interval width between values
STDEV = 10;
epsilon = normrnd(0,STDEV,num,1)/100;
K_eps = Km*(1+epsilon);
% K_eps = Kl+(Ku-Kl)*rand(num,1); % Method 2: uniform distribution
minK = min(K_eps);
maxK = max(K_eps);
dx = width; % Interval width for empirical PDF
centers = minK : dx : maxK; % Centre of Interval
edges = minK-dx/2 : dx : maxK+dx/2; % Edges of Interval
N = histcounts(K_eps,edges); % Count ocurrences of K
figure ('Name','Random Distribution')
subplot(1,2,1)
histogram(K_eps,centers) % Plot random values vs occurrency
axis([minK-dx maxK+dx 0 1.075*max(N)])
xlabel('Stiffness Value, \itK','FontSize',11)
ylabel('Frequency of Occurrence, \itn_{\rmi}','FontSize',11)
N = N(1:end); % Frequency of Occurrence n
Ni = [0 N 0];
Pri = Ni/(num*dx); %Rescale ordinate values to Pri
xi = minK-dx : dx : maxK+dx;
subplot(1,2,2)
plot(xi,Pri,'.r',xi,Pri,'b')
hold on
axis([minK-dx maxK+dx 0 1.075*max(Pri)])
xlabel('Stiffness Value, \itK','FontSize',11)
ylabel('Approximation of \itf_{\rmK}\rm(\itu\rm)','FontSize',11)
text(0.33,9.2,'MonteCarlo Simulation of Stiffness Uncertainty','FontSize' ,11)
legend('\itf_{\rmi}','Empirical pdf','location','eastoutside')
hold off
% edgeStart = 200;
% edgeEnd = 600;
% eds = edges(edgeStart:edgeEnd);
% % Empirical prob that K will fall between edgeStart and edgeEnd
% Pr = sum(Pri(edgeStart:edgeEnd))*dx;
% T_prob = table(edges(edgeStart),edges(edgeEnd),Pr);
% T_prob.Properties.VariableNames = {'IntervalStart' 'IntervalEnd' ...
% 'Probability'}
% Calculate Natural Frequencies
Fcell{length(K_eps)} = []; % Preallocate memory for cell array
Kcell{length(K_eps)} = []; % Preallocate memory for cell array;
syms x__s y__s u__s ...
x__r y__r u__r ...
x__c y__c u__c ...
zeta__1 eta__1 u__1 ...
zeta__2 eta__2 u__2 ...
zeta__3 eta__3 u__3
q = [x__s y__s u__s ...
x__r y__r u__r ...
x__c y__c u__c ...
zeta__1 eta__1 u__1 ...
zeta__2 eta__2 u__2 ...
zeta__3 eta__3 u__3 ];
alpha = ds(1);
m__s = ds(2);
m__ws = ds(3);
r__s = ds(4);
m__r = ds(5);
m__wr = ds(6);
r__r = ds(7);
m__c = ds(8);
m__wc = ds(9);
r__c = ds(10);
m__p = ds(11);
m__wp = ds(12);
r__p = ds(13);
k__r = ds(14);
k__s = ds(15);
k__c = ds(16);
k__cu = ds(18);
k__ru = ds(19);
k__su = ds(20);
k__m = ds(21);
Omega__c = ds(22);
radius = [ r__s 0 0 r__r 0 0 r__c 0 0 r__p ]; % Array of all radius
MM = [m__c 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 m__c 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 m__wc 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 m__r 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 m__r 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 m__wr 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 m__s 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 m__s 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 m__ws 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 m__p 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 m__p 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 m__wp 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 m__p 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 m__p 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 m__wp 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 m__p 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 m__p 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 m__wp;];
%% ~~ Create Cell Array for K ~~~ %
for i = 1:(length(K_eps))
k__p = K_eps(i); % Random Value
%% N = 3
q = q(1:18);
KK = [-m__c * Omega__c ^ 2 + k__c + 3 * k__p 0 0 0 0 0 0 0 0 -k__p 0 0 k__p / 0.2e1 k__p * sqrt(0.3e1) / 0.2e1 0 k__p / 0.2e1 -k__p * sqrt(0.3e1) / 0.2e1 0; 0 -m__c * Omega__c ^ 2 + k__c + 3 * k__p 0 0 0 0 0 0 0 0 -k__p 0 -k__p * sqrt(0.3e1) / 0.2e1 k__p / 0.2e1 0 k__p * sqrt(0.3e1) / 0.2e1 k__p / 0.2e1 0; 0 0 3 * k__p + k__cu 0 0 0 0 0 0 0 -k__p 0 0 -k__p 0 0 -k__p 0; 0 0 0 k__m * cos(pi / 0.6e1 + alpha) ^ 2 - cos(pi / 0.3e1 + alpha) ^ 2 * k__m - cos(alpha) ^ 2 * k__m - (m__r * Omega__c ^ 2) + 0.2e1 * k__m + k__r sin(pi / 0.6e1 + alpha) * cos(pi / 0.6e1 + alpha) * k__m - k__m * sin(pi / 0.3e1 + alpha) * cos(pi / 0.3e1 + alpha) - cos(alpha) * sin(alpha) * k__m -cos(pi / 0.6e1 + alpha) * k__m + sin(pi / 0.3e1 + alpha) * k__m - sin(alpha) * k__m 0 0 0 cos(alpha) ^ 2 * k__m - k__m cos(alpha) * sin(alpha) * k__m sin(alpha) * k__m -cos(pi / 0.6e1 + alpha) * sin(alpha) * k__m cos(pi / 0.6e1 + alpha) * cos(alpha) * k__m cos(pi / 0.6e1 + alpha) * k__m sin(pi / 0.3e1 + alpha) * sin(alpha) * k__m -sin(pi / 0.3e1 + alpha) * cos(alpha) * k__m -sin(pi / 0.3e1 + alpha) * k__m; 0 0 0 sin(pi / 0.6e1 + alpha) * cos(pi / 0.6e1 + alpha) * k__m - k__m * sin(pi / 0.3e1 + alpha) * cos(pi / 0.3e1 + alpha) - cos(alpha) * sin(alpha) * k__m -k__m * cos(pi / 0.6e1 + alpha) ^ 2 + cos(pi / 0.3e1 + alpha) ^ 2 * k__m + cos(alpha) ^ 2 * k__m - (m__r * Omega__c ^ 2) + k__m + k__r -sin(pi / 0.6e1 + alpha) * k__m - cos(pi / 0.3e1 + alpha) * k__m + k__m * cos(alpha) 0 0 0 cos(alpha) * sin(alpha) * k__m -cos(alpha) ^ 2 * k__m -k__m * cos(alpha) -sin(pi / 0.6e1 + alpha) * sin(alpha) * k__m sin(pi / 0.6e1 + alpha) * cos(alpha) * k__m sin(pi / 0.6e1 + alpha) * k__m -cos(pi / 0.3e1 + alpha) * sin(alpha) * k__m cos(pi / 0.3e1 + alpha) * cos(alpha) * k__m cos(pi / 0.3e1 + alpha) * k__m; 0 0 0 -cos(pi / 0.6e1 + alpha) * k__m + sin(pi / 0.3e1 + alpha) * k__m - sin(alpha) * k__m -sin(pi / 0.6e1 + alpha) * k__m - cos(pi / 0.3e1 + alpha) * k__m + k__m * cos(alpha) 0.3e1 * k__m + k__ru 0 0 0 sin(alpha) * k__m -k__m * cos(alpha) -k__m sin(alpha) * k__m -k__m * cos(alpha) -k__m sin(alpha) * k__m -k__m * cos(alpha) -k__m; 0 0 0 0 0 0 k__m * cos(pi / 0.6e1 + alpha) ^ 2 - cos(pi / 0.3e1 + alpha) ^ 2 * k__m - cos(alpha) ^ 2 * k__m - (m__s * Omega__c ^ 2) + 0.2e1 * k__m + k__s -sin(pi / 0.6e1 + alpha) * cos(pi / 0.6e1 + alpha) * k__m + k__m * sin(pi / 0.3e1 + alpha) * cos(pi / 0.3e1 + alpha) + cos(alpha) * sin(alpha) * k__m cos(pi / 0.6e1 + alpha) * k__m - sin(pi / 0.3e1 + alpha) * k__m + sin(alpha) * k__m cos(alpha) ^ 2 * k__m - k__m -cos(alpha) * sin(alpha) * k__m sin(alpha) * k__m sin(pi / 0.3e1 + alpha) * sin(alpha) * k__m sin(pi / 0.3e1 + alpha) * cos(alpha) * k__m -sin(pi / 0.3e1 + alpha) * k__m -cos(pi / 0.6e1 + alpha) * sin(alpha) * k__m -cos(pi / 0.6e1 + alpha) * cos(alpha) * k__m cos(pi / 0.6e1 + alpha) * k__m; 0 0 0 0 0 0 -sin(pi / 0.6e1 + alpha) * cos(pi / 0.6e1 + alpha) * k__m + k__m * sin(pi / 0.3e1 + alpha) * cos(pi / 0.3e1 + alpha) + cos(alpha) * sin(alpha) * k__m -k__m * cos(pi / 0.6e1 + alpha) ^ 2 + cos(pi / 0.3e1 + alpha) ^ 2 * k__m + cos(alpha) ^ 2 * k__m - (m__s * Omega__c ^ 2) + k__m + k__s -sin(pi / 0.6e1 + alpha) * k__m - cos(pi / 0.3e1 + alpha) * k__m + k__m * cos(alpha) -cos(alpha) * sin(alpha) * k__m -cos(alpha) ^ 2 * k__m k__m * cos(alpha) cos(pi / 0.3e1 + alpha) * sin(alpha) * k__m cos(pi / 0.3e1 + alpha) * cos(alpha) * k__m -cos(pi / 0.3e1 + alpha) * k__m sin(pi / 0.6e1 + alpha) * sin(alpha) * k__m sin(pi / 0.6e1 + alpha) * cos(alpha) * k__m -sin(pi / 0.6e1 + alpha) * k__m; 0 0 0 0 0 0 cos(pi / 0.6e1 + alpha) * k__m - sin(pi / 0.3e1 + alpha) * k__m + sin(alpha) * k__m -sin(pi / 0.6e1 + alpha) * k__m - cos(pi / 0.3e1 + alpha) * k__m + k__m * cos(alpha) 0.3e1 * k__m + k__su -sin(alpha) * k__m -k__m * cos(alpha) k__m -sin(alpha) * k__m -k__m * cos(alpha) k__m -sin(alpha) * k__m -k__m * cos(alpha) k__m; -k__p 0 0 cos(alpha) ^ 2 * k__m - k__m cos(alpha) * sin(alpha) * k__m sin(alpha) * k__m cos(alpha) ^ 2 * k__m - k__m -cos(alpha) * sin(alpha) * k__m -sin(alpha) * k__m -0.2e1 * cos(alpha) ^ 2 * k__m - (m__p * Omega__c ^ 2) + 0.2e1 * k__m + k__p 0 -0.2e1 * sin(alpha) * k__m 0 0 0 0 0 0; 0 -k__p -k__p cos(alpha) * sin(alpha) * k__m -cos(alpha) ^ 2 * k__m -k__m * cos(alpha) -cos(alpha) * sin(alpha) * k__m -cos(alpha) ^ 2 * k__m -k__m * cos(alpha) 0 0.2e1 * cos(alpha) ^ 2 * k__m - (m__p * Omega__c ^ 2) + k__p 0 0 0 0 0 0 0; 0 0 0 sin(alpha) * k__m -k__m * cos(alpha) -k__m sin(alpha) * k__m k__m * cos(alpha) k__m -0.2e1 * sin(alpha) * k__m 0 0.2e1 * k__m 0 0 0 0 0 0; k__p / 0.2e1 -k__p * sqrt(0.3e1) / 0.2e1 0 -cos(pi / 0.6e1 + alpha) * sin(alpha) * k__m -sin(pi / 0.6e1 + alpha) * sin(alpha) * k__m sin(alpha) * k__m sin(pi / 0.3e1 + alpha) * sin(alpha) * k__m cos(pi / 0.3e1 + alpha) * sin(alpha) * k__m -sin(alpha) * k__m 0 0 0 -0.2e1 * cos(alpha) ^ 2 * k__m - (m__p * Omega__c ^ 2) + 0.2e1 * k__m + k__p 0 -0.2e1 * sin(alpha) * k__m 0 0 0; k__p * sqrt(0.3e1) / 0.2e1 k__p / 0.2e1 -k__p cos(pi / 0.6e1 + alpha) * cos(alpha) * k__m sin(pi / 0.6e1 + alpha) * cos(alpha) * k__m -k__m * cos(alpha) sin(pi / 0.3e1 + alpha) * cos(alpha) * k__m cos(pi / 0.3e1 + alpha) * cos(alpha) * k__m -k__m * cos(alpha) 0 0 0 0 0.2e1 * cos(alpha) ^ 2 * k__m - (m__p * Omega__c ^ 2) + k__p 0 0 0 0; 0 0 0 cos(pi / 0.6e1 + alpha) * k__m sin(pi / 0.6e1 + alpha) * k__m -k__m -sin(pi / 0.3e1 + alpha) * k__m -cos(pi / 0.3e1 + alpha) * k__m k__m 0 0 0 -0.2e1 * sin(alpha) * k__m 0 0.2e1 * k__m 0 0 0; k__p / 0.2e1 k__p * sqrt(0.3e1) / 0.2e1 0 sin(pi / 0.3e1 + alpha) * sin(alpha) * k__m -cos(pi / 0.3e1 + alpha) * sin(alpha) * k__m sin(alpha) * k__m -cos(pi / 0.6e1 + alpha) * sin(alpha) * k__m sin(pi / 0.6e1 + alpha) * sin(alpha) * k__m -sin(alpha) * k__m 0 0 0 0 0 0 -0.2e1 * cos(alpha) ^ 2 * k__m - (m__p * Omega__c ^ 2) + 0.2e1 * k__m + k__p 0 -0.2e1 * sin(alpha) * k__m; -k__p * sqrt(0.3e1) / 0.2e1 k__p / 0.2e1 -k__p -sin(pi / 0.3e1 + alpha) * cos(alpha) * k__m cos(pi / 0.3e1 + alpha) * cos(alpha) * k__m -k__m * cos(alpha) -cos(pi / 0.6e1 + alpha) * cos(alpha) * k__m sin(pi / 0.6e1 + alpha) * cos(alpha) * k__m -k__m * cos(alpha) 0 0 0 0 0 0 0 0.2e1 * cos(alpha) ^ 2 * k__m - (m__p * Omega__c ^ 2) + k__p 0; 0 0 0 -sin(pi / 0.3e1 + alpha) * k__m cos(pi / 0.3e1 + alpha) * k__m -k__m cos(pi / 0.6e1 + alpha) * k__m -sin(pi / 0.6e1 + alpha) * k__m k__m 0 0 0 0 0 0 -0.2e1 * sin(alpha) * k__m 0 0.2e1 * k__m;];
Kcell{i} = KK; % Save each KK in Kcell array
[frqs_KM] = eigenprob_KM (q_std,cell2mat(Kcell(i)),MM,ds(22));
Fcell{i} = real(frqs_KM);
end % for
Fmat = cell2mat(Fcell);
for i = 1 : length(q_std)
figure('Name','Natural Frequency with Uncertainty');
NatFrqs = Fmat (i,:);
s = scatter ( epsilon,NatFrqs', 20,'b','filled');
s.LineWidth = 0.6;
s.MarkerEdgeColor = 'b';
s.MarkerFaceColor = [0 0.5 0.5];
Fmat = cell2mat(Fcell);
xlabel('\bf \epsilon \sim N(\mu,\sigma)','FontSize',14)
ylabel('\bf Natural Frequency, \omega_{n}','FontSize',14)
set(gca,'yaxislocation','right');
title(sprintf('Designated Natural Frequency [Hz] = %f', ...
round(frqs_std(i),2)));
ax = gca;
ax.YAxisLocation = 'origin'; % setting y axis location to origin
grid on
hold on
end % for
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Functions ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%% ~~~~~~~~~~~~~~~ Calculate Eigen-Problem for N Planets ~~~~~~~~~~~~~~~~ %
function [q,Kcell,MM,radius] = gear_eigen (randomvar,data)
%% ~~~~~~~~~~~~~~~~~~~~~~ Define Simulation data ~~~~~~~~~~~~~~~~~~~~~~~~ %
syms x__s y__s u__s ...
x__r y__r u__r ...
x__c y__c u__c ...
zeta__1 eta__1 u__1 ...
zeta__2 eta__2 u__2 ...
zeta__3 eta__3 u__3
q = [x__s y__s u__s ...
x__r y__r u__r ...
x__c y__c u__c ...
zeta__1 eta__1 u__1 ...
zeta__2 eta__2 u__2 ...
zeta__3 eta__3 u__3 ];
alpha = data(1);
m__s = data(2);
m__ws = data(3);
r__s = data(4);
m__r = data(5);
m__wr = data(6);
r__r = data(7);
m__c = data(8);
m__wc = data(9);
r__c = data(10);
m__p = data(11);
m__wp = data(12);
r__p = data(13);
k__r = data(14);
k__s = data(15);
k__c = data(16);
k__p = data(17);
k__cu = data(18);
k__ru = data(19);
k__su = data(20);
k__m = data(21);
Omega__c = data(22);
radius = [ r__s 0 0 r__r 0 0 r__c 0 0 r__p ]; % Array of all radius
Kcell{length(randomvar)} = []; % Preallocate memory for cell array;
%% ~~ Create Cell Array for K ~~~ %
for i = 1:(randomvar)
%% N = 3
q = q(1:18);
KK = [-m__c * Omega__c ^ 2 + k__c + 3 * k__p 0 0 0 0 0 0 0 0 -k__p 0 0 k__p / 0.2e1 k__p * sqrt(0.3e1) / 0.2e1 0 k__p / 0.2e1 -k__p * sqrt(0.3e1) / 0.2e1 0; 0 -m__c * Omega__c ^ 2 + k__c + 3 * k__p 0 0 0 0 0 0 0 0 -k__p 0 -k__p * sqrt(0.3e1) / 0.2e1 k__p / 0.2e1 0 k__p * sqrt(0.3e1) / 0.2e1 k__p / 0.2e1 0; 0 0 3 * k__p + k__cu 0 0 0 0 0 0 0 -k__p 0 0 -k__p 0 0 -k__p 0; 0 0 0 k__m * cos(pi / 0.6e1 + alpha) ^ 2 - cos(pi / 0.3e1 + alpha) ^ 2 * k__m - cos(alpha) ^ 2 * k__m - (m__r * Omega__c ^ 2) + 0.2e1 * k__m + k__r sin(pi / 0.6e1 + alpha) * cos(pi / 0.6e1 + alpha) * k__m - k__m * sin(pi / 0.3e1 + alpha) * cos(pi / 0.3e1 + alpha) - cos(alpha) * sin(alpha) * k__m -cos(pi / 0.6e1 + alpha) * k__m + sin(pi / 0.3e1 + alpha) * k__m - sin(alpha) * k__m 0 0 0 cos(alpha) ^ 2 * k__m - k__m cos(alpha) * sin(alpha) * k__m sin(alpha) * k__m -cos(pi / 0.6e1 + alpha) * sin(alpha) * k__m cos(pi / 0.6e1 + alpha) * cos(alpha) * k__m cos(pi / 0.6e1 + alpha) * k__m sin(pi / 0.3e1 + alpha) * sin(alpha) * k__m -sin(pi / 0.3e1 + alpha) * cos(alpha) * k__m -sin(pi / 0.3e1 + alpha) * k__m; 0 0 0 sin(pi / 0.6e1 + alpha) * cos(pi / 0.6e1 + alpha) * k__m - k__m * sin(pi / 0.3e1 + alpha) * cos(pi / 0.3e1 + alpha) - cos(alpha) * sin(alpha) * k__m -k__m * cos(pi / 0.6e1 + alpha) ^ 2 + cos(pi / 0.3e1 + alpha) ^ 2 * k__m + cos(alpha) ^ 2 * k__m - (m__r * Omega__c ^ 2) + k__m + k__r -sin(pi / 0.6e1 + alpha) * k__m - cos(pi / 0.3e1 + alpha) * k__m + k__m * cos(alpha) 0 0 0 cos(alpha) * sin(alpha) * k__m -cos(alpha) ^ 2 * k__m -k__m * cos(alpha) -sin(pi / 0.6e1 + alpha) * sin(alpha) * k__m sin(pi / 0.6e1 + alpha) * cos(alpha) * k__m sin(pi / 0.6e1 + alpha) * k__m -cos(pi / 0.3e1 + alpha) * sin(alpha) * k__m cos(pi / 0.3e1 + alpha) * cos(alpha) * k__m cos(pi / 0.3e1 + alpha) * k__m; 0 0 0 -cos(pi / 0.6e1 + alpha) * k__m + sin(pi / 0.3e1 + alpha) * k__m - sin(alpha) * k__m -sin(pi / 0.6e1 + alpha) * k__m - cos(pi / 0.3e1 + alpha) * k__m + k__m * cos(alpha) 0.3e1 * k__m + k__ru 0 0 0 sin(alpha) * k__m -k__m * cos(alpha) -k__m sin(alpha) * k__m -k__m * cos(alpha) -k__m sin(alpha) * k__m -k__m * cos(alpha) -k__m; 0 0 0 0 0 0 k__m * cos(pi / 0.6e1 + alpha) ^ 2 - cos(pi / 0.3e1 + alpha) ^ 2 * k__m - cos(alpha) ^ 2 * k__m - (m__s * Omega__c ^ 2) + 0.2e1 * k__m + k__s -sin(pi / 0.6e1 + alpha) * cos(pi / 0.6e1 + alpha) * k__m + k__m * sin(pi / 0.3e1 + alpha) * cos(pi / 0.3e1 + alpha) + cos(alpha) * sin(alpha) * k__m cos(pi / 0.6e1 + alpha) * k__m - sin(pi / 0.3e1 + alpha) * k__m + sin(alpha) * k__m cos(alpha) ^ 2 * k__m - k__m -cos(alpha) * sin(alpha) * k__m sin(alpha) * k__m sin(pi / 0.3e1 + alpha) * sin(alpha) * k__m sin(pi / 0.3e1 + alpha) * cos(alpha) * k__m -sin(pi / 0.3e1 + alpha) * k__m -cos(pi / 0.6e1 + alpha) * sin(alpha) * k__m -cos(pi / 0.6e1 + alpha) * cos(alpha) * k__m cos(pi / 0.6e1 + alpha) * k__m; 0 0 0 0 0 0 -sin(pi / 0.6e1 + alpha) * cos(pi / 0.6e1 + alpha) * k__m + k__m * sin(pi / 0.3e1 + alpha) * cos(pi / 0.3e1 + alpha) + cos(alpha) * sin(alpha) * k__m -k__m * cos(pi / 0.6e1 + alpha) ^ 2 + cos(pi / 0.3e1 + alpha) ^ 2 * k__m + cos(alpha) ^ 2 * k__m - (m__s * Omega__c ^ 2) + k__m + k__s -sin(pi / 0.6e1 + alpha) * k__m - cos(pi / 0.3e1 + alpha) * k__m + k__m * cos(alpha) -cos(alpha) * sin(alpha) * k__m -cos(alpha) ^ 2 * k__m k__m * cos(alpha) cos(pi / 0.3e1 + alpha) * sin(alpha) * k__m cos(pi / 0.3e1 + alpha) * cos(alpha) * k__m -cos(pi / 0.3e1 + alpha) * k__m sin(pi / 0.6e1 + alpha) * sin(alpha) * k__m sin(pi / 0.6e1 + alpha) * cos(alpha) * k__m -sin(pi / 0.6e1 + alpha) * k__m; 0 0 0 0 0 0 cos(pi / 0.6e1 + alpha) * k__m - sin(pi / 0.3e1 + alpha) * k__m + sin(alpha) * k__m -sin(pi / 0.6e1 + alpha) * k__m - cos(pi / 0.3e1 + alpha) * k__m + k__m * cos(alpha) 0.3e1 * k__m + k__su -sin(alpha) * k__m -k__m * cos(alpha) k__m -sin(alpha) * k__m -k__m * cos(alpha) k__m -sin(alpha) * k__m -k__m * cos(alpha) k__m; -k__p 0 0 cos(alpha) ^ 2 * k__m - k__m cos(alpha) * sin(alpha) * k__m sin(alpha) * k__m cos(alpha) ^ 2 * k__m - k__m -cos(alpha) * sin(alpha) * k__m -sin(alpha) * k__m -0.2e1 * cos(alpha) ^ 2 * k__m - (m__p * Omega__c ^ 2) + 0.2e1 * k__m + k__p 0 -0.2e1 * sin(alpha) * k__m 0 0 0 0 0 0; 0 -k__p -k__p cos(alpha) * sin(alpha) * k__m -cos(alpha) ^ 2 * k__m -k__m * cos(alpha) -cos(alpha) * sin(alpha) * k__m -cos(alpha) ^ 2 * k__m -k__m * cos(alpha) 0 0.2e1 * cos(alpha) ^ 2 * k__m - (m__p * Omega__c ^ 2) + k__p 0 0 0 0 0 0 0; 0 0 0 sin(alpha) * k__m -k__m * cos(alpha) -k__m sin(alpha) * k__m k__m * cos(alpha) k__m -0.2e1 * sin(alpha) * k__m 0 0.2e1 * k__m 0 0 0 0 0 0; k__p / 0.2e1 -k__p * sqrt(0.3e1) / 0.2e1 0 -cos(pi / 0.6e1 + alpha) * sin(alpha) * k__m -sin(pi / 0.6e1 + alpha) * sin(alpha) * k__m sin(alpha) * k__m sin(pi / 0.3e1 + alpha) * sin(alpha) * k__m cos(pi / 0.3e1 + alpha) * sin(alpha) * k__m -sin(alpha) * k__m 0 0 0 -0.2e1 * cos(alpha) ^ 2 * k__m - (m__p * Omega__c ^ 2) + 0.2e1 * k__m + k__p 0 -0.2e1 * sin(alpha) * k__m 0 0 0; k__p * sqrt(0.3e1) / 0.2e1 k__p / 0.2e1 -k__p cos(pi / 0.6e1 + alpha) * cos(alpha) * k__m sin(pi / 0.6e1 + alpha) * cos(alpha) * k__m -k__m * cos(alpha) sin(pi / 0.3e1 + alpha) * cos(alpha) * k__m cos(pi / 0.3e1 + alpha) * cos(alpha) * k__m -k__m * cos(alpha) 0 0 0 0 0.2e1 * cos(alpha) ^ 2 * k__m - (m__p * Omega__c ^ 2) + k__p 0 0 0 0; 0 0 0 cos(pi / 0.6e1 + alpha) * k__m sin(pi / 0.6e1 + alpha) * k__m -k__m -sin(pi / 0.3e1 + alpha) * k__m -cos(pi / 0.3e1 + alpha) * k__m k__m 0 0 0 -0.2e1 * sin(alpha) * k__m 0 0.2e1 * k__m 0 0 0; k__p / 0.2e1 k__p * sqrt(0.3e1) / 0.2e1 0 sin(pi / 0.3e1 + alpha) * sin(alpha) * k__m -cos(pi / 0.3e1 + alpha) * sin(alpha) * k__m sin(alpha) * k__m -cos(pi / 0.6e1 + alpha) * sin(alpha) * k__m sin(pi / 0.6e1 + alpha) * sin(alpha) * k__m -sin(alpha) * k__m 0 0 0 0 0 0 -0.2e1 * cos(alpha) ^ 2 * k__m - (m__p * Omega__c ^ 2) + 0.2e1 * k__m + k__p 0 -0.2e1 * sin(alpha) * k__m; -k__p * sqrt(0.3e1) / 0.2e1 k__p / 0.2e1 -k__p -sin(pi / 0.3e1 + alpha) * cos(alpha) * k__m cos(pi / 0.3e1 + alpha) * cos(alpha) * k__m -k__m * cos(alpha) -cos(pi / 0.6e1 + alpha) * cos(alpha) * k__m sin(pi / 0.6e1 + alpha) * cos(alpha) * k__m -k__m * cos(alpha) 0 0 0 0 0 0 0 0.2e1 * cos(alpha) ^ 2 * k__m - (m__p * Omega__c ^ 2) + k__p 0; 0 0 0 -sin(pi / 0.3e1 + alpha) * k__m cos(pi / 0.3e1 + alpha) * k__m -k__m cos(pi / 0.6e1 + alpha) * k__m -sin(pi / 0.6e1 + alpha) * k__m k__m 0 0 0 0 0 0 -0.2e1 * sin(alpha) * k__m 0 0.2e1 * k__m;];
MM = [m__c 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 m__c 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 m__wc 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 m__r 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 m__r 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 m__wr 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 m__s 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 m__s 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 m__ws 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 m__p 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 m__p 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 m__wp 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 m__p 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 m__p 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 m__wp 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 m__p 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 m__p 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 m__wp;];
Kcell{i} = KK; % Save each KK in Kcell array
end % for
end % function
%% ~~~~~ Function for Calculating Natural Frequencies & Modes (MDoF) ~~~~ %
%~~~~~~~~~~~~~~~~~~~~~~~ WITH 'K' AND 'M' MATRICES ~~~~~~~~~~~~~~~~~~~~~~~%
function [natural_freq] = eigenprob_KM(q,K,M,w)
%% Compute Eigenvectors 'V' and Eigenvalues 'D'
[V,D] = eig(K,M);
D = diag(sqrtm(D));
frq = (D);
[frq_sorted,index] = sort(frq,'ascend'); % Sort frq values
mode_shape_sorted = V(:,index); % Mode Shapes for each Natural Frequency
natural_freq = frq_sorted/(2*pi); % Natural Frequencies for each DoF
end % function
