MATLAB Answers

Study of oscillations and how to solve in MATLAB

11 views (last 30 days)
Khang Ngo
Khang Ngo on 30 May 2021
Edited: Khang Ngo on 4 Jun 2021
I have known the oscillation of any body due to elastic force can be described by the differential equation:
y''+b*y'+(w0^2)*y =F*cos(w*t)
In which, y is oscillation displacement, b is damped coefficient, w0 is angular frequency of free oscillation, w is angular frequency of stimulating force.
And I have tried for the case of stimulated oscillation where (w0 = 10; b = 0.1 ; F = 10; w= 10.0, 5.0, 3.0, 0.0; t = 150s ) % many values of w, (with initial conditions y(0) = 5; y’(0) = 0)
I have written this one but I couldn't find something is wrong. And I got this figure below.
Can someone help me with the problem
>> syms y(t)
Dy = diff(y, t);
D2y = diff(Dy, t);
w0 = 10; b = 0.1; F = 10; w = [10.0, 5.0, 3.0, 0.0]; t = 150;
SOL1 = dsolve(D2y==F*cos(w(1)*t)-b*Dy-y*w0^2, y(0)==5, Dy(0)==0);
[time1, Y1]=fplot(SOL1, [0, 150]);
plot(time1, Y1, 'k-'), hold on
xlabel('t'), ylabel('y(t)'), grid on
SOL2 = dsolve(D2y==F*cos(w(2)*t)-b*Dy-y*w0^2, y(0)==5, Dy(0)==0);
[time2, Y2]=fplot(SOL2, [0, 150]);
plot(time2, Y2, 'r--')
SOL3 = dsolve(D2y==F*cos(w(3)*t)-b*Dy-y*w0^2, y(0)==5, Dy(0)==0);
[time3, Y3]=fplot(SOL3, [0, 150]);
plot(time3, Y3, 'b-.')
SOL4 = dsolve(D2y==F*cos(w(4)*t)-b*Dy-y*w0^2, y(0)==5, Dy(0)==0);
[time4, Y4]=fplot(SOL4, [0, 150]);
plot(time4, Y4, 'g:', 'linewidth', 2)
legend('w = 10.0', 'w=5.0', 'w=3.0', 'w=0.0')
  1 Comment
William Rose
William Rose on 31 May 2021
@Khang Ngo, I get the same graph as you, when I run your code.
Your system is equivalent to a mass-dashpot-spring system in which m=1, b=b, and k=w0^2. Since w0=10, k=100. The theoretical solution is
where yss(t)=steady state response to the driving force, and yh(t)= homogeneous solution =response to the initial conditions.
The homogeneous solution, yh(t), is the same in all four cases considered here, because the only thing that changes is the frequency of the inhomogeneous driving force. The homogeneous solution is
where and =10.000 and, from the initial conditions, C=5 and D=0. So we have
.
The steady state solution for driving force , when ω>0, is
where and . I obtained that result by standard calculus, as you can verify.
In the special case where : B=0, , and therefore simplifies to .
The theoretical solution above and the Matlab symbolic solution are different, except in the special case where w=0. I also computed the result by numerical integration, using ode45(). The numerical integration result matches the theoretical solution and does not match the symbolic solution. I conclude the symbolic solution is not working right in this problem. I don't know why.

Sign in to comment.

Answers (2)

Sulaymon Eshkabilov
Sulaymon Eshkabilov on 31 May 2021
Hi,
Everything is working ok with Symbolic MATH. You have made a mistake in your code. Here is corrected code with one numerical solution obtained with ode45().
clearvars
syms y(t)
Dy = diff(y, t);
D2y = diff(Dy, t);
w0 = 10; b = 0.1; F = 10; w = [10.0, 5.0, 3.0, 0.0]; %t = 15; ERR!!!
SOL1 = dsolve(D2y==F*cos(w(1)*t)-b*Dy-y*w0^2, y(0)==5, Dy(0)==0);
[time1, Y1]=fplot(SOL1, [0, 150]);
plot(time1, Y1, 'k-'), hold on
xlabel('t'), ylabel('y(t)'), grid on
[t, S] = ode45(@(t, u)([u(2); F*cos(w(1)*t)-b*u(2)-u(1)*w0^2]), [0, 150], [5; 0]);
plot(t, S(:,1), 'r-.'), legend('Symbolic', 'Numerical ODE45')
%%
...
Good luck.

Sulaymon Eshkabilov
Sulaymon Eshkabilov on 4 Jun 2021
for all w's:
clearvars
syms y(t)
Dy = diff(y, t);
D2y = diff(Dy, t);
w0 = 10; b = 0.1; F = 10; w = [10.0, 5.0, 3.0, 0.0]; %t = 15; ERR!!!
figure(1)
subplot(211)
for ii = 1:numel(w)
SOL1(ii,:) = dsolve(D2y==F*cos(w(ii)*t)-b*Dy-y*w0^2, y(0)==5, Dy(0)==0);
[time1, Y1]=fplot(SOL1(ii,:), [0, 150]);
plot(time1, Y1), hold all
end
xlabel('t'), ylabel('y(t)'), grid on
title('Symbolic Solution')
subplot(212)
for ii = 1:numel(w)
[t, S] = ode45(@(t, u)([u(2); F*cos(w(ii)*t)-b*u(2)-u(1)*w0^2]), [0, 150], [5; 0]);
tall(ii)={t};
SOL(ii) = {S(:,1)}; dS(ii) = {S(:,2)};
plot(t, S(:,1)), hold all
end
title('Numerical with ODE45')
xlabel('t'), ylabel('y(t)'), grid on
You may need to add legends within loop as well to show which 'w' produces what.
Good luck
  1 Comment
Khang Ngo
Khang Ngo on 4 Jun 2021
yah, thank you so much bro
I will add legends by myself, thanks

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!