11 views (last 30 days)

Show older comments

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

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

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

Start Hunting!