Question

matlab

INSTRUCTIONS Consider the spring-mass damper that can be used to model many dynamic systems Applying Newton's Second Law to a free-body diagram of the mass m yields the following ODE m습+8 +kx=F(t) (1) dt2 Where F() is a forcing function, Consider the case where the forcing function itself is a damped oscillation: Where F F(t)-Ae-Bt COS(wt) (2) For this activity we'll see how we can formulate this ODE for a solution using MATLAB that could be used to study the sensitivity of the system to changes in any of the parameter values. Task 1: Define two new variables yi and y,--. The resulting relationship between yi and y2 gives us the first equation in the system of first order ODEs that is equivalent to equation (1): dt yi dt Note this equation will ALWAYS be the first ODE in the system of first order ODEs written to be equivalent to ANY higher order ODE. Next use the definitions of y and yz to substitute for any terms involving x in equation (1) such that the only derivative term in the equation is Solve the equation for dy2. This is the dy2 dt second equation in the system. Task 2: MATLAB ODE solvers, as well as the ones we write, expect a function defining the system of ODEs to be in a specific format that outputs a column vector of evaluations of the derivatives given inputs t and a vector of y values for the current time step. Write (by hand) a MATLAB anonymous function to define this system of ODEs in this format. You can assume the parameters of the system have values in the workspace OR include them as function inputs in the appropriate position in the input list.

Answer 1

The system of ODEs is

$\dot{y}_1=y_2$

F(t) kB 7m

Write the following in a file called Dydx.m

function DX=Dydx(t,x)
global k A B m omega
DX=[x(2) ; A*exp(-B*t).*cos(omega*t)-k*x(1)/m-B*x(2)/m];
end

_________________________

The following in a separate file.

global m B k A omega
m=input('Enter the mass: ');
B=input('Damping: ');
k=input('Spring Constant: ');
x0=input('Initial position: ');
A=input('Enter force amplitude: ');
omega=input('Omega: ');
Tspan=[0 10];
[t, x]=ode45(@(t,x)Dydx(t,x),Tspan,[x0; 0]);
plot(t,x(:,1))
xlabel('Time, t','fontsize',14)
hold on
plot(t,x(:,2))
legend('Displacement','Velocity');
D=B^2-4*k*m;
if D>0
fprintf('\nSystem is overdamped.\n');
elseif D==0
fprintf('\nSystem is critically damped.\n');
else
fprintf('\nSytem is underdamped.\n');
end