Involves matlab Problem 2: Consider the IVP that defined the displacement of a m

Question

Involves matlab

Problem 2: Consider the IVP that defined the displacement of a mass in the mass-spring system where M = 2, c = 0.1, r(t) = cost, but the spring stiffness is not constant and determined by the equation k = 104(1.1-exp(-ly!)). Physically, this case corresponds to a spring that becomes softer at small deformation rates (small ly'l) and stiffer at large deformation rates (large ly'l) a. Write the ODE that defines displacement y(t) in this system. b. Transform this equation into a fundamental system and write this system using matrix notation. c. Develop a MATLAB computer code for numerical solving the IVP based on the d. Solve the IVP numerically with initial conditions y'(0) 0 and y'(0) 0 and e. Plot solution of the problem in the form of function y(t). RK4 method. integration step size 0.01 for the time interval from O to 10.

Explanation / Answer

F = 36; % Newtons (Applied force)
K = 12; % N/m (Spring Constant)
m = 4; % Kg (mass)4
% Case -1
B = 24.33 % N.s/m (Damping coefficient)
% To plot the responses:
% Define time vector
t = 0:0.01:10;
x = 0.3240*exp(-5.54*t) -3.3240*exp(-0.54*t)+3;
plot(t,x)
xlabel('Time (sec.)')
ylabel('Displacement (meters)')
Title(' Mass-Spring-Damper System : Case -i')

F = 36; % Newtons (Applied force)
K = 12; % N/m (Spring Constant)
m = 4; % Kg (mass)4
% Case -2
B = 13.8565 % N.s/m (Damping coefficient)
% To plot the responses:
% Define empty vectors:
T=[];
X=[];
% Define time vector
for t = 0:0.01:10;
x=-(3+5.1963*t)*exp(-1.7321*t)+3;
X=[X x];
T=[T t];
end
plot(T,X)
xlabel('Time (sec.)')
ylabel('Displacement (meters)')
Title(' Mass-Spring-Damper System : Case -i')

grid

% Mechanical Systems Example
F = 36; % Newtons (Applied force)
K = 12; % N/m (Spring Constant)
m = 4; % Kg (mass)4
% Case - 3
B = 8 % N.s/m (Damping coefficient)
% To plot the responses:
% Define empty vectors:
T=[];
X=[];
% Define time vector
for t = 0:0.01:10;
x=-3*exp(-t)*(cos(sqrt(2)*t)+1/sqrt(2)*sin(sqrt(2)*t))+3;
X=[X x];
T=[T t];
end
plot(T,X)
xlabel('Time (sec.)')
ylabel('Displacement (meters)')
Title(' Mass-Spring-Damper System : Case -i')
grid

clear all
F = 36; % Newtons (Applied force)
K = 12; % N/m (Spring Constant)
m = 4; % Kg (mass)4
% All cases, B is a vector:
% emty vectos for all cases of x(t)
% x1: case 1, X2: case -2, X3:case -3;
X1=[]; X2=[]; X3=[]; T=[];
for t = 0:0.01:10;

x1 = 0.3240*exp(-5.54*t) -3.3240*exp(-0.54*t)+3;
x2=-(3+5.1963*t)*exp(-1.7321*t)+3;
x3=-3*exp(-t)*(cos(sqrt(2)*t)+1/sqrt(2)*sin(sqrt(2)*t))+3;
X1=[X1 x1];
X2=[X2 x2];
X3=[X3 x3];
T=[T t];
end
plot(T,X1,'r',T,X2,'g',T,X3,'b')
xlabel('Time (sec.)')
ylabel('Displacement (meters)')
Title(' Mass-Spring-Damper System :All cases')
grid

x=0:0.1:10;
y1 = x.*2;
y2 = x.*3;
plot(x,y1)
hold on
plot(x,y2,'-.')
grid on
xlabel('x')
ylabel('y')
title('X Vs Y')
legend('y1','y2')
% to illustrate the figure command
x=[1,2,3,4,5,6,7,8,9,10];
y=[12 23 4 5 65 67 89];
figure(1)
plot(x)
grid
title('plot of X')
figure(2)
plot(y)
grid
title('plot of Y')
tspan=[0 4];
y0=[0.02;0];
[t,y]=ode45('unforced1',tspan,y0);
plot(t,y(:,1));
grid on
xlabel(‘time’)
ylabel(‘Displacement’)
title(‘Displacement Vs Time’)
hold on;
tspan=[0 5];
y0=[0.02;0];
[t,y]=ode45('forced',tspan,y0);
plot(t,y(:,1));
grid on
xlabel(‘time’)
ylabel(‘Displacement’)
title(‘Displacement Vs Time’)
hold on;
tspan = [0 5];
y0 = [1.57;0];
[t,y] = ode45('linear',tspan,y0)
plot(t,y(:,1))
grid on;
xlabel(‘Time’)
ylabel(‘Theta’)
title(‘Theta Vs Time’)
hold on;

tspan=[0 2]
y0=[0;0]
[t,y]=ode45('coulomb',tspan,y0)
plot(t,y(:,1))
grid
xlabel(‘Time’)
ylabel(‘Theta’)
title(‘Theta Vs Time’)
hold on;
function yp = airdrag(t,y)
m = 10;
cd = 0.2;
g = 9.81;
w = m*g;
yp = zeros(4,1);
yp(1) = y(2);
yp(2) = ((-cd/m)*y(2)*(y(2)^2+y(4)^2)^(0.5));
yp(3) = y(4);
yp(4) = ((-w/m)-(cd/m)*y(4)*(y(2)^2+y(4)^2)^(0.5));
tspan=[0 5];
y0=[0;100;0;10]
[t,y]=ode45('airdrag',tspan,y0);
plot(y(:,1),y(:,3))
grid
xlabel(‘X-Displacement’)
ylabel(‘Y-Displacement’)
title(‘X vs Y Displacement’)
hold on;

tspan=[0 10];
y0=[1.57;0]
[t,y]=ode45('pendulum1',tspan,y0);
plot(t,y(:,1));grid
xlabel(‘Time’)
ylabel(‘Theta’)
title(‘Theta Vs Time’)
hold on;

function yp = pendulum2(t,y)
m=10;
g=981;
l=5;
Ca=14;
yp=[y(2);((-g/l)*(y(1))-(Ca/m)*(y(2)))];

tspan=[0 10];
y0=[1.57;0]
[t,y]=ode45('pendulum2',tspan,y0);
plot(t,y(:,1));grid
xlabel(‘Time’)
ylabel(‘Theta’)
title(‘Theta Vs Time’)
hold on;

function yp= pend(t,y)
M1=5;
M2=5;
g=9.81;
l1=1;
l2=1;
w2=M2*9.81;
w1=M1*9.81;
yp=zeros(4,1);
52
yp(1)=y(2);
yp(2)=-(w1+w2)*sin(y(1))+M2*l2*(y(4)^2)*sin(y(3)-y(1));
yp(3)=y(4);
yp(4)=-w2*sin(y(3))-M2*l1*(y(2)^2)*sin(y(3)-y(1));
Similarly, to store the coefficient matrix a separate function file is written which is stored as
‘MM.m’.
% the following function contains the mass matrix.
%it is separately stored in a file named, mass.m
function m = mass(t,y)
M1=5;
M2=5;
g=9.81;
l1=1;
l2=1;
m1=[1 0 0 0];
m2=[0 (M1+M2)*l1 0 M2*l2*cos(y(3)-y(1))];
m3=[0 0 1 0];
m4=[0 M2*l1*cos(y(3)-y(1)) 0 M2*l2];
m=[m1;m2;m3;m4];
orresponding functions depending on the value

using ode113.
%it then plots the first variable.
tspan=[0 10]
y0=[0.5233;0;0.5233;0]
options=odeset('mass','M(t,y)')
[t,y]=ode113('indmot_ode',tspan,y0,options)
subplot(2,1,1)
plot(t,y(:,1))
grid
xlabel('Time')
ylabel('Theta1')
53
subplot(2,1,2)
plot(t,y(:,3))
grid
xlabel('Time')
ylabel('Theta2')

for i = 1:50
omega(i)=2*pi*f(i); %omega in terms of frequency
omega2(i)=omega(i)*omega(i); % squaring omega
a11=-omega2(i)*m+k; % representing the left hand…
a12=omega(i)*c; % matrix of the single matrix…
a21=-omega(i)*c; % equation
a22=-omega2(i)*m+k;
a=[a11 a12;a21 a22];
b=inv(a);
c1=[0;0;fi];
d(1,i)=b(1,:)*c1;
d(2,i)=b(2,:)*c1;
d(3,i)=b(3,:)*c1;
d(4,i)=b(4,:)*c1;
x(1,i)=sqrt(abs(d(1,i))^2+abs(d(3,i))^2);
x(2,i)=sqrt(abs(d(2,i))^2+abs(d(4,i))^2);
p(1,i)=atan(d(1,i)/d(3,i))*180/pi;
if p(1,i)<0 % to check whether the angle is
negative or not.
p(1,i)=180+p(1,i);
else

p(1,i)=p(1,i);
end

p(2,i)=atan(d(2,i)/d(4,i))*180/pi;

if p(2,i)<0
if d(4,i)<0
p(2,i) = -180 + p(2,i)
else
p(2,i)=180+p(2,i);
end

else
p(2,i)=p(2,i);
end
end

Related Questions

Navigate

• Browse (All)
• Browse I
• Subjects

---

---

Integrity-first tutoring: explanations and feedback only — we do not complete graded work. Learn more.