# 
# MATH 3150 Spring 99
# PROJECT 1
# HEAT FLOW AND FOURIER SERIES
# 
# 
# 
# This Maple project has four parts. 
# 
# A. Heat flow driven by temperature gradients
# B. Fourier sine series
# C. Fourier cosine series
# D. The complete Fourier Series.
# 
# Each Part has  two problems for you to do, or eight all together,
# numbered  A.1, A.2, B.1, B.2, C.1, C.2, D.1, D.2 .
# 
# This Maple work sheet has a lot of extra junk in it -- namely all my
# comments and explanations. 
# 
# DO  NOT HAND ALL THAT IN.
# 
# Instead,  just copy over the commands you want to another worksheet
# and do the problems there. Label each problem as you go, A.1 etc.,
# otherwise the grader may get lost and be unable to find your marks.
# 
# We need the packages  <plots> and <linalg>, so we load them up now. 
# The commands are ended with colons instead of semicolons in order to
# avoid a lot of extraneous output.
> with(plots):
> with(linalg):
# 
# A. Heat flow driven by temperature gradients
# 
# Recall that heat flows in the direction  in which the temperature
# decreases most rapidly at a rate proportional to the slope of the
# temperature in that direction.  If T (x,y) is the temperature,this
# translates into mathematics as F = -K *grad (T ) = K*grad (-T), where
# F is the heat flow field and K is the conductivity. Let us illustrate
# this with some plots of the heat flow superimposed on contour plots of
# the temperature. In these problems K is always equal to 1.
# 
# Example: T(x,y) =cos(Pi*y)*sin(2*Pi*x)
# 
# We set this up with T as an expression rather than a function because
# it is a bit more convenient to do it this way in Maple.  The Pi's are
# in there to make it look pretty in a domain which is a 1x1 box.
#  
> T:=cos(Pi*y)*sin(2*Pi*x);

                      T := cos(Pi y) sin(2 Pi x)

# The Maple command for grad is "grad" , with the list [x,y] denoting
# the variables with respect to which  the partials are taken.  Notice
# the - on the T.
> F:=grad(-T,[x,y]);

     F := [-2 cos(Pi y) cos(2 Pi x) Pi, sin(Pi y) Pi sin(2 Pi x)]

# We now make a contour plot of the temperature and a field plot (little
# arrows) of the heat flow  in a 1x1 box.  We follow the commands by :
# so they aren't  displayed.
> Tplot:=contourplot(T,x=0..1,y=0..1):
> Fplot:=fieldplot(F,x=0..1,y=0..1):
# Now we make the display with the two supperimposed  on each other.
> display([Tplot,Fplot]);

# There are hot spots in the lower left and upper right and cold spots
# in the upper left and lower right.  Notice how the heat flows downhill
# and perpindicular to the contours, although if the scales on the x and
# y axis are not the same the arrows may not appear perpindicular. Also
# note that the arrows on the top and bottom edges do not point across
# the  edges, i.e. , there is no heat flow across the edges, they are
# insulated. On the other hand, the left and right edges are at  0
# degrees, so heat flows from the edges to the cold spots, and from the
# hot spots to the edges.
# 
# Question :What kinds of boundary conditions (Dirichlet or Neumann) are
# satisfied on the four different edges?
# Answer:   top:         Neumann boundary condition with  dT/dn=0.
#                     bottom:  Neumann boundary condition with  dT/dn=0.
#                     left:          Dirichlet boundary condition with
# T=0.
#                     right:       Dirichlet boundary condition with
# T=0.
#  
# Answer:  top:           Dirichlet boundary condition with
# T=-sin(2*Pi*x)             
#                    bottom:   Dirichlet boundary condition with
# T=sin(2*Pi*x) 
#                     left:          Dirichlet boundary condition with
# T=0.
#                     right:       Dirichlet boundary condition with
# T=0.
# 
# Both of these answers are correct. The advantage of taking the first
# interpretation is that another solution that satisifies the same
# boundary conditions can be added to this one and the sum of the two
# solutions will still satisfy the boundary conditions because 0+0=0.
# This trivial fact is actually very important. 
# 
# If we impose the boundary conditions which T satisfies, then we can
# check if T is a steady state solution of the heat equation with these
# boundary conditions.   Recall  that  dT/dt=-div (F). Thus dT/dt = 0 
# if, and only if, div(F)=0.  So
# 
# Question: Is T a steady state solution of the heat equation with the
# boundary conditions found 
#                     in the first questions?   (The Maple function for
# div is diverge)
# 
> diverge(F,[x,y]);

                                               2
                     5 cos(Pi y) sin(2 Pi x) Pi

# Answer: No, div(F) is not identically zero.
# 
# OK, here are two similar  problems for you to do. Do the plots on
# x=0..1, y=0..1.
# 
# 
# 
# Problem A.1  T(x,y) = sinh(Pi*y)*sin(Pi*x)
# part 1. Do a field plot of the heat flow superimposed on a contour
# plot of the temperature.
# part 2. What kinds of boundary conditions (Dirichlet or Neumann) are
# satisfied on the four 
#              different edges? (Watch that top edge)
# part 3. Is T a steady state solution of the heat equation with the
# boundary conditions of part 2 ?
# 
# Problem A.2 T(x,y)=cosh(Pi*y)*cos(Pi*x)
# part 1. Do a field plot of the heat flow superimposed on a contour
# plot of the temperature.
# part 2. What kinds of boundary conditions (Dirichlet or Neumann) are
# satisfied on the four
#               different edges? (Watch that top edge)
# part 3. Is T a steady state solution of the heat equation with the
# boundary conditions of part 2 ?
# 
# B: Fourier sine series
# 
# All we are going to do here is compute some Fourier sine series and
# watch them converge.  In all of the problems L=1, so that the 
# expansion will be in terms of sin(n*Pi*x). Maple doesn't know n from
# spaghetti, so we have to tell it that n is an integer. That way 
# Maple will know sin(n*Pi) is  0.
> assume(n,integer);
# Now let us find the Fourier sine series for the function 
# 
# f(x) =    1       0    <=   x <  1/2
# f(x) = -1      1/2 <=   x <  1.
# 
# Recall that bn=(2/L) int(f(x)*sin(n*Pi*x/L),x=0..L). Because L=1 in
# these problems, this becomes
# bn=(2) int(f(x)*sin(n*Pi*x),x=0..1). 
# In order to compute the coefficients we need to split the integral
# into 2 parts.
> bn:=2*(int(1*sin(n*Pi*x),x=0..1/2)+int(-1*sin(n*Pi*x),x=1/2..1));# 

                                            n~
                cos(1/2 n~ Pi) - 1     -(-1)   + cos(1/2 n~ Pi)
       bn := -2 ------------------ - 2 ------------------------
                      n~ Pi                     n~ Pi

# Boy, is that ugly or what? Lets see if we can make it a bit prettier.
# 
> bn:=simplify(bn);

                                                    n~
                         2 cos(1/2 n~ Pi) - 1 - (-1)
                bn := -2 -----------------------------
                                     n~ Pi

# OK, that's a bit better. The twiddle (~) after the n marks n as an
# integer. It comes from the assume we did earlier. We really don't want
# bn as given because we can't vary n to use it in a sum. We want  b to
# be a function of n, i.e. , b(n). We do it this way:
> b:=unapply(bn,n);

                                                      n~
                           2 cos(1/2 n~ Pi) - 1 - (-1)
             b := n~ -> -2 -----------------------------
                                       n~ Pi

# For instance, to compute the first 10 coefficents,  just do:
> for i from 1 to 10 do b(i) od;

                                  0


                                   1
                                4 ----
                                   Pi


                                  0


                                  0


                                  0


                                    1
                               4/3 ----
                                    Pi


                                  0


                                  0


                                  0


                                    1
                               4/5 ----
                                    Pi

# OK, we want to see how well the sine series converges to the function.
# That is, we will approximate f by using the first n  terms in  the
# series. These approximations are called partial sums and depend on
# both  x and the number of terms  n; it is a function of  both x and n.
#  We do it this way 
> Psum:=(x,n)->sum(b(k)*sin(k*Pi*x),k=1..n);

                                   n
                                 -----
                                  \
               Psum := (x, n) ->   )   b(k) sin(k Pi x)
                                  /
                                 -----
                                 k = 1

# OK, we are not going to plot f(x) because it is a hassle , but if we
# take enough terms in the partial sum we should get a reasonable
# acsimile thereof.  So let's plot Psum with 2,10,22, and 102 terms.
> plot([Psum(x,2),Psum(x,10),Psum(x,22),Psum(x,102)],x=-1..2);

# Not too shoddy. By the way, I plotted from x= -1 to  x=2 (not from 0
# to 1) because I wanted you to see that the partial sums are odd
# functions of period 2. Here they are also of period 1, but this is a a
# special case.
# 
# Notice  the bat ears at the jumps. The super technical EE term for
# this is "ringing", but among mathematicians it goes by the name 
# "Gibb's phenomenon". There are ways to reduce  them , techniques of
# apodization - cutting off the feet - but that is not our concern here.
#  What you should notice is that the ringing 
# DOES    NOT     GO     TO     ZERO  
# as the number of terms increases. Yet at any fixed x (say x=.49999)
# the errors will get small as we take more and more terms. We say the
# series converges pointwise (i.e. at any point),  but, because the
# maximum error is always large, it does not converge uniformly.
# 
# OK, now do the same thing for the following functions. Remember that
# L=1, and experiment with the partial sums you plot in order to obtain
# pleasing results. (But don't show your experiments).
# 
# Problem B.1  
#    f(x) =    x                    0    <=   x <  3/4
#    f(x) = -3*(x-1)      3/4 <=   x <  1.
# 
# Problem B.2
#   f(x)=x^2                    0   <= x    < 1
# 
# C: Fourier cosine series.
# 
# Now we want to express f(x) as f(x) = a0 + sum( an*sin(n*Pi*x/L,
# n=1..infinity). Recall that 
# 
# a0=(1/L)*int(f(x),x=0..L)
# an=(2/L)*int(f(x)*cos(n*Pi*x/L),x=0..L)
# 
# We will use L=1 in all of the problems. The procedure is just as
# before, except for the  a0 term. Let's illustrate this for the
# function
#  
#  f(x) =  1       0    <=   x <  1/2
#  f(x) =  0      1/2 <=   x <  1.
> a0:=int(1,x=0..1/2);

                              a0 := 1/2

> an:=2*int(1*cos(n*Pi*x),x=0..1/2);

                                sin(1/2 n~ Pi)
                        an := 2 --------------
                                    n~ Pi

> a:=unapply(an,n);

                                  sin(1/2 n~ Pi)
                     a := n~ -> 2 --------------
                                      n~ Pi

> Psum:=(x,n)->a0+sum(a(k)*cos(k*Pi*x),k=1..n);

                                  /  n                   \
                                  |-----                 |
                                  | \                    |
           Psum := (x, n) -> a0 + |  )   a(k) cos(k Pi x)|
                                  | /                    |
                                  |-----                 |
                                  \k = 1                 /

> plot([Psum(x,1),Psum(x,5),Psum(x,21),Psum(x,101)],x=-1..2);

# 
# Got those same old rabbit ears at the breaks. Please notice that the
# series is an even function of period 2.
# 
# OK , do these sort of plots for the following. Remember that L=1.
# 
# Problem C.1  
#    f(x) =    x                    0    <=   x <  3/4
#    f(x) = -3*(x-1)      3/4 <=   x <  1.
# 
# Problem C.2
#   f(x)=x*(1-x)               0   <= x    < 1
# 
# 
# D. The complete Fourier Series.
# 
# Well, nothing left to do but put the previous two together.  We are
# looking at functions of period 2L whose Fourier series representation
# is of the form
# 
# f(x) = a0 + sum( an*cos(n*Pi*x) + bn*sin(n*Pi*x), n=1..infinity)
# 
# where
# a0 =(1/2L)* int(f(x), x=0..2*L)
# an =(2/2L)*int(f(x)*cos(n*Pi*x/L),x=0..2L)
# bn =(2/2L)*int(f(x)*sin(n*Pi*x/L),x=0..2L)
# 
# We will assume L=1 for all of these problems. Let us illustrate with
# the function
# 
# f(x) =  1       0    <=   x   <  1/2
# f(x) =  0      1 /2 <=   x   <  2
> a0:=(1/2)*int(1,x=0..1/2);

                              a0 := 1/4

> an:=int(1*cos(n*Pi*x),x=0..1/2);

                               sin(1/2 n~ Pi)
                         an := --------------
                                   n~ Pi

> bn:=int(1*sin(n*Pi*x),x=0..1/2);

                              cos(1/2 n~ Pi) - 1
                      bn := - ------------------
                                    n~ Pi

> a:=unapply(an,n);

                                 sin(1/2 n~ Pi)
                      a := n~ -> --------------
                                     n~ Pi

> b:=unapply(bn,n);

                                cos(1/2 n~ Pi) - 1
                   b := n~ -> - ------------------
                                      n~ Pi

> Psum:=(x,n)->a0+sum(a(k)*cos(k*Pi*x)+b(k)*sin(k*Pi*x),k=1..n);

  Psum :=

                       /  n                                        \
                       |-----                                      |
                       | \                                         |
        (x, n) -> a0 + |  )   (a(k) cos(k Pi x) + b(k) sin(k Pi x))|
                       | /                                         |
                       |-----                                      |
                       \k = 1                                      /

> plot([Psum(x,3),Psum(x,10),Psum(x,20),Psum(x,100)],x=-1..4);

# Got them old mouse ear blues.  Notice that the Fourier series
# representation has period 2, but it is neither even nor odd about 0.
# OK,  do the same with the following..
# 
# 
# Problem D.1
# 
# f(x) =  x*(1-x)       0    <=   x   <  1
# f(x) =  0                    1    <=   x   <  2
# 
# Problem D.2
# 
# f(x) = (x-1/2)^2   0  <=   x    <  2
# 
# FINIS