how should I add to cord
2 views (last 30 days)
Show older comments
syms c m k w w_n c_c zeta x X p;
w_n = sqrt(k/m); % natural frequency of undamped oscillation
c_c = 2*m*w_n; % critical damping
zeta = c/ c_c; % damping factor
p = atan(c*w/(k-m*w^2)); % Φ
X = F_0/(2*zeta*k); % amplitude of oscillation
% the differential equation and its complete solution
% including the transient term
% x" + 2*zeta*x' + (w_n)^2*x = F_0*sin(wt)/m
eqn = F_0/k - sin(w*t - p)/sqrt((1-(w/w_n)^2)+(2*zeta*w/w_n)^2)...
+ X*exp(-zeta*w_n*t)*sin(sqrt(1-zeta^2)*w_n*t+p);
% Solve the differential equation
sol = dsolve(eqn);
% Display the solution
disp('Solution:');
disp(sol);
% Plot the solution
figure;
plot(sol, [0, 10], 'LineWidth', 2);
title('Solution of the Mass-Damper-Spring System with Harmonic Excitation');
xlabel('Time');
ylabel('Displacement');
I've written the code up to this point. I need to draw the attached picture. How should I use the plot?
1 Comment
Dyuman Joshi
on 2 Dec 2023
There are undefined parameters in your code.
You are using dsolve, but have not defined a differential equation. And I don't see a relation between the differential equation in the comments and the equation defined in the code just below.
Answers (1)
Sulaymon Eshkabilov
on 2 Dec 2023
How you are trying to solve your given exercise is not correct.
Here is one function that I wrote to simulate the 1 DOF spring-mass-damper system that is what you're trying to plot.
M=2; S=50;
omegaN=sqrt(S/M);
omega=0:.2:10;
r=omega./omegaN;
zeta=0:0.2:1;
zeta_1DOF_SMD(M, zeta, S, omega)
function zeta_1DOF_SMD(M, zeta, S, omega)
% HELP. Study damping ratio (zeta) influence on Amplitude change
% M - Mass of the system
% D - Damping coefficient of the system ranges (row matrix): 0...
% S - Stiffness of the system
%{
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Written by Sulaymon L. ESHKABILOV, PhD %
% April, 2013 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%}
close all
if nargin<1
M=2.5; S=25; omegaN=sqrt(S/M); omega=0:.2:10; r=omega./omegaN; zeta=0:0.2:1;
end
omegaN=sqrt(S/M);
r=omega/omegaN;
Mag=zeros(length(zeta), length(r));
Theta=ones(length(zeta), length(r));
for ii=1:length(zeta)
for k=1:length(r)
Mag(ii,k)=1/sqrt((1-r(k)^2)^2+(2*zeta(ii)*r(k))^2);
Theta(ii,k)=atan2(2*zeta(ii)*r(k), (1-r(k)^2));
end
end
%
close all
Labelit = {};
Colorit = 'bgrmkgbckmbgrygr';
Lineit = '--:-:--:-:--:----:----:--';
Markit = 'oxs+*^v<p>.xsh+od+*^v';
for m=1:length(zeta)
Stylo = [Colorit(m) Lineit(m) Markit(m)];
Labelit{m} = ['\zeta= ' num2str(zeta(m))];
semilogy(r, Mag(m,:), Stylo, 'linewidth', 1.5), hold on
end
legend(Labelit{:})
title('\zeta influence on Response Amplitude change'),
xlabel('r = \omega/\omega_n'), ylabel('Normalized Amplitude'), axis tight; grid on; ylim([0, 5])
hold off
figure
for m=1:length(zeta)
Stylo = [Colorit(m) Lineit(m) Markit(m)];
Labelit{m} = ['\zeta= ' num2str(zeta(m))];
plot(r, Theta(m,:), Stylo, 'linewidth', 1.5)
hold on
end
legend(Labelit{:}), axis tight, grid on, xlim([0, 2]),
title('\zeta influence on Phase change')
xlabel('r = \omega/\omega_n '), ylabel('Phase, \theta(r)')
end
0 Comments
See Also
Categories
Find more on Calculus in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!