E12 Assignment 04

Description: 

 

Problems: 
1. Some Laplace Transform Functions and Properties

For the following problems you may use the properties of the Laplace Transform (link to table), and the Laplace transform of an exponential.

No integration should be necessary.

a)  Derive the Laplace Transform of a sinusoid  (you needn't do any integrals to solve this one).  Recall Euler's identity for sine.

 

b)  Derive the Laplace Transform of a decaying sinusoid  (you needn't do any integrals to solve this one). 

 

c) Given the Laplace Transform of the sine function

use the differentiation property of the Laplace transform to find the Laplace transform of the cosine.   In other words show that wwithout taking the derivative in the time domain.

Solution: 

a)

b)



... or use complex shift property

c)

2. Derive Convolution Property of Laplace Transform

Derive the convolution property of the Laplace Transform (assume f1(t) and f2(t) are causal (i.e., they are zero for t<0)).  The convolution property states that the Laplace Transform of the convolution of two function is simply the product of the Laplace Transforms of the functions - this makes convolution easy.

Hint: Write the definition of the convolution (in terms of integrals) and then use it in the definition of the Laplace Transform.  This will give you a double integral - you can switch order of integration and with some judicious change of variables (and pulling constants out of integrals as appropriate) the answer develops.

Solution: 

See here for explanation of the steps in the solution

3. E12 Spring Mass Dashpot - Laplace

 

Consider the system shown.

a) Show that the differential equation relating postion, x(t), to applied force, fa(t) is

b) Write an expression for X(s) if fa(t)=δ(t).

c) Determine an expression for the unit impulse response (all initial conditions are zero at t=0-) if m=1 kg, k=25 N/m, b=11.125 N-s/m.

d) Find the unit impulse response if m=1 kg, k=25 N/m, b=10 N-s/m.

e) Find the unit impulse response if m=1 kg, k=25 N/m, b=6 N-s/m.

f) Plot all results on a single graph with a legend and axis labels for t=0 to 3 seconds.

Solution: 

% Plot the impulse response.
t=linspace(-1,3,1000);    % Create time vector from -0.1 to 2 in 1000 steps.

y_over = (exp(-3.125*t)-exp(-8*t)).*heaviside(t)/4.875;
y_crit = t.*exp(-5*t).*heaviside(t);
y_under = exp(-3*t).*sin(4*t).*heaviside(t)/4;

plot(t,y_over,t,y_crit,'g-.',t,y_under,'r:');

title('Impulse Response, Position of Mass, E12 Assignment'); % Plot title
xlabel('Time (s)');                  %   x-axis label
ylabel('Position (m)');              %   y-axis label
legend('Over','Crit','Under');

4. E12 Spring Mass Dashpot - RK2

 

Consider the system shown. All initial conditions are zero at t=0-, and the input is a unit impulse, fa(t)=δ(t).

a) Write an expression for X(s) and use it to find the initial condition on position at t=(0+), i.e. x(0+).

b) Write an expression (in the Laplace Domain) for the velocity of the object, V(s), and use it to find the initial condition on velocity at t=(0+), i.e. v(0+).

c) Write a second order Runge-Kutta program to solve for x(t) for t=0 to 3 seconds if m=1 kg, k=25 N/m, b=6 N-s/m. Set step size so h=0.05 seconds.

d) Plot the simulation results along with the exact solution. Make sure graph has legend and axis labels.

e) What would you have to change in your simulation if the value of the friction coefficient, "b," changed?

Solution: 

a)

b) Velocity is derivative of postion

c,d)

% Plot the impulse response.
h=0.05;
t=0:h:3;    % Create time vector from -0.1 to 2 in 1000 steps.

% y_over = (exp(-3.125*t)-exp(-8*t)).*heaviside(t)/4.875;
% y_crit = t.*exp(-5*t).*heaviside(t);
y_under = exp(-3*t).*sin(4*t).*heaviside(t)/4;

m=1; b=6; k=25;

q0 = [0; 1];         % Initial Condition (vector)
A = [0 1; -k/m -b/m];  % A Matrix

qstar = zeros(2,length(t));  % Preallocate array (good coding practice)

qstar(:,1) = q0;             % Initial condition gives solution at t=0.
for i=1:(length(t)-1)
  k1 = A*qstar(:,i);                % Approx for y gives approx for deriv
  q1 = qstar(:,i)+k1*h/2;           % Intermediate value
  k2 = A*q1;                        % Approx deriv at intermediate value.
  qstar(:,i+1) = qstar(:,i) + k2*h; % Approx solution at next value of q
end
plot(t,y_under,'x',t,qstar(1,:),'+:');        % ystar = first row of qstar
xlabel('Time (s)');                  %   x-axis label
ylabel('Position (m)');              %   y-axis label
legend('Exact','Approximate');

e) The only thing that needs to change is the value of "b" in the code for the simulation. That is part of the beauty of simulation - it doesn't care if the problem is over, under, or critically damped (or even the order of the problem - for an 8th order problem only the A matrix and initial condition would need to be changed). The exact response would change as well, but the question only asked about the simulation.

5. E12 Spring Mass Dashpot - Simulink

 

Consider the system shown. All initial conditions are zero at t=0-, and the input is a unit impulse, fa(t)=δ(t).

a) Make a simulink model to solve for x(t) for t=0 to 3 seconds if m=1 kg, k=25 N/m, b=6 N-s/m. Include the diagram with your homework.

b) Plot the simulation results and the exact solution on two separate graphs and quantitatively compare to the exact result (or if you are ambitious, send the results to MATLAB and plot the simulation results along with the exact result).

c) What would you have to change in your simulation if the value of the friction coefficient, "b," changed?

Solution: 

a) diagram is below. Initial condition on velocity integrator (left) is v(0+)=1/m=1, initial condition on position is x(0+)=0. We form the second derivative of x as the sum of the first derivative, and x itself.

Input to leftmost integrator is acceleration. Output of left integrator is velocity. Output of right integrator is position.

b)

>> t=x.time;
>> plot(t,x.data,'x-',t, exp(-3*t).*sin(4*t).*heaviside(t)/4,'+:')
>> title('Impulse Response'); xlabel('Time (s)'); ylabel('x(t) (m)');
>> legend('Simulink','Exact');

c) The only thing that needs to change is the value of "b" in the code for the simulation. That is part of the beauty of simulation - it doesn't care if the problem is over, under, or critically damped. The exact response would change as well, but the question only asked about the simulation.

6. Partial Fractions, Laplace

Find the Inverse Laplace Transform of the following function using partial fraction expansion.

a)

b)

c)

d) Check your results with Matlab. 

...Turn in the Matlab code and output and demonstrate that it agrees with your hand calculations.
Solution: 

a)

 

b)

c)

Since we have a repeated root in the denominator we cross multiply

Equate coefficients of like powers of s, and solve

 

d)

>> %Part a
>> [r,p,k]=residue([0 2 5],[1 5 6])
r =
    1.0000
    1.0000
p =
   -3.0000
   -2.0000
k =
     []

>> %Part b
>> [r,p,k]=residue([3 0 12],[1 7 6 0])
r =
     4
    -3
     2
p =
    -6
    -1
     0
k =
     []

>> %Part c
>> [r,p,k]=residue([1 0 1],[1 2 1 0])
r =
     0
    -2
     1
p =
    -1
    -1
     0
k =
     []

 

7. Pulse response of 2nd order system

A system is defined by the differential equation

a) The graph below shows f(t).  Express f(t) as the difference of two step functions (multiplied by a constant).

b) Find F(s)

c) Find an expression for y(t).  Assume all initial conditions at t=0- are zero.

d) Plot f(t) and y(t) for t=-1 to t=10 seconds with Matlab.  Include a legend on your plot.  Include your code and your graph.

Solution: 

a)

b)

c)

As an intermediate step, let's solve for a function of the form

We can find the coefficients from cross multiplication:

 

The last step was done with the "generic decaying oscillatory" in the Laplace Transform table.  Looking back at our expression for y(t), we find that it is simply x(t) minus a delayed version of x(t) (In other words Y(s) is just X(s) minus X(s) delayed (multiplied by e-1.5s)).

d)

t=linspace(-1,10,1000);     %Time vector
%Define output (continue over two lines).
y=2*(1-exp(-2*t).*(cos(5*t)+2/5*sin(5*t))).*(t>=0)-...
    2*(1-exp(-2*(t-1.5)).*(cos(5*(t-1.5))+2/5*sin(5*(t-1.5)))).*((t-1.5)>=0);
% Plot input and output
plot(t,2*((t>=0)-((t-1.5)>=0)),'Linewidth',2);
hold on

% Label Graphs
plot(t,y,'r','Linewidth',2);
xlabel('Time (S)');
legend('Input','Output');
axis([-1 10 -1 3])
hold off

8. Complicated Partial Fraction, Laplace (by hand)

Find the Inverse Laplace Transform of the following function by hand using partial fraction expansion (without MATLAB's "residue" command); you may use MATLAB or other tool to solve simultaneous equations involved in Inverse Transform.   Laplace Transform table.

Solution: 

Cross multiply

Equate like powers of s

This looks like a pretty full 5x5 set of equations – so either use a computer, or find A and B by cover up method and solve for others

>> X=[1 0 1 1 0; 6 1 4 4 1; 22 2 14 4 4; 48 10 20 0 4; 40 0 0 0 0];
>> inv(X)*[5 25 76 132 40]'
ans =
    1.0000
    2.0000
    3.0000
    1.0000
    1.0000