# Math 2250
# Maple Lab 8: Earthquake project
# F2008
# 
# Name _____________________________________   Class Time __________
# 
# Project 8.  Solve problems L8.1 to L8.5. The problem headers:
# 
# _______     PROBLEM L8.1. EARHQUAKE MODEL FOR A BUILDING.
# _______     PROBLEM L8.2. TABLE OF NATURAL FREQUENCIES AND PERIODS.
# _______     PROBLEM L8.3. UNDETERMINED COEFFICIENTS STEADY-STATE SOL
# _______     PROBLEM L8.4. PRACTICAL RESONANCE.
# _______     PROBLEM L8.5. EARTHQUAKE DAMAGE.
# 
# SIX FLOOR Model.
# Refer to the textbook of Edwards-Penney, section 7.4, page 437.
# Consider a building with six floors each weighing 50 tons. 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.
# 
> 
# PROBLEM L8.1. BUILDING MODEL FOR AN EARTHQUAKE.
# Model the 6-floor problem in Maple. 
# Define the 6 by 6 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.
# 
# Sanity check:  Mass m=3125, and the 6x6 matrix contains fraction 16/5.
# 
# Then find the eigenvalues of the matrix A to six digits, using the 
# Maple command "eigenvals(A)."
# Sanity check: All six eigenvalues should be negative.
# 
> # 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]]);
> # with(linalg): evalf(eigenvals(A));
> 
> 
> # Problem L8.1
> # Define k, m and the 6x6 matrix A.
> # with(linalg): evalf(eigenvals(A));
> 
# PROBLEM L8.2. TABLE OF NATURAL FREQUENCIES AND PERIODS.
# Refer to figure 7.4.17, page 437.
# 
# Find the natural angular frequencies  omega=sqrt(-lambda) for the
# six story building and also the corresponding periods
# 2PI/omega, accurate to six digits.  Display the answers in a table .
# Compare with answers in Figure 7.4.17, page 437, for the 7-story case.
# 
> # Sample code for a 4x3 table, 4-story building.
> # 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 L8.2
> # ev:=[fill this in]: 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  L8.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  6-vector of excursions from equilibrium
# of the corresponding 6 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.
# 
# Define  b:=0.25*w*w*vector([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]):
> #  A:=matrix([ [-20,10,0,0], [10,-20,10,0],
> [0,10,-20,10],[0,0,10,-10]]);
> # Au:=evalm(A+u*diag(1,1,1,1));
> #  c:=linsolve(Au,-b): 
> #  evalf(c[1],2);
> 
> 
> # PROBLEM L8.3
> # Define w, u, b, A, Au, c  
> # evalf(c[1],2);
> 
# PROBLEM L8.4. PRACTICAL RESONANCE.
# Consider the forced equation  x'=Ax+cos(wt)b  of L8.3 above with
# b:=0.25*w*w*vector([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.
# 
# 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 L8.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]):
> # A:=matrix([ [-20,10,0,0], [10,-20,10,0], [0,10,-20,10],
> [0,0,10,-10]]);
> #  Au:=evalm(A+u*diag(1,1,1,1));
> #  C:=ww -> subs(w=ww,linsolve(Au,-b)):
> #  plot(Max(C(2*Pi/r)),r=0..6,0..10,numpoints=150);
> 
> # PROBLEM L8.4. WARNING: Save your file often!!!
> # w:='w': Max:= c -> norm(c,infinity): u:=w*w:   
> # Define b
> # Define A
> # Define Au
> # Define C  
> # plot(Max(C(2*Pi/r)),r=0..6,0..10,numpoints=150);
> 
> 
# 
# PROBLEM L8.5. EARTHQUAKE DAMAGE.
# The maximum amplitude plot of L8.4 can be used to detect the 
# 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 L8.4.
# 
# (a) Replot the amplitudes in L8.4 for periods 1.5 to 5.5 and
# amplitudes 5 to 10. 
# There will be five spikes. 
# (b) Create five zoom-in plots, one for each spike, choosing a
# T-interval that shows the full spike. 
# (c) Determine from the five 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 L8.5. WARNING: Save your file often!!
> #(a) Re-plot the five spikes.
> # plot(Max(C(2*Pi/r)),r=1.5..5.5,5..10,numpoints=150);
> #(b) Plot five zoom-in graphs.
> # one:=1.79..1.83:plot(Max(C(2*Pi/r)),r=one,5..10,numpoints=150);
> # two:=???:plot(Max(C(2*Pi/r)),r=two,5..10,numpoints=150);
> # three:=???:plot(Max(C(2*Pi/r)),r=three,5..10,numpoints=150);
> # four:=???:plot(Max(C(2*Pi/r)),r=four,5..10,numpoints=150);
> # five:=???:plot(Max(C(2*Pi/r)),r=five,5..10,numpoints=150);
> #(c) Print period ranges.
> # PeriodRanges:=[one,two,three,four,five];
>