Simple Undamped Forced Vibration Problem

I am trying to replicate a solution in Matlab for the following problem

x¨+k*x/m=Fo*sin*wo*t/m

using 2 first order linear differential equations in Matlab as shown below

tspan=[0 4];
y0=[.02;1]; %Initial Conditions for y(1)=x and y(2)= x dot
[t,y]=ode45(@forced1,tspan,y0); %Calls forced1.m
plot(t,y(:,2)); %y(:,1) represents the displacement and y(:,2) the velocity
grid on
xlabel('time')
ylabel('Displacement')
title('Displacement Vs Time')
hold on;

function yp = forced1(t,y)

m=20;
k=800;
f=8;
w=8;
yp = [y(2);(((f/m)*sin(w*t))-((k/m)*y(1)))];

The problem is I dont know whether Matlab considers both the complementary and particular solution. The theoretical solution is given as

x=A*sin(wn*t)+B*cos(wn*t)+(Fo*sin(wo*t)/k)/(1−(wo/wn)^2)

where the 3rd term is the particular solution assumed of the form xp=C*sin(wo*t). I am not sure how to implement this correctly in Matlab

Any ideas?

Thanks

if true
  % code
end

Answer 1

[Mischa Kim]

bugatti79, your code looks fine. In other words, you correctly implemented the differential equation and, yes, MATLAB does return the correct solution: general plus particular.

[bugatti79]

Hi Mischa

I think the bit I am misunderstanding is how to implement the code correctly to reflect the third term in some way which has wo and wn etc?

[Mischa Kim]

If you solve the DE numerically (as you do right now) MATLAB returns numeric values. You could, of course, integrate the DE twice, once with F0=0 and then subtract the two solutions from each other. But again, you will only get numeric values.

I believe what you want is to go symbolic. E.g.

 syms x(t) k m F0 w0
 DxDt   = diff(x);
 D2xDt2 = diff(x,2);
 % Define differential equation
 my_DE  = D2xDt2 + k*x/m == F0*sin(w0*t)/m;
 % Set initial conditions
 x0     = x(0)    == 0.02;
 Dx0    = DxDt(0) == 1;
 % Solve differential equation and display
 x_sol  = dsolve(my_DE, x0, Dx0);
 display(['x(t) = ',char(10),char(x_sol),char(10)])
 pretty(x_sol)

Does that help?

[bugatti79]

Hi Mischa,

Just looking at the problem now. I get an error message stating

Error in ==> AlgebraicDE at 2
 syms x(t) k m F0 w0