# MATH 2250 # PROJECT 5: EARTHQUAKES # APRIL 1999 # # Your final project for Math 2250 this semester is the Earthquake # project on pages 302-304 of Edwards-Penney. You are mostly on your # own for this project, but here is a small example of a spring system # worked out on Maple, so that you can get an idea about useful commands # to use. This file is available from the web, at the 2250 maple # projects home page: # http://www.math.utah.edu/~korevaar/2250spring99.html Remember, as # always, that after you save it from your browser to your own # directory you should then open it as a ``Maple Text'' document from # MAPLE. # # Let's start with example 1 on page 293 of Edwards-Penney. Initially # it is an unforced system with two masses and two springs, as you can # see from the description on page 293. We can write the system as # Mx''=Kx, where M is the ``mass matrix'', K is the ``spring matrix'', # and x is the displacement vector. Following the book's notation, we # enter > with(linalg):with(plots):with(DEtools): #tools for project > M:=matrix([[2,0],[0,1]]); > K:=matrix([[-150,50],[50,-50]]); > A:=evalm(inverse(M)&*K); [2 0] M := [ ] [0 1] [-150 50] K := [ ] [ 50 -50] [-75 25] A := [ ] [ 50 -50] Then the system can also be written as x''=Ax, and the eigenvectors of A determine fundamental modes, and the corresponding negative eigenvalues are the (opposites) of the squares of the corresponding angular frequencies: > eigenvects(A); [-25, 1, {[1, 2]}], [-100, 1, {[1, -1]}] # Therefore, the natural frequencies of this system are the 5 and 10, # and the two fundamental modes correspond to the masses moving in # parallel directions (with amplitude ratio of two and angular frequency # 5), and in opposite directions (with equal amplitudes and angular # frequency 10). # Now, let's consider the forced system with force vector equal to # cos(wt)[0,50], i.e. the second mass is being forced periodically. In # other words, the system Mx''=Kx + F, where F=cos(wt)[0,50] discussed # on page 299. We follow the method described at the top of page 299 to # find a particular solution to the forced oscillation problem. In # fact, this turns into example 3. > F0:=evalm(inverse(M)&*vector([0,50])); > #The F0 in the normalized equation (30), page 299 > Iden:=array(1..2,1..2,identity); > #the 2 by 2 identity matrix > Aleft:=omega->evalm(A + omega^2*Iden); > #the matrix function on the left side of (32) > c:=omega->evalm(-inverse(Aleft(omega))&*F0); > #the vector c(omega) in (32) F0 := [0, 50] Iden := array(identity, 1 .. 2, 1 .. 2, []) 2 Aleft := omega -> evalm(A + omega Iden) c := omega -> evalm(-(inverse(Aleft(omega)) `&*` F0)) > c(omega); #see equation (35) page 299 [ 2 ] [ 1 -75 + omega ] [1250 --------------------------, -50 --------------------------] [ 2 4 2 4] [ 2500 - 125 omega + omega 2500 - 125 omega + omega ] # The vector c(w) above, times the oscillation cos(wt), is a particular # solution to the forced oscillation problem we are considering. If we # assume that our actual problem has a small amount of damping, then we # expect that this particular solution is very close to the steady state # solution to the damped problem. See the dsiscussion on page 300. We # can study resonance phenomena for these slightly damped problems by # plotting the maximum amplitude of the (almost) steady state solutions. # Use ``norm'' to measure this maximum amplitude: > norm(c(omega)); 1 max(1250 ------------------------------, 2 4 | 2500 - 125 omega + omega | | 2 | | -75 + omega | 50 | -------------------------- |) | 2 4 | | 2500 - 125 omega + omega | > plot(norm(c(omega)),omega=0..15,y=0..15, > numpoints=200,color=`black`); # This is the picture on page 300. Notice the peaks at angular # frequency 5 and 10, corresponding to resonance with the two # fundamental modes. # # # # # # We can get a plot of resonance as a function of period by # recalling that 2*Pi/T=omega: > res:=T->norm(c(2*Pi/T)); > plot(res(T),T=0.1..3,y=0..15,numpoints=200,color=`black`); Pi res := T -> norm(c(2 ----)) T # # COMMENTS FOR THE EARTHQUAKE PROJECT: # (1) Students are often confused by the forcing term used on page 303, # namely > E*(omega)^2*cos(omega*t)*b; 2 E omega cos(omega t) b # where b is the vector [1,1,1,1,1,1,1]. They ask, ``how can the # earthquake be forcing all seven stories, it seems like it's just # shaking the bottom one.'' Well, the students are correct, but so is # Edwards-Penney. The authors talk about an ``opposite inertial force'' # being the reason for this forcing term and here's one way to think # about it. Maybe your instructor can help you more if it's still # confusing. Anyhow, think of the ground as the zeroth story. In the # rest frame it is shaking with oscillation Ecos(wt). And so its # acceleration is its second time derivative, namely -E*w^2*cos(wt). If # you write down the inhomogeneous system of EIGHT second order DE's for # the accelerations of stories zero thru seven, the forcing (well, # accelerating) term is -E*w^2*cos(wt)*[1,0,0,0,0,0,0,0], as you would # expect. Call the solution 8-vector to this system y(t), then see what # the shaking looks like to someone on the ground by letting # x(t)=y(t)-E*cos(wt)*[1,1,1,1,1,1,1,1]. Then the zeroth story # component of x(t) will be identically zero, and the other seven # components will satisfy equation (2) on the bottom of page 303, # exactly as the authors claim. # # (2) For large matrices the eigenvect command won't work well unless # you enter at least one decimal number; if all entries are rational # numbers (expressed without decimal points), Maple tries to find the # eigenvalues and eigenvectors algebraically and exactly, instead of # numerically, and often fails. Make sure at least one of your matrix # entries has a decimal point in it.