# Math 2250
# Earthquake project
# April 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. THREE FLOORS.
# 
> with(DEtools):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 augment, stackmatrix, diag, transpose,
# evalm.
#  A1:=diag(-20,-20,-20,-10):
#  A2:=augment(vector([0,0,0,0]),stackmatrix(diag(1,1,1),matrix([[0,0,0]
# ]))):
#  A3:=transpose(A2):
#  A:=evalm(A1+10*A2+10*A3);
#  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 printf.
#  ev:=[-38.3,-33.4,-26.2,-17.9]:  n:=4:
#  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 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 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 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.
#   w:='w': 
#   b:=0.25*w*w*vector([1,1,1,1]):
#   A1:=diag(-20,-20,-20,-10): 
#  
# A2:=augment(vector([0,0,0,0]),stackmatrix(diag(1,1,1),matrix([[0,0,0]]
# ))):
#   A:=evalm(A1+10*A2+10*transpose(A2));
#   c:=linsolve(evalm(A+w*w*diag(1,1,1,1)),-b): 
#   evalf(c[1]);
# 
> # 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 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 3.4.
# Let Max(c) denote the maximum modulus of the components of vector c.
# Plot  g(T)=Max(c(w)) with w=T/(2*Pi)
# 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."
#   with(linalg):
#   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);
# 
> # PROBLEM 3.4
# 
# 
# 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. Determine
# approximate ranges for the period  T  such that some floor in the
# building will
# undergo excursions from equilibrium in excess of 5 feet.
# 
> # PROBLEM 3.5
> # Plot it.
> # Print period ranges ???.
# 
# 
# 3.6. THREE FLOORS.
# Consider a building with only three 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 3-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.
# 
> # PROBLEM 3.6
> # Define k=??? and m=???, then matrix A=???.
> # Amplitude plot for T=0..6,C=4..10
> # From the graphics, T-ranges are ???? 
> # which give amplitude spikes above 4 feet. These are
> # determined by mouse-clicks on the graph, so they
> # are approximate values only. 
>