Study of oscillations and how to solve in MATLAB

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')

[William Rose]

@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.

Answer 1

[Sulaymon Eshkabilov]

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.

[William Rose]

Thank you, @Sulaymon Eshkabilov!

[Sulaymon Eshkabilov]

Most Welcome!

[Khang Ngo]

Sorry, but is that all code, or just the first (w) ?

If not, can you show me more clearly

><

Thank you

Answer 2

[Sulaymon Eshkabilov]

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

[Khang Ngo]

yah, thank you so much bro

I will add legends by myself, thanks