Autumn Schuelka
Wednesday 8:35am
Palou-Genevive Toutain
MAT 275 MATLAB Lab 5
1(a)
function LAB05ex1
m = 1; % mass [kg]
k = 4; % spring constant [N/m]
omega0=sqrt(k/m);
y0=0.1; v0=0; % initial conditions
[t,Y]=ode45(@f,[0,10],[y0,v0],[],omega0); % solve for 0<t<10
y=Y(:,1); v=Y(:,2); % retrieve y, v from Y
figure(1); plot(t,y,'b+-',t,v,'ro-'); % time series for y and v
grid on;
%------------------------------------------------------
function dYdt= f(t,Y,omega0)
y = Y(1); v= Y(2);
dYdt = [v; -omega0^2*y];
end
end
The blue curve represents the displacement that the weight travels
attached to the spring. It only moves from 0.1 to -0.1 meters and it starts
at the position of 0.1m. The red curve represents the velocity of the mass,
this is obvious because it begins at zero and then correlates with the
changes of position. When the position is at its equilibrium y=0, then the
velocity is at a maximum or minimum value. When the position is at its
maximum or minimum, the velocity is zero. This means that the blue curve is
y=y(t).
1(b)
If we look at the graph, the distance between two peaks is a period. One
period is at about 3.25 and the next peak is at about 6.25. The distance
between the two peaks is approximately 3.
Analytically we can calculate ω0 by substituting the values of k and m in.
𝜔0=√𝑘
𝑚=√4
1=2
Then we can determine the period
𝑇= 2𝜋
𝜔0=2𝜋
2=𝜋≈3.14
1(c)
In this example the mass will never come to rest. According to the graph, v
and y never come to zero at the same time so this system doesn’t stop. Also,
assuming ideal conditions, since there is no friction in the system, it will
never stop.
1(d)
The amplitude is the distance from the equilibrium position to the peak
position. According to the graph, the equilibrium position is 0 and the max
peak is at 0.1, this means the amplitude of y is 0.1.
1(e)
The maximum velocity corresponds with the position where the displacement is
at equilibrium. To find the values of t where the velocity is maximum it is
easier to zoom into where the displacement is zero. The maximum velocity in
magnitude is 0.2m/s and this occurs at t=
4,3
4,5
4,7
4,9
4 𝑎𝑛𝑑 11
4 . These exact
values can be determined because they occur distance apart starting at 𝜋
4.
1(f)
Since 𝜔0=√𝑘
𝑚 as k increases, 𝜔0 increases, and as k decreases 𝜔0 decreases.
But as m increases, 𝜔0 decreases, and as m decreases, 𝜔0 increases. Compared
to the period, 𝑇= 2𝜋
𝜔0 as 𝜔0 increases, T decreases, and as 𝜔0 decreases, T
increases.
Overall, if k is larger (the spring is more stiff), the period is shorter,
but with a larger m (larger mass), the period is longer.
With m=5 and k=4, the period is much longer (approx. 7) than the original.
function LAB05ex1
m = 5; % mass [kg]
k = 4; % spring constant [N/m]
omega0=sqrt(k/m);
y0=0.1; v0=0; % initial conditions
[t,Y]=ode45(@f,[0,10],[y0,v0],[],omega0); % solve for 0<t<10
y=Y(:,1); v=Y(:,2); % retrieve y, v from Y
figure(1); plot(t,y,'b+-',t,v,'ro-');
grid on;
%------------------------------------
function dYdt= f(t,Y,omega0)
y = Y(1); v= Y(2);
dYdt = [v; -omega0^2*y];
end
end
With m=1 and k=16, the period is much smaller (approx. 2.5) than the
original.
function LAB05ex1
m = 1; % mass [kg]
k = 16; % spring constant [N/m]
omega0=sqrt(k/m);
y0=0.1; v0=0; % initial conditions
[t,Y]=ode45(@f,[0,10],[y0,v0],[],omega0); % solve for 0<t<10
y=Y(:,1); v=Y(:,2); % retrieve y, v from Y
figure(1); plot(t,y,'b+-',t,v,'ro-'); % time series for y and v
grid on;
%------------------------------------------------------
function dYdt= f(t,Y,omega0)
y = Y(1); v= Y(2);
dYdt = [v; -omega0^2*y];
end
end
2(a)
Modified function LAB05ex1 to graph E
function LAB05ex1
m = 1; % mass [kg]
k = 4; % spring constant [N/m]
omega0=sqrt(k/m);
y0=0.1; v0=0; % initial conditions
[t,Y]=ode45(@f,[0,10],[y0,v0],[],omega0); % solve for 0<t<10
y=Y(:,1); v=Y(:,2); % retrieve y, v from Y
figure(1); plot(t,y,'b+-',t,v,'ro-'); % time series for y and v
grid on;
E=m*v.^2/2+k*y.^2/2; % this defines the equation
figure(2); plot(t,E); % this opens a new plot
ylim([0,0.04]) % this changes the y axis limits
%------------------------------------------------------
function dYdt= f(t,Y,omega0)
y = Y(1); v= Y(2);
dYdt = [v; -omega0^2*y];
end
end
We expect the graph of E to be constant because energy should be conserved
since there is no damping or external forces on the system. In reality it is
not but it is very close. The graph before using ylim was very jumpy but
when observing the y axis scale, it is obvious that the change was very
small, after using the ylim command a more constant looking graph can be
shown.
Graph of E before using ylim Graph of E after using ylim
2(b)
By showing that 𝑑𝐸
𝑑𝑡 =0, this means that E is constant. From the definition of
kinetic energy we can show that 𝑣′=𝑦′′ = −𝜔02𝑦= −𝑘
𝑚𝑦. Then we can substitute
values into the original equation of 𝐸=1
2𝑚𝑣2+1
2𝑘𝑦2 after we differentiate
with respect to v and y 𝐸=1
2𝑚2𝑣𝑣′+1
2𝑘2𝑦𝑦′
simplify =𝑚𝑣𝑣′+𝑘𝑦𝑦′
Substitute v for y’ and v’ for −𝑘
𝑚𝑦
=𝑚𝑦′(−𝑘
𝑚𝑦)+𝑘𝑦𝑦′
M cancels to give −𝑘𝑦𝑦′+𝑘𝑦𝑦′=0
This is expected because if the derivative of E is zero, then the slope of
its graph is zero meaning that E is constant.
2(c)
The phase plot is an ellipse and is closed which indicates a periodic system
as expected. The maximum and minimums of the v and y correspond with the
sinusoidal graph of v and y. This graph never comes close to the origin
because this would mean that both v and y are zero, meaning that the
velocity and position are both zero which would result in a stopped spring.
This would never happen because this system is in ideal conditions and there
is no friction to slow its motion.
function LAB05ex1
m = 1; % mass [kg]
k = 4; % spring constant [N/m]
omega0=sqrt(k/m);
y0=0.1; v0=0; % initial conditions
[t,Y]=ode45(@f,[0,10],[y0,v0],[],omega0); % solve for 0<t<10
y=Y(:,1); v=Y(:,2); % retrieve y, v from Y
figure(1); plot(t,y,'b+-',t,v,'ro-'); % time series for y and v
grid on;
E=m*v.^2/2+k*y.^2/2;
figure(2); plot(t,E);
ylim([0,0.04])
figure(3) % this opens a new plot
plot(y,v); % this plots v vs y
xlabel('y'); ylabel('v'); grid on; % this labels the axes and adds a grid
%------------------------------------------------------
function dYdt= f(t,Y,omega0)
y = Y(1); v= Y(2);
dYdt = [v; -omega0^2*y];
end
end
Phase Plot v vs y
3(a)
The values of t where |y(t)|<0.01 are greater than approximately 3.8. To be
exact we can plot a horizontal line at y=0.1 and y=-0.1 so that we can be
sure and the graph reveals that when t>3.81 then | y(t)|<0.01.
3(b)
The maximum velocity in an undamped system corresponds with the position
where the displacement is at equilibrium. Since this system is damped, the
maximum magnitude of velocity is where the first peak in the red graph
occurs.
After zooming in on the first
peak of the velocity curve
which is actually a local
minimum, we can see that the
value is approximately
0.14197m/s in magnitude and
this occurs at t=0.714775
This is a smaller velocity than
the undamped system because of
the damping effect.
3(c)
We expect that as the value c increases, the motion of the mass will
decrease because if the effects of damping are increased (represented by c)
then the number of oscillations will decrease, resulting in a slowing of the
system at a faster rate.
When c=2 the system is classified at underdamped because there is greater
than one oscillation occurring. In this example there are about two
oscillations.
When c=4 the system could be either critically damped, or overdamped. It is
not easy to tell yet because the critically damped system is classified as
the exact moment where the damping on the system restricts to no
oscillations. It looks like the system begins to go straight to the
equilibrium position and then just stays there.
When c=6 the system is classified as overdamped because when c=4, the system
already showed that there were no oscillations so increasing the damping
will only result in an overdamped system. The graph also shows that it takes
the system longer to reach the equilibrium position than when c=4.
When c=8 the system is classified as overdamped. It can be observed that it
takes the system even longer to reach the equilibrium position then when
c=6.
3(d)
The critical value of c can be found by solving the homogeneous equation to
the differential equation 𝑦′′+2𝑝𝑦′+𝜔0
2𝑦=0
𝑟2+2𝑝𝑟+𝜔0
2=0
𝑟= −2𝑝±√4𝑝2−4𝜔0
2
2
𝑟=−𝑝±√𝑝2−𝜔0
2
The repeated root occurs only if 𝑝2−𝜔0
2=0 so that both solutions of r will
be -p
Since 𝑝= 𝑐
2𝑚 then this can be substituted into the variables under the root
to continue to solve for the desired value of c
𝑐
2𝑚2−𝜔0
2=0
Then solve for c 𝑐=√4𝑚𝜔0
From here we can substitute value for m and 𝜔0 from the initial conditions
𝑐=√4∗1∗2
𝑐=4
So the graph that shows c=4 is a graph showing the critically damped system
because there are no oscillations and this is confirmed where it is the
smallest value of c where this occurs.
4(a)
When the Energy of the original system is plotted, at first it can be
observed that the energy of the system is not conserved because the graph
seems to be steadily decreasing.
function LAB05ex1a
m = 1; % mass [kg]
k = 4; % spring constant [N/m]
c = 1; % friction coefficient [Ns/m]
omega0 = sqrt(k/m); p = c/(2*m);
y0 = 0.1; v0 = 0; % initial conditions
[t,Y]=ode45(@f,[0,10],[y0,v0],[],omega0,p); % solve for 0<t<10
y=Y(:,1); v=Y(:,2); % retrieve y, v from Y
figure(1); plot(t,y,'b+-',t,v,'ro-'); % time series for y and v
E=m*v.^2/2+k*y.^2/2; % this is the energy equation
figure(2); plot(t,E); % this plots energy
ylim([-.01,0.03]) % this alters the y axis
grid on
%------------------------------------------------------
function dYdt= f(t,Y,omega0,p)
y = Y(1); v= Y(2);
dYdt = [v;-omega0^2*y-2*p*v]; % fill-in dv/dt
When the ylim command is added, then a clearer representation is shown that
the energy of the system decreases to zero.
Original graph Graph using ylim command
4(b)
To show 𝑑𝐸
𝑑𝑡 <0 when c>0 we need to use a similar method as before in question
2. From the definition of kinetic energy we can show that
𝑣′=𝑦′′ = −𝜔02𝑦−2𝑝𝑣= −𝑘
𝑚𝑦−𝑐
𝑚𝑣. Then we can substitute values into the
original equation of 𝐸=1
2𝑚𝑣2+1
2𝑘𝑦2 after we differentiate with respect to v
and y 𝐸=1
2𝑚2𝑣𝑣′+1
2𝑘2𝑦𝑦′
simplify =𝑚𝑣𝑣′+𝑘𝑦𝑦′
Substitute v for y’ and v’ for −𝑘
𝑚𝑦−𝑐
𝑚𝑣
=𝑚𝑦′(−𝑘
𝑚𝑦−𝑐
𝑚𝑦′)+𝑘𝑦𝑦′
Simplify =𝑚𝑦′(−𝑘𝑦−𝑐𝑦′
𝑚)+𝑘𝑦𝑦′
=𝑦′(−𝑘𝑦−𝑐𝑦′)+𝑘𝑦𝑦′
=(−𝑘𝑦𝑦′−𝑐2𝑦′)+𝑘𝑦𝑦′
𝑑𝐸
𝑑𝑡 =−𝑐2𝑦′
When c>0 then the term -c2y’ is negative, showing that 𝑑𝐸
𝑑𝑡 <0.
Similarly, when c<0, then the term -c2y’ is positive showing that 𝑑𝐸
𝑑𝑡 >0.
4(c)
function LAB05ex1a
m = 1; % mass [kg]
k = 4; % spring constant [N/m]
c = 1; % friction coefficient [Ns/m]
omega0 = sqrt(k/m); p = c/(2*m);
y0 = 0.1; v0 = 0; % initial conditions
[t,Y]=ode45(@f,[0,10],[y0,v0],[],omega0,p); % solve for 0<t<10
y=Y(:,1); v=Y(:,2); % retrieve y, v from Y
figure(1); plot(t,y,'b+-',t,v,'ro-'); % time series for y and v
E=m*v.^2/2+k*y.^2/2;
figure(2); plot(t,E);
ylim([-.01,0.03])
grid on
figure(3) % this opens a new plot
plot(y,v); % this plots v vs y
xlabel('y'); ylabel('v'); grid on; % this labels the axes and adds a grid
%------------------------------------------------------
function dYdt= f(t,Y,omega0,p)
y = Y(1); v= Y(2);
dYdt = [v;-omega0^2*y-2*p*v]; % fill-in dv/dt
Phase plot v vs y
This phase plot shows a circular graph, but it is not enclosed so the
behavior cannot be classified as periodic. The values of x and y never
return to the original position, or any position twice. Since the system is
damped, both v and y decrease. They graph does reach the origin because
there is a point where both v and y are zero, this occurs at the end.
