how should I add to cord

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?

[Dyuman Joshi]

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.

Answer 1

[Sulaymon Eshkabilov]

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