Lab 5 -- A Swinging Skyscraper
Lab 5 (click here for pdf) concerns the modeling of the swaying of a skyscraper. The theoretical justification for the model of skyscraper movement is given in the attached handout at the end of the lab. Be sure to read through the lab handout, and then get started.
Background:
This lab's premise: You have a set of actual data showing the position of a skyscraper at various points in time.
Objectives:
1. Determine which of the following two equations best models skyscraper movement:
y'' + p*y' + q*y = 0; [ Equation A ]
y'' + p*y' + q*y = k*y^3 [ Equation B ]
2. Use the given models of skyscraper movement to come up with a system of ODEs that accurately fits the actual data. These means you'll adjust the constants (p, q, and/or k) in your given equations until you arrive at a solution that actual closely fits the actual data.
3. Discern how well the model fits the data, and use your models to predict whether or not there is a danger that the skyscraper will collapse.
Getting started
Let's first consider objective (1). You have Equation A and Equation B (see above), and it wants to know which model is superior. Hint: Consider that one equation is just a special case of the other. Also look at the attached handout's information regarding the P-delta effect (due to gravity) and take that into consideration when making your response.
Now let's consider objective (2). In this step you will be trying to find, by trial and error, values of p, q, and k that make these differential equations fit the actual data closely. The biggest chore in this lab is to get started doing that. Here are some steps that will help you get started:
(1) "Take k = 0 and find values of p and q which make the computed solution match the y1 data as closely as possible."
We'll start by considering how to fit our model to the first set of data (the set given by [t, y1(t)] in the table). We will be following the steps given in the handout.
a.) Find the initial values of y and y'. The initial value of y is y(0) = .14 [get this from your actual data]. Estimate the initial velocity, v0, by taking the average velocity over the first time interval in the table of actual data. For the first set of data, y1(t), this estimate is
v0 = y'(0) ~= [y(2)-y(0)]/(2-0) = (.38-.14)/2 = .12
t y1(t) y2(t) 0 0.14 0.50 2 0.38 0.40 4 0.21 -0.12 6 -0.18 -0.49 8 -0.37 -0.33 10 -0.15 0.18 12 0.21 0.48 14 0.34 0.25 16 0.10 -0.23 18 -0.24 -0.45 20 -0.31 -0.18 22 -0.05 0.27 b.) Enter the actual data for time, y1, and y2 into vectors by entering these commands:
s = 0:2:22;
y1 = [.14, .38, .21, ..., -.05]
y2 = [.50, -.12, -.49, ..., .27]
And let's go ahead and plot the actual data for the first set of values:
cla;
plot(s, y1);
Continuing to follow the steps described in the handout, we notice that we will need to make an M-file so we can use it. This M-file is to model the system y'' + p*y' + q*y = k*y^3, which is equivalent to the system
y' = v
v' = -q*y - p*v + k*y^3.
To make an M-file that ODE23 will use to solve this differential equation, enter this for your M-file, yy.m:
-------------------------------------
function g = yy(t,x)
k = ...;
q = ...;
p = ...;
A = [0, 1; -q, -p];
g = A*x + [0; k*x(1)^3];
-------------------------------------
(Initially you will want to set p = k = 0, and use some random value for q, as per handout instructions.) Save this and make sure that it is in a directory in your Matlab path. But don't close this file. You're going to adjust the values of k, p, and q in this file many times.
Now use ode23 and this M-file to generate a vector [t,y]of a solutions to the ODE we're modeling. Here is the command we use:
[t, y] = ode23('yy', [0,22], [.14, .12]);
This command uses 0 <= t <= 22 and y(0) = .14, v0 = .12.
Now let's plot the ode23-computed solution, [t, y]:
hold on;
plot(t, y(:,1));
You might find that the initial velocity of the ode23 solution doesn't quite match the actual data. If you do, then just re-enter the "[t, y] = ode23(...)" command, but using your own value for v0 instead of 0.12. Then you would clear the graph and re-plot your actual values along with your new vector [t, y].
c.) Now we start adjust p, q, and k. Remember that p = k = 0 initially and
Step 1: In your m-file, adjust q (stiffness) until the period matches with the actual data.
Step 2: Readjust v0 until initial velocity looks correct.
Step 3: Adjust p (damping) until the amplitude is correct.
After making these adjustments, you should have a set of data that matches the actual data for y1. Record your values of p, q, and v0.
(2) "Using the same p and q, compare the computed solution to the y2 data. If the agreement is reasonably close, then leave k at zero. If the approximation is bad, try increasing k until a good approximation is obtained."
Now you want to model the y2 data. Keep p and q the same, but use values for y(0) and v0 that correspond to y2's data rather than y1's data. Note that for y2(0) = .50. An initial guess for v0 is (y2(2) - y2(0)) / 2 = (.5 - .4) / 2 = .05. You may want to adjust this value, though. Use these initial data and enter
cla;
[t, y] = ode23('yy', [0,22], [.50, your v0 goes here ]);
plot(s,y2);
plot(t,y);
Adjust v0 until it the initial velocity seems to match with the actual data. If the agreement is close, then leave k = 0. Otherwise, increase k until the approximation is good.
Regarding questions 1-3 from the first page:
- Which of the following equations forms a better model?
y'' + p*y' + q*y = 0,
y'' + p*y' + q*y = k*y^3
- Review the last page of the lab and / or use any mathematical reasoning you can to say which one is better.- Having decided which equation is better, determine approximate values of the parameters.
- You did this above.- Comment on how well the models with your choice of parameters agree with the data. Is there any danger that the building will fall down in a heavy wind storm?
- First, comment on how well your model agrees with the actual data.
- To see whether or not the building will collapse under your model, use the parameters you found and get the phase plane portrait of the system using pplane. Plot a bunch of orbits, and see what happens to the displacement of the building as t approaches infinity. Since y represents displacement, look at the behavior of y over time and give a reason why we can tell whether the building is falling over based on where y is going as time goes on.That should be it for this lab.
|
e-mail: renn@purdue.edu |