# Math 2250
# Maple Lab 7: Earthquake project
# Fall 2006
#
# Name _____________________________________ Class Time __________
#
# Project 7. Solve problems L7.1 to L7.6. The problem headers:
#
# _______ PROBLEM L7.1. BUILDING MODEL FOR AN EARTHQUAKE.
# _______ PROBLEM L7.2. TABLE OF NATURAL FREQUENCIES AND PERIODS.
# _______ PROBLEM L7.3. UNDETERMINED COEFFICIENTS STEADY-STATE
# PERIODIC SOLUTION.
# _______ PROBLEM L7.4. PRACTICAL RESONANCE.
# _______ PROBLEM L7.5. EARTHQUAKE DAMAGE.
# _______ PROBLEM L7.6. SIX FLOORS.
#
> with(linalg):
#
# 1. BUILDING MODEL FOR AN EARTHQUAKE.
# Refer to the textbook of Edwards-Penney, section 7.4, page 437.
# Consider a building with 7 floors.
#
# Let the mass in slugs of each story be m=1000.0 and let the spring
# constant be k=10000.0 (lbs/foot).
# Define the 7 by 7 mass matrix M and Hooke's matrix K for this system
# and convert Mx''=Kx into
# the system x''=Ax where A is defined by textbook equation (1) , page
# 437.
#
# PROBLEM L7.1
# Find the eigenvalues of the matrix A to six digits, using the Maple
# command "eigenvals(A)."
# Justify in particular that all seven eigenvalues are negative by
# direct computation.
#
# # Sample Maple code for a model with 4 floors.
# # Use maple help to learn about evalf and eigenvals.
# A:=matrix([ [-20,10,0,0], [10,-20,10,0], [0,10,-20,10],
# [0,0,10,-10]]);
# evalf(eigenvals(A));
#
> # Problem L7.1
>
#
# 2. TABLE OF NATURAL FREQUENCIES AND PERIODS.
# Refer to figure 7.4.17, page 437.
#
# PROBLEM L7.2.
# Find the natural angular frequencies omega=sqrt(-lambda) for the
# seven story building and also the corresponding periods
# 2PI/omega, accurate to six digits. Display the answers in a table .
# The answers appear in Figure 7.4.17, page 437.
#
# # Sample code for a 4x3 table.
# # Use maple help to learn about nops and printf.
# ev:=[-10,-1.206147582,-35.32088886,-23.47296354]: n:=nops(ev):
# Omega:=lambda -> sqrt(-lambda):
# format:="%10.6f %10.6f %10.6f\n":
# seq(printf(format,ev[i],Omega(ev[i]),2*evalf(Pi)/Omega(ev[i])),i=1..n
# );
#
> # Problem L7.2
>
#
# 3. UNDETERMINED COEFFICIENTS STEADY-STATE PERIODIC SOLUTION.
# Consider the forced equation x'=Ax+cos(wt)b where b is a constant
# vector. The earthquake's
# ground vibration is accounted for by the extra term cos(wt)b, which
# has period T=2Pi/w.
# The solution x(t) is the 7-vector of excursions from equilibrium
# of the corresponding 7 floors.
# Sought here is not the general solution, which certainly contains
# transient terms, but rather the
# steady-state periodic solution, which is known from the theory to have
# the form x(t)=cos(wt)c
# for some vector c that depends only on A and b.
#
# PROBLEM L7.3.
# Define b:=0.25*w*w*vector([1,1,1,1,1,1,1]): in Maple and find the
# vector c in the undetermined
# coefficients solution x(t)=cos(wt)c. Vector c depends on w. As
# outlined in the textbook, vector c
# can be found by solving the linear algebra problem -w^2 c = Ac + b;
# see page 433. Don't print c,
# as it is too complex; instead, print c[1] as an illustration.
#
# # Sample code for defining b and A, then solving for c in the
# 4-floor case.
# # See maple help to learn about vector and linsolve.
# w:='w': u:=w*w: b:=0.25*u*vector([1,1,1,1]):
# Au:=matrix([ [-20+u,10,0,0], [10,-20+u,10,0], [0,10,-20+u,10],
# [0,0,10,-10+u]]);
# c:=linsolve(Au,-b):
# evalf(c[1],2);
#
> # PROBLEM L7.3
>
#
# 4 PRACTICAL RESONANCE.
# Consider the forced equation x'=Ax+cos(wt)b of 3.3 above with
# b:=0.25*w*w*vector([1,1,1,1,1,1,1]).
# Practical resonance can occur if a component of x(t) has large
# amplitude compared to the vector
# norm of b. For example, an earthquake might cause a small 3-inch
# excursion on level ground, but
# the building's floors might have 50-inch excursions, enough to destroy
# the building.
#
# PROBLEM L7.4.
# Let Max(c) denote the maximum modulus of the components of vector c.
# Plot g(T)=Max(c(w)) with w=(2*Pi)/T
# for periods T=0 to T=6, ordinates Max=0 to Max=10, the vector c(w)
# being the answer produced in 3.3 above.
# Compare your figure to the textbook Figure 7.4.18, page 438.
#
# # Sample maple code to define the function Max(c), 4-floor building.
# # Use maple help to learn about norm, vector, subs and linsolve.
# with(linalg):
# w:='w': Max:= c -> norm(c,infinity); u:=w*w:
# b:=0.25*w*w*vector([1,1,1,1]):
# Au:=matrix([ [-20+u,10,0,0], [10,-20+u,10,0], [0,10,-20+u,10],
# [0,0,10,-10+u]]);
# C:=ww -> subs(w=ww,linsolve(Au,-b)):
# plot(Max(C(2*Pi/r)),r=0..6,0..10,numpoints=150);
#
> # PROBLEM L7.4. WARNING: Save your file often!!!
>
#
#
# 5. EARTHQUAKE DAMAGE.
# The maximum amplitude plot of 3.4 can be used to detect the likelihood
# of earthquake damage for a given
# ground vibration of period T. A ground vibration (1/4)cos(wt),
# T=2*Pi/w, will be assumed, as in 3.4.
#
# PROBLEM L7.5.
# (a) Replot the amplitudes in 3.4 for periods 1.14 to 4 and amplitudes
# 5 to 10. There will be four spikes.
# (b) Create four zoom-in plots, one for each spike, choosing a
# T-interval that shows the full spike.
# (c) Determine from the four zoom-in plots approximate intervals for
# the period T such that
# some floor in the building will undergo excursions from equilibrium in
# excess of 5 feet.
#
# # Example: Zoom-in on a spike for amplitudes 5 feet to 10 feet,
# periods 1.97 to 2.01.
# with(linalg): w:='w': Max:= c -> norm(c,infinity); u:=w*w:
# Au:=matrix([ [-20+u,10,0,0], [10,-20+u,10,0], [0,10,-20+u,10],
# [0,0,10,-10+u]]);
# b:=0.25*w*w*vector([1,1,1,1]):
# C:=ww -> subs(w=ww,linsolve(Au,-b)):
# plot(Max(C(2*Pi/r)),r=1.97..2,01,5..10,numpoints=150);
#
> # PROBLEM L7.5. WARNING: Save your file often!!
> #(a) Plot four spikes on separate graphs
> #(b) Plot four zoom-in graphs.
> #(c) Print period ranges.
>
#
#
# 6. SIX FLOORS.
# Consider a building with six floors each weighing 50 tons. Assume each
# floor corresponds to a restoring
# Hooke's force with constant k=5 tons/foot. Assume that ground
# vibrations from the earthquake are modeled by
# (1/4)cos(wt) with period T=2*Pi/w (same as the 7-floor model above).
#
# PROBLEM L7.6.
# Model the 6-floor problem in Maple. Plot the maximum amplitudes
# against the period 0 to 6 and amplitude
# 4 to 10. Determine from the graphic the period ranges which cause the
# amplitude plot to spike above 4 feet.
# Sanity check: m=3125, and the 6x6 matrix contains fraction 16/5.
# There are five spikes.
#
> # PROBLEM L7.6. WARING: Save your file often!!
> # Define k, m and the 6x6 matrix A.
> # Amplitude plot for T=0..6,C=4..10
> # Plot five zoom-in graphs
> # From the graphics, five T-ranges give amplitude
> # spikes above 4 feet. These are determined by left
> # mouse-clicks on the graph, so they are approximate values only.
>