# Math 2250
# Earthquake project
# November 2001
# 
# Enter your name here: ?????????????????????????????
# 
# Project 3. Please fill in all areas below marked ????? and solve
# problems 3.1 to 3.6. The problem headers:
# 
#  _______     PROBLEM 3.1. BUILDING MODEL FOR AN EARTHQUAKE.
#  _______     PROBLEM 3.2. TABLE OF NATURAL FREQUENCIES AND PERIODS.
#  _______     PROBLEM 3.3. UNDETERMINED COEFFICIENTS STEADY-STATE 
#                            PERIODIC SOLUTION.
#  _______     PROBLEM 3.4. PRACTICAL RESONANCE.
#  _______     PROBLEM 3.5. EARTHQUAKE DAMAGE.
#  _______     PROBLEM 3.6. FIVE FLOORS.
# 
> restart:with(plots):with(linalg):
# 
# 3.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 3.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 3.1
> 
# 
# 3.2. TABLE OF NATURAL FREQUENCIES AND PERIODS.
# Refer to figure 7.4.17, page 437.
# 
# PROBLEM 3.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 .
# 
# # Sample code for a  4x3  table.
# # Use maple help to learn about nops, printf.
#  ev:=[-38.3,-33.4,-26.2,-17.9]: 
#  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
# ops(ev));
# 
> # Problem 3.2
> 
# 
# 3.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 contains the general
# solution to
# to homogeneous problem,  but rather the particular periodic solution
# with angular freqency w, 
# 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.  The idea is that in any real problem
# there would be a small amount
# of damping, which would lead to transient homogeneous solution, and
# steady periodic solution close to
# the one you are finding for this undamped problem.  
# 
# PROBLEM 3.3.
# Define  b:=(1/4)*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 evalm, diag, vector.
#   w:='w': 
#   b:=0.25*w*w*vector([1,1,1,1]):
#   A:=matrix([ [-20,10,0,0], [10,-20,10,0], [0,10,-20,10],
# [0,0,10,-10]]);
#   c:=linsolve(evalm(A+w*w*diag(1,1,1,1)),-b): 
#   evalf(c[1],6);
# 
> # PROBLEM 3.3
> 
# 
# 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 occurs when any component of  x(t)  has large
# amplitude compared to the vector
# norm of  the forcing amplitudes in 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 3.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, diag, vector, subs, evalm and
# linsolve.
#   with(linalg):
#   w:='w': Max:= c -> norm(c,infinity); 
#   B:=w*w*diag(1,1,1,1): b:=0.25*w*w*vector([1,1,1,1]):
#   C:=ww -> subs(w=ww,linsolve(evalm(A+B),-b)):
#   plot(Max(C(2*Pi/r)),r=0..6,0..10);
# 
> restart:with(plots):with(linalg):
> # PROBLEM 3.4. WARNING: Save your file often!!!
> 
# 
# 
# 3.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 3.5.
# Replot the amplitudes in 3.4 for periods 0 to 6 and amplitudes 5 to
# 10. There will be three
# spikes. Zoom-in on each spike, choosing a T-interval that shows the
# full spike. Determine
# from the three 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.4 to 1.5.
# plot(Max(C(2*Pi/r)),r=1.4..1.5,5..10);
# 
> # PROBLEM 3.5. WARNING: Save your file often!!
> # Plot three spikes on one graph.
> # Plot three zoom-in graphs.
> # Print period ranges.
# 
# 
# 3.6. FIVE FLOORS.
# Consider a building with five floors each weighing 20 tons. Assume
# each floor corresponds to a restoring
# Hooke's force with constant k=4 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 3.6.
# Model the 5-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.
# 
> restart:with(plots):with(linalg):
> # PROBLEM 3.6. WARNING: Save your file often!!
> # Define k=??? and m=???, then matrix A=???.
> 
> # Amplitude plot for T=0..6,C=4..10
> # Plot 4 zoom-in graphs
> # From the graphics, T=??..??, ??..??, ??..??, ??..??
> # give amplitude spikes above 4 feet. These are
> # determined by left mouse-clicks on the graph, so they
> # are approximate values only. 
> 
>