clc
clear
disp(' HOLZER MATRIX FOR N DOF ');
disp('');
n = input('Enter the number of degrees of freedom (n): ');
if ~isnumeric(n) || n<1 || floor(n) ~= n
error('n must be a integer greater than or equal to 1');
end
start_freq = input('Enter the starting frequency (rad/s): ');
end_freq = input('Enter the ending frequency (rad/s): ');
step_freq = 0.01;
StiffnessK = input(['TYPES OF SYSTEM\1. n DOF, 2 support \n2. n DOF, 1 support \n3.
n DOF, 0 support \n' ...
'Enter the case number (1, 2 or 3): ']);
X1 = double(1);
frequencies = start_freq:step_freq:end_freq;
X_values = zeros(length(frequencies), n+1); % Store X values for plotting
%disp(X_values);
X_values(:, 1) = X1; % Set X1 for all frequencies
natural_frequency = NaN; % To store the natural frequency if found
% Loop through each frequency to calculate X values
for f_idx = 1:length(frequencies)
w = frequencies(f_idx); % Current frequency (rad/s)
X = X1; % Initialize X as X1 for each frequency
switch StiffnessK
case 1
disp('DEFINITE SYSTEM for N DOF');
disp('');
disp('You have chosen for n mass, n+1 spring system');
m = input('Enter an array of masses (size n): ');
k = input('Enter an array of stiffness values (size n+1): ');
for i = 1:n
% Summed Forces for Ratio calculation
for j = 1:n
if j == 1
SumForce = ((m(j) * w*w * X1) - (k(j) * X1));
Ratio = SumForce / k(i);
X_values(f_idx, i+1) = X1 - Ratio;
else
SumForce = SumForce + (m(j) * w*w * X_values(f_idx, j));
Ratio = SumForce / k(i);
X_values(f_idx, i+1) = X_values(f_idx, i) - Ratio; % X3, X4, ...
end
end
end
% Check if we have reached the natural frequency
if abs(X_values(f_idx, i+1)) < 1e-3
natural_frequency = w;
break;
end
% Stop if natural frequency is found
if ~isnan(natural_frequency)
break;
end
if isnan(natural_frequency)
disp('Natural frequency not found in the specified range.');
end
% Plot X vs Frequency
figure;
hold on;
for i = n+1:n+1
plot(frequencies, X_values(:, n+1), 'DisplayName', sprintf('X%d', i));
end
xlabel('Frequency (rad/s)');
ylabel('Displacement (X)');
title('Displacement X vs Frequency');
legend('show');
grid on;
hold off;
case 2
disp('DEFINITE SYSTEM for N DOF');
disp('');
disp(['You have chosen for n mass, n spring system. ' ...
'Please enter k1=0 reference to 2 support system.']);
m = input('Enter an array of masses (size n): ');
k = input('Enter an array of stiffness values (size n+1): ');
for i = 1:n
% Summed Forces for Ratio calculation
for j = 1:n
if j == 1
SumForce = ((m(j) * w*w * X1) - (k(j) * X1));
Ratio = SumForce / k(i);
X_values(f_idx, i+1) = X1 - Ratio;
else
SumForce = SumForce + (m(j) * w*w * X_values(f_idx, j));
Ratio = SumForce / k(i);
X_values(f_idx, i+1) = X_values(f_idx, i) - Ratio; % X3, X4, ...
end
end
end
% Check if we have reached the natural frequency
if abs(X_values(f_idx, i+1)) < 1e-3
natural_frequency = w;
break;
end
% Stop if natural frequency is found
if ~isnan(natural_frequency)
break;
end
if isnan(natural_frequency)
disp('Natural frequency not found in the specified range.');
end
% Plot X vs Frequency
figure;
hold on;
for i = n+1:n+1
plot(frequencies, X_values(:, n+1), 'DisplayName', sprintf('X%d', i));
end
xlabel('Frequency (rad/s)');
ylabel('Displacement (X)');
title('Displacement X vs Frequency');
legend('show');
grid on;
hold off;
case 3
disp('SEMIDEFINITE SYTEM of N DOF');
didsp('');
disp(['You have chosen for n mass, (n-1) spring system. ' ...
'Please enter k1=0 and k(n+1)=0, reference to 2 support system. ']);
m = input('Enter an array of masses (size n): ');
k = input('Enter an array of stiffness values (size n+1): ');
for i = 1:n
% Summed Forces for Ratio calculation
for j = 1:n
if j == 1
SumForce = ((m(j) * w*w * X1) - (k(j) * X1));
Ratio = SumForce / k(i);
X_values(f_idx, i+1) = X1 - Ratio;
else
SumForce = SumForce + (m(j) * w*w * X_values(f_idx, j));
Ratio = SumForce / k(i);
X_values(f_idx, i+1) = X_values(f_idx, i) - Ratio; % X3, X4, ...
end
end
end
% Check if we have reached the natural frequency
if abs(SumForces) < 1e-3
natural_frequency = w;
break;
end
% Stop if natural frequency is found
if ~isnan(natural_frequency)
break;
end
if isnan(natural_frequency)
disp('Natural frequency not found in the specified range.');
end
% Plot SumofForces vs Frequency
figure;
hold on;
for i = n+1:n+1
plot(frequencies, SumForce, 'DisplayName', sprintf('SumOfForces %d', i));
end
xlabel('Frequency (rad/s)');
ylabel('Sum of Forces');
title('Sum of Forces vs Frequency');
legend('show');
grid on;
hold off;
end
end
