# E12 Assignment 10

Description:

Feel free to use MATLAB for any problem larger than 2x2 and to check all solutions with MATLAB (or your calculator, abacus...).

Problems:
1. Oscillation of Four Story Building (four masses with springs)

The system shown is a very crude model of a building with no damping. The masses can only move in a horizontal direction.

a) draw free body diagrams

b) write equations of the form

c) If m1=m2=m3=m4=1, and k=1, find mode shapes and associated frequencies (you can use the code from class).  Note: these mode shapes are similar to those of the air pressure in wind instruments.

Solution:

a)

b)

c) Pertinent equation is:

`A=[-2 1 0 0; 1 -2 1 0; 0 1 -2 1; 0 0 1 -1];`

```----------------
Eigenvectors and values
v =
-0.4285   -0.6565    0.5774    0.2280
0.6565    0.2280    0.5774    0.4285
-0.5774    0.5774   -0.0000    0.5774
0.2280   -0.4285   -0.5774    0.6565
d =

-3.5321         0         0         0
0   -2.3473         0         0
0         0   -1.0000         0
0         0         0   -0.1206

----------------
Frequencies
omega =
1.8794
1.5321
1.0000
0.3473```
2. Four Story Building, excitation of lowest frequency mode

For the previous problem

a) find initial conditions such that only the lowest frequency mode is excited.

b) simulate the system and plot the trajectories of the masses.  (you can use the code from class)

Solution:

a) Choose as initial conditions a multiple of the eigenvector associated with the lowest frequency (in this case the 4th eigenvector).  We can choose this eigenvector in Matlab

`x0=v(:,4)           %Initial condition (size must match that of "A" array).`

b)

Note: we can get animation to work properly with the following lines (in place of the original lines in code).

```% Animate (animates objects as beads on a string.
.
.
.
yRope=[0 yBeads];   % The y deflection of the rope (0 at left end)
.
.
.
set(ropeLine,'Ydata',[0 x(:,i)']);    % Change y values of rope
```

3. Oscillations in LC Circuit (Double Pendulum Analog)

For the circuit shown L1=1,  L2=2, C=1.

a) Determine the circuit equation in the form:

b) Find frequencies of vibrations, and mode shapes (do by hand, you can check with Matlab).

c) Draw the force-current analog of the system.

Solution:

a)

b)

Solve for ω2.

Find v1.

Find v2.

Check with Matlab

```>> [v,d]=eig([-3/2 1/2; 1/2 -3/2])
v =
0.7071    0.7071
-0.7071    0.7071
d =
-2     0
0    -1```

c)

4. Oscillations, 3 Masses with Springs

For the system shown  the masses can only move in a horizontal direction.

a) Draw free body diagrams

b) Write equations of the form

c) If m1=m2=m3=1, and k=1, find mode shapes and associated frequencies (you can use the code from class).  To think about: are these approximately what you expected?

d) Find the weights (i.e., the "γ" terms) associated with each eigenvector and plot trajectories for the initial conditions .

Solution:

a)

b)

c) Pertinent line is:

`A=[-2 1 0; 1 -2 1; 0 1 -2];`

```----------------
Eigenvectors and values

v =
0.5000   -0.7071   -0.5000
-0.7071    0.0000   -0.7071
0.5000    0.7071   -0.5000

d =
-3.4142         0         0
0   -2.0000         0
0         0   -0.5858

----------------
Frequencies

omega =
1.8478
1.4142
0.7654```

Lowest frequency mode is smoothest.  Increasing frequency with increasing number of extrema.

d) only modes one and three (with mass 1 and mass 3 in same direction) are excited.

```----------------
Mode amplitudes

gam =
1.0000
0.0000
-1.0000```

5. Oscillations, 3 Masses with Springs, Middle mass is heavier

For the system shown

As you are doing this problem, think about the how the system would behave as the middle mass gets to be much heavier than the other two.   This should allow you to figure out in advance approximately what the mode shapes and frequencies should be.

a)   If m1=m3=1, m2=10, and k=1, find mode shapes and associated frequencies (you can use the code from class).   Try to figure out approximately what the mode shapes and frequencies will be before solving.  Recall that for a mass-spring system the natural frequency is (k/m)½.

b) Find the weights  (i.e., the "γ" terms) associated with each eigenvector and plot trajectories for the initial conditions .

Solution:

a) Because the middle mass is so heavy, it won't generally move much.

There will be two modes with similar frequencies in which the outer masses are in and out of phase with each other, but m2 barely moving.  These frequencies will be close.  Theses are modes 2 and 3 in the solution below.   Since the middle mass is not moving it is as if m1 and m3 have an effective spring constant of 2*k attached (this is because there are two springs attached, and m3 is practically fixed).  The frequency will be  (2*k/m1)½=1.41.

The other mode has m2 moving, and m1 and m3 just going along without affecting the dynamics of the problem very much.  Therefore we expect m1 and m3 to be in phase with m2, but with an amplitude about half as large.  If we ignore m1 and m3, then the mass m2 has 2 springs attached each of magnitude keff=k/2 (this is the equivalent spring constant of two springs in series between the mass and the wall).  So the total spring constant is approximately k, and the frequency should be (2*keff/m2)½=(1/10)½=0.32

Pertinent line is:

`A=[-2 1 0; 0.1 -0.2 0.1; 0 1 -2];`

```----------------
Eigenvectors and values
v =
-0.4215    0.7071    0.7052
-0.8029    0.0000   -0.0740
-0.4215   -0.7071    0.7052
d =
-0.0950         0         0
0   -2.0000         0
0         0   -2.1050

----------------
Frequencies
omega =
0.3082
1.4142
1.4509```

b)

```----------------
Initial Conditions
x0 =
1
0
1

----------------
Mode amplitudes
gam =
-0.1239
-0.0000
1.3440```

6. State Space as Eigen Problem

This problem demonstrates that zero-input state space problems can be framed as an eigenvector problem.

Consider the state space system (with no input) given by:

Assume the solution is of the form q(t)=cvest where c and s are scalars, and v is a 2x1 vector.

a) Put q(t) (and its derivative) into the differential equation, and show that an eigenvalue problem results (A-sI)v=0 (or Av=sv), where I is the identity matrix.

b) Find the eigenvalues, s1 and s2.

c) Find the eigenvectors, v1 and v2.

Just as in E11, we form the generic solution as the weighted sum of these solutions: where the constants c1 and c2 are determined from initial conditions.  However, in this case we easily represent the initial conditions in matrix form.

or q(0)=Vc  where q(0) is the (vector) initial condition, V is a square matrix in which each column is an eigenvector, and c is a column vector whose elements are the unknown coefficients (or weights) of the individual solutions.

d) If , find c1 and c2.

e)  If  , solve for y(t) and plot it for 7 seconds.

f) Choose a set of initial conditions such that only the faster exponential is present in the solution.

Solution:

a)

b)

c)

With MATLAB

```>> A=[-5 -2; 2 0]
A =
-5    -2
2     0
>> [v,d]=eig(A)
v =
-0.8944    0.4472
0.4472   -0.8944
d =
-4     0
0    -1```

Eigenvalues are d(1,1) and d(2,1)  or -4, -1.   This agrees with our results.
Eigvenvectors are v(:,1) and v(:,2) or [-0.9;0.45], [0.45,-0.9].  This agrees with our results (within a multiplicative constant).

d)

```>> c=inv(v)*[1; -1]
c =
-0.7454
0.7454```

Note you may get different values for this part of the problem depending on your eigenvectors.

e)

```>> % Purely numerical solution
>> t=linspace(0,7);
>> s=diag(d)
s =
-4
-1
>> % Solve for state variables
>> q = c(1)*v(:,1)*exp(s(1)*t) + c(2)*v(:,2)*exp(s(2)*t);
>> % Find y and plot it.
>> y = [1 3/4]*q;
>> plot(t,y)```

```>> % or symbolic solution|
>> syms t
>> [v,d]=eig([-5 -2;2 0]);
>> s=diag(d);
>> q0=[1; -1];
>> c=inv(v)*q0;
>> q = c(1)*v(:,1)*exp(s(1)*t) + c(2)*v(:,2)*exp(s(2)*t)
>> y=[1 3/4]*q
y =
(5*exp(-4*t))/12 - exp(-t)/6
>> pretty(y)

5 exp(-4 t)   exp(-t)
----------- - -------
12           6```

f) The first eigenvalue represents the faster decay, so any initial conditions that is a multiple of the first eigenvector will not include the other in the solution.

```>> % Check with MATLAB
>> q0 = [20; -10]  % This is a multiple of the eigenvector
q0 =
20
-10
>> c=inv(v)*q0
c =
-22.3607
0
>> % Note that c(2) (contribution from slower exponential) is zero.```