(*^
::[	Information =

	"This is a Mathematica Notebook file.  It contains ASCII text, and can be
	transferred by email, ftp, or other text-file transfer utility.  It should
	be read or edited using a copy of Mathematica or MathReader.  If you 
	received this as email, use your mail application or copy/paste to save 
	everything from the line containing (*^ down to the line containing ^*)
	into a plain text file.  On some systems you may have to give the file a 
	name ending with ".ma" to allow Mathematica to recognize it as a Notebook.
	The line below identifies what version of Mathematica created this file,
	but it can be opened using any other version as well.";

	FrontEndVersion = "Macintosh Mathematica Notebook Front End Version 2.2";

	MacintoshStandardFontEncoding; 
	
	fontset = title, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeTitle, center, M7, bold, e8,  24, "Times"; 
	fontset = subtitle, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeTitle, center, M7, bold, e6,  18, "Times"; 
	fontset = subsubtitle, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeTitle, center, M7, italic, e6,  14, "Times"; 
	fontset = section, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, grayBox, M22, bold, a20,  18, "Times"; 
	fontset = subsection, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, blackBox, M19, bold, a15,  14, "Times"; 
	fontset = subsubsection, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, whiteBox, M18, bold, a12,  12, "Times"; 
	fontset = text, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  12, "Times"; 
	fontset = smalltext, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  10, "Times"; 
	fontset = input, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeInput, M42, N23, bold, L-4,  12, "Courier"; 
	fontset = output, output, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, L-4,  12, "Courier"; 
	fontset = message, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, R32768, L-4,  12, "Courier"; 
	fontset = print, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, L-4,  12, "Courier"; 
	fontset = info, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, B32768, L-4,  12, "Courier"; 
	fontset = postscript, PostScript, formatAsPostScript, output, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeGraphics, M7, l34, w282, h287,  12, "Courier"; 
	fontset = name, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, italic,  10, "Geneva"; 
	fontset = header, inactive, noKeepOnOnePage, preserveAspect, M7,  12, "Times"; 
	fontset = leftheader, inactive, L2,  12, "Times"; 
	fontset = footer, inactive, noKeepOnOnePage, preserveAspect, center, M7,  12, "Times"; 
	fontset = leftfooter, inactive, L2,  12, "Times"; 
	fontset = help, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  10, "Times"; 
	fontset = clipboard, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  12, "Times"; 
	fontset = completions, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  12, "Times"; 
	fontset = special1, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  12, "Times"; 
	fontset = special2, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  12, "Times"; 
	fontset = special3, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  12, "Times"; 
	fontset = special4, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  12, "Times"; 
	fontset = special5, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  12, "Times"; 
	paletteColors = 128; currentKernel; 
]
:[font = title; inactive; preserveAspect; startGroup]
Linear Methods of Applied Mathematics

Solving the Heat Equation
;[s]
4:0,1;38,2;39,1;64,0;65,-1;
3:1,25,18,Times,1,24,0,0,0;2,24,17,Times,1,23,0,0,0;1,12,8,Times,1,11,0,0,0;
:[font = text; inactive; preserveAspect; plain; bold]
(c) Copyright 1994-1997 by Evans M. Harrell II and James V. Herod. 
     All rights reserved.
:[font = text; inactive; Cclosed; preserveAspect; plain; bold; startGroup]
Notes for the instructor.
;[s]
2:0,1;25,0;26,-1;
2:1,13,9,Times,1,12,0,0,0;1,11,8,Times,1,9,0,0,0;
:[font = text; inactive; preserveAspect; fontSize = 14; endGroup]
This contains calculations and examples which correlate with chapter 9 of the WWW text by Harrell and Herod.  Students can be encouraged to cut and paste from this notebook to do homework.
:[font = subsubsection; inactive; preserveAspect; startGroup]
Instructions
:[font = text; inactive; preserveAspect; endGroup]
This notebook uses Mathematica  to perform calculations for Harrell and Herod's hypertext book, Linear Methods of Applied Mathematics.   The student needs only a basic knowledge of Mathematica  to use the notebook, which is designed both to show how to work problems in the text and to provide a template for using Mathematica  to work other problems of the student's own design.  

Calculations will be performed when the reader presses enter in a given calculation cell (bold print).  It is best to activate the cells in order, so that Mathematica  will be able to call on operators and functions defined in earlier cells.  Red color coding is used to warn the reader when a given calculation relies on an earlier one, and magenta color coding is used to warn the reader when a calculation requires Mathematica's vector analysis package to be loaded.
;[s]
17:0,0;19,1;30,0;96,1;133,0;181,1;192,0;315,1;326,0;538,1;549,0;626,2;642,0;725,3;745,0;801,1;812,0;853,-1;
4:9,13,9,Times,0,12,0,0,0;6,13,9,Times,2,12,0,0,0;1,13,9,Times,0,12,65535,0,0;1,13,9,Times,0,12,65535,0,65535;
:[font = text; inactive; preserveAspect; fontSize = 14]
A useful substitution, which we shall often make, is:
:[font = input; preserveAspect]
TrigId = {Cos[Pi n_] -> (-1)^n, Sin[Pi n_] -> 0};
:[font = text; inactive; preserveAspect; fontSize = 14]
Please activate this cell before proceeding.  Also, for the multidimensional calculations we shall invoke the vector calculus analysis with:
:[font = input; preserveAspect]
Needs["Calculus`VectorAnalysis`"]
:[font = section; inactive; Cclosed; preserveAspect; startGroup]
The Laplace Operator
:[font = text; inactive; preserveAspect; fontSize = 14]
The Laplace operator and other vector derivatives are recognized by Mathematica's vector analysis package.  The creators of Mathematica  have chosen to designate the Laplace operator as Laplacian (though many writers prefer the spelling  Laplacean).
;[s]
15:0,0;68,1;79,0;124,1;135,0;185,2;186,3;192,4;193,3;195,0;237,2;238,3;244,4;245,3;247,0;250,-1;
5:5,16,12,Times,0,14,0,0,0;2,16,12,Times,2,14,0,0,0;2,14,10,Times,0,13,0,0,0;4,13,10,Courier,0,12,0,0,0;2,13,10,Courier,1,12,0,0,0;
:[font = text; inactive; preserveAspect; fontSize = 14]
Mathematica  knows the definition of the Laplace operator symbolically as well as operationally:
;[s]
2:0,1;11,0;97,-1;
2:1,16,12,Times,0,14,0,0,0;1,16,12,Times,2,14,0,0,0;
:[font = input; preserveAspect; startGroup]
Laplacian[f[x,y,z]] == Div[Grad[f[x,y,z]]]
;[s]
6:0,1;9,0;23,1;26,0;27,1;31,0;43,-1;
2:3,12,10,Courier,1,12,0,0,0;3,12,10,Courier,1,12,65535,0,65535;
:[font = output; output; inactive; preserveAspect; endGroup]
True
;[o]
True
:[font = text; inactive; preserveAspect; fontSize = 14]
Unless and until we instruct Mathematica  otherwise, Cartesian variables named {x,y,z} are assumed.
;[s]
3:0,0;29,1;40,0;100,-1;
2:2,16,12,Times,0,14,0,0,0;1,16,12,Times,2,14,0,0,0;
:[font = subsection; inactive; Cclosed; preserveAspect; startGroup]
Examples
:[font = input; preserveAspect]
Grad[x y/z]
;[s]
2:0,1;4,0;12,-1;
2:1,12,10,Courier,1,12,0,0,0;1,12,10,Courier,1,12,65535,0,65535;
:[font = input; preserveAspect]
Div[{y/z, x/z, -((x y)/z^2)}]
;[s]
2:0,1;3,0;30,-1;
2:1,12,10,Courier,1,12,0,0,0;1,12,10,Courier,1,12,65535,0,65535;
:[font = input; preserveAspect]
Curl[{y/z, x/z, -((x y)/z^2)}]
;[s]
2:0,1;4,0;31,-1;
2:1,12,10,Courier,1,12,0,0,0;1,12,10,Courier,1,12,65535,0,65535;
:[font = text; inactive; preserveAspect; fontSize = 14]
(The curl will not be important in this chapter.)
:[font = input; preserveAspect; endGroup; endGroup]
Laplacian[x y/z]
;[s]
2:0,1;9,0;17,-1;
2:1,12,10,Courier,1,12,0,0,0;1,12,10,Courier,1,12,65535,0,65535;
:[font = section; inactive; preserveAspect]
Separation of variables and examples
:[font = text; inactive; preserveAspect]
As in the "flat" text, we separate variables for a rectangle with homogeneous boundary conditions on three sides.  Suppose that 0 < x < a and 0 < y < b.  On this rectangle, u[x,y] solves Laplace's equation in two variables.  
;[s]
1:0,0;226,-1;
1:1,16,12,Times,0,14,0,0,0;
:[font = text; inactive; preserveAspect]
We can define Laplace's equation as a logical operator, in case we simply wish to test whether a function is harmonic (= solution of Laplace's equation):
;[s]
1:0,1;154,-1;
2:0,13,9,Times,0,12,0,0,0;1,16,12,Times,0,14,0,0,0;
:[font = input; preserveAspect]
IsItHarmonic[u_] := Laplacian[u] == 0
;[s]
3:0,0;20,1;29,0;38,-1;
2:2,12,10,Courier,1,12,0,0,0;1,12,10,Courier,1,12,65535,0,65535;
:[font = subsection; inactive; Cclosed; preserveAspect; startGroup]
Examples
:[font = input; preserveAspect; startGroup]
IsItHarmonic[x^2 + y]
;[s]
2:0,1;12,0;22,-1;
2:1,12,10,Courier,1,12,0,0,0;1,12,10,Courier,1,12,65535,0,65535;
:[font = output; output; inactive; preserveAspect; endGroup]
False
;[o]
False
:[font = input; preserveAspect; startGroup]
IsItHarmonic[Sin[x]Sinh[y] + 5]
;[s]
2:0,1;12,0;32,-1;
2:1,12,10,Courier,1,12,0,0,0;1,12,10,Courier,1,12,65535,0,65535;
:[font = output; output; inactive; preserveAspect; endGroup; endGroup]
True
;[o]
True
:[font = subsection; inactive; preserveAspect]
Separation of variables 
:[font = text; inactive; preserveAspect; fontSize = 14]
Separation of variables works much as it did for the heat equation and the wave equation.  As with our earlier equations, we begin with the ansatz  that u is a product solution 
                 u[x,y] = X[x] Y[y]
and substitute into Laplace's equation; it is again convenient to divide the result by u.  We use first information from the PDE, then the information from homogeneous boundary conditions, and lastly information from inhomogeneous boundary conditions, which are treated much like the initial conditions for the earlier PDE's.
;[s]
3:0,0;140,1;146,0;540,-1;
2:2,16,12,Times,0,14,0,0,0;1,16,12,Times,2,14,0,0,0;
:[font = input; preserveAspect; startGroup]
Simplify[Laplacian[X[x] Y[y]]/(X[x] Y[y]) == 0]
;[s]
3:0,0;9,1;18,0;48,-1;
2:2,12,10,Courier,1,12,0,0,0;1,12,10,Courier,1,12,65535,0,65535;
:[font = output; output; inactive; preserveAspect; endGroup]
Derivative[2][X][x]/X[x] + Derivative[2][Y][y]/Y[y] == 0
;[o]
X''[x]   Y''[y]
------ + ------ == 0
 X[x]     Y[y]
:[font = text; inactive; preserveAspect; fontSize = 14]
Evidently, if the y-independent quantity X''/X equals the x-independent quantity -Y''/Y, both must be a constant, and we have the familiar ordinary differential equations,
	      X''  =  mu X
	      -Y'' =  mu Y

(notice the different sign).

Rather than tackling non-homogeneous boundary conditions on all four edges of a rectangle all at once, we begin by setting three of the four boundary functions to 0.  For definiteness, we also make a specific choice of the fourth while developing the ideas.  The specific example we analyze is:
;[s]
3:0,0;172,1;211,0;538,-1;
2:2,16,12,Times,0,14,0,0,0;1,15,11,Courier,0,14,0,0,0;
:[font = subsection; inactive; preserveAspect; startGroup]
Model Problem IX.1. Let us solve Laplace's equation with boundary conditions which are homogeneous on three of the four sides of a rectangle:

	u(0,y) = 0								 (9.7)
	u(1,y) = 1								 (9.8) 	
	u(x,0) = 0								 (9.9) 	
	u(x,2) = 0.								(9.10)
:[font = text; inactive; preserveAspect; fontSize = 14]
A good rule of thumb in this subject is to deal with the most homogeneous boundaries first.  It is the variable y which has two homogeneous boundary conditions in this case.  It satisfies essentially the same eigenvalue problem as we have seen in previous chapters:
	      -Y'' =  mu Y
	      Y[0] = 0, Y[2] = 0

;[s]
5:0,0;43,1;90,0;266,2;311,0;313,-1;
3:3,16,12,Times,0,14,0,0,0;1,16,12,Times,2,14,0,0,0;1,15,11,Courier,0,14,0,0,0;
:[font = input; dontPreserveAspect]
 HeatOp[f_] := HeatOp[f,1]
;[s]
3:0,0;15,1;21,0;27,-1;
2:2,12,10,Courier,1,12,0,0,0;1,12,10,Courier,1,12,65535,0,0;
:[font = text; inactive; Cclosed; preserveAspect; fontSize = 14; startGroup]
In this cell you may try out the wave operator.
:[font = input; dontPreserveAspect; fontSize = 14; startGroup]
HeatOp[t^3 Sin[x]^2]
;[s]
2:0,1;6,0;21,-1;
2:1,14,11,Courier,1,14,0,0,0;1,14,11,Courier,1,14,65535,0,0;
:[font = output; output; inactive; preserveAspect; fontSize = 14; endGroup]
-2*t^3*Cos[x]^2 + 3*t^2*Sin[x]^2 + 2*t^3*Sin[x]^2
;[o]
    3       2      2       2      3       2
-2 t  Cos[x]  + 3 t  Sin[x]  + 2 t  Sin[x]
:[font = text; inactive; preserveAspect; fontSize = 14]
That function was not a solution.  The following one is:
:[font = input; dontPreserveAspect; fontSize = 14; startGroup]
HeatOp[Exp[-t] Sin[x]] == 0
;[s]
2:0,1;6,0;28,-1;
2:1,14,11,Courier,1,14,0,0,0;1,14,11,Courier,1,14,65535,0,0;
:[font = output; output; inactive; preserveAspect; fontSize = 14; endGroup; endGroup; endGroup]
True
;[o]
True
:[font = section; inactive; Cclosed; preserveAspect; startGroup]
Solving the heat equation with the method of separation of variables - derivation
:[font = text; inactive; dontPreserveAspect; fontSize = 14]
To simplify matters a bit, we assume that k = 1 and x ranges from 0 to 1.  The more general solution can be obtained by scaling and shifting the variables t and x.  

We try to solve the wave equation with the ansatz that u[t,x] is of the special form T[t] X[x], known as a product solution.
:[font = input; preserveAspect]
u[t_,x_] := T[t] X[x]
:[font = input; preserveAspect; startGroup]
HeatOp[u[t,x]]
;[s]
2:0,1;6,0;15,-1;
2:1,12,10,Courier,1,12,0,0,0;1,12,10,Courier,1,12,65535,0,0;
:[font = output; output; inactive; dontPreserveAspect; endGroup]
X[x]*Derivative[1][T][t] - T[t]*Derivative[2][X][x]
;[o]
X[x] T'[t] - T[t] X''[x]
:[font = text; inactive; dontPreserveAspect; fontSize = 14]

This simplifies if we divide through by T[t] X[x]:
:[font = input; preserveAspect; startGroup]
%/(T[t] X[x])
:[font = output; output; inactive; preserveAspect; endGroup]
(X[x]*Derivative[1][T][t] - T[t]*Derivative[2][X][x])/
 
  (T[t]*X[x])
;[o]
X[x] T'[t] - T[t] X''[x]
------------------------
       T[t] X[x]
:[font = input; dontPreserveAspect; startGroup]
Simplify[%]
:[font = output; output; inactive; dontPreserveAspect; endGroup]
Derivative[1][T][t]/T[t] - Derivative[2][X][x]/X[x]
;[o]
T'[t]   X''[x]
----- - ------
T[t]     X[x]
:[font = text; inactive; preserveAspect; fontSize = 14]
If u solves the wave equation, this expression must be zero. (We could quibble about what happens when T or X = 0, but that would either mean that u is the uninteresting "trivial" solution or else the problem occurs only at isolated values of x and t.)  Thus 
                        T'/T = X''/X.
Since the left side is independent of x and the right side is independent of t, both sides must in fact equal a constant. We don't know its value yet, so let us call it - mu. (It turns out to be negative.)
Thus X''[x] = - mu X[x]:
:[font = input; preserveAspect; startGroup]
 DSolve[- X''[x]== mu X[x],X[x],x]
:[font = output; output; inactive; dontPreserveAspect; endGroup]
{{X[x] -> E^(-I*mu^(1/2)*x)*C[1] + E^(I*mu^(1/2)*x)*C[2]}}
;[o]
           -I Sqrt[mu] x         I Sqrt[mu] x
{{X[x] -> E              C[1] + E             C[2]}}
:[font = text; inactive; dontPreserveAspect; fontSize = 14; startGroup]

We now remind Mathematica about Euler's interesting formula for complex exponentials:
:[font = input; preserveAspect; startGroup]
 % /. {E^(I Sqrt[mu] x) -> Cos[Sqrt[mu] x] + 
               I Sin[Sqrt[mu] x],
         E^(- I Sqrt[mu] x) -> Cos[Sqrt[mu] x] -
                  I Sin[Sqrt[mu] x]}
:[font = output; output; inactive; preserveAspect; endGroup; endGroup]
{{X[x] -> C[1]*(Cos[mu^(1/2)*x] - I*Sin[mu^(1/2)*x]) + 
 
     C[2]*(Cos[mu^(1/2)*x] + I*Sin[mu^(1/2)*x])}}
;[o]
{{X[x] -> C[1] (Cos[Sqrt[mu] x] - I Sin[Sqrt[mu] x]) + 
 
     C[2] (Cos[Sqrt[mu] x] + I Sin[Sqrt[mu] x])}}
:[font = text; inactive; preserveAspect; fontSize = 14]

Evidently, we could reexpress this as 
:[font = input; preserveAspect]
Clear[X,D1,D2]
:[font = input; preserveAspect]
X[x_] := D1 Cos[Sqrt[mu] x] + 
        D1 Sin[Sqrt[mu] x]
:[font = input; preserveAspect; startGroup]
X[0]
:[font = output; output; inactive; preserveAspect; endGroup]
D1
;[o]
D1
:[font = text; inactive; preserveAspect; fontSize = 14]

This means that in order to satisfy the BC at x=0, we must have:
:[font = input; preserveAspect]
X[x_] := Sin[Sqrt[mu] x]
:[font = input; preserveAspect; startGroup]
Solve[X[L] == 0, mu]
:[font = message; inactive; preserveAspect]
Solve::ifun: 
   Warning: Inverse functions are being used by Solve, so
     some solutions may not be found.
:[font = output; output; inactive; preserveAspect; endGroup]
{{mu -> 0}}
;[o]
{{mu -> 0}}
:[font = input; preserveAspect]
Clear[X]
:[font = input; preserveAspect]
 X[x_,n_Integer] := Sin[n Pi x/L]

:[font = text; inactive; preserveAspect; fontSize = 14]

Since the second derivative of this is - (n Pi/L)^2 Sin[n Pi x/L], we see that
:[font = input; preserveAspect]
mu[n_Integer] := -(n Pi/L)^2
:[font = text; inactive; preserveAspect; fontSize = 14]
Now let's look at the other part of the PDE which we "separated off":
T'[t] = - mu T[t]

As we saw, there are many possible values of \mu which are consistent with the boundary condition for the function X, and the same possibilities occur in equation (6.3) (multiply \mu by k here if k is not set to 1). Therefore
         T[t_,n_] = A[n] Exp[- (n \pi)^2 t].
         
We cannot go further with the function T without being given the initial condition, that is, information about what happens at t=0.

The general solution we come up with is a linear combination of the particular solutions we get by separating the equation. It is not yet obvious that this is a completely general solution, but it is. The most general linear combination is of the form:
:[font = input; preserveAspect; startGroup]
u[t_,x_] = Sum[c[n] Exp[- (n Pi) t ] Sin[n Pi x], \
       {n,1,infinity}] 
:[font = message; inactive; preserveAspect]
General::spell1: 
   Possible spelling error: new symbol name "infinity"
     is similar to existing symbol "Infinity".
:[font = output; output; inactive; preserveAspect; endGroup]
Sum[(c[n]*Sin[n*Pi*x])/E^(n*Pi*t), {n, 1, infinity}]
;[o]
    c[n] Sin[n Pi x]
Sum[----------------, {n, 1, infinity}]
         n Pi t
        E
:[font = text; inactive; preserveAspect; fontSize = 14; endGroup]
The coefficients A[n] have to be determined from the initial condition, u(0,x); they are the Fourier sine coefficients of u(0,x).  (See the examples in the next section.)
:[font = section; inactive; Cclosed; preserveAspect; startGroup]
Automating the solution to the heat equation
The general solution of the heat equation with 
						(i)  homogeneous Dirichlet boundary conditions at x=a and x=b, 
and 
      (ii)  initial condition
      
	           u[0,x] == f[x]
can be automated with the Mathematica commands in this section.
;[s]
4:0,0;45,1;258,2;269,1;296,-1;
3:1,19,14,Times,1,18,0,0,0;2,16,12,Times,0,14,0,0,0;1,16,12,Times,2,14,0,0,0;
:[font = text; inactive; preserveAspect; fontSize = 14]
We break the solution into two steps, the first of which calculates the Fourier sine coefficients for the initial condition.
:[font = input; preserveAspect]
B[f_,m_, {a_,b_}] := (2/(b-a)) Integrate[                  \
         (Sin[m Pi intvar/(b-a)] f  /. x -> intvar),       \
         {intvar, a, b}]            /. TrigId
:[font = input; preserveAspect]
HeatSol[f_,k_,t_,{x_,a_,b_}, infinity_] :=                     \
         Sum[B[f,n, {a,b}] Exp[-(n Pi) k t] Sin[n Pi x/(b-a)], \
         {n,1,infinity}]
:[font = text; inactive; preserveAspect; fontSize = 14]
These commands allow a lot of flexibility.  It will often be more convenient to simplify them a standard form, assuming that the variables have been scaled so that k == 1 and 0 < x < 1:
:[font = input; preserveAspect]
B[f_,m_] := 2 Integrate[                                      \
            (Sin[m Pi intvar] f  /. x -> intvar),             \
            {intvar, 0,1}]            /. TrigId
:[font = input; preserveAspect; endGroup]
HeatSol[f_,t_,x_, infinity_] :=                               \
            Sum[B[f,n] Exp[- (n Pi) t ] Sin[n Pi x],          \
            {n,1,infinity}]
:[font = section; inactive; preserveAspect]
Examples - heat equation
:[font = subsection; inactive; Cclosed; preserveAspect; startGroup]
1.  k = 1, 0 < x < 1, and the rod is initially at a constant temperature.  (Model Problem VIII.1)
:[font = text; inactive; preserveAspect; fontSize = 14]
Model Problem VIII.1.  Find the future temperature within a rod of length 1, k=1, insulated along its length, and with ends in thermal contact with ice water, if at time 0 the rod is at a uniform temperature u(0,x) = 50.
 
;[s]
2:0,1;20,0;223,-1;
2:1,16,12,Times,0,14,0,0,0;1,16,12,Times,1,14,0,0,0;
:[font = input; preserveAspect; startGroup]
HeatSol[50,t,x, 7]
;[s]
2:0,1;7,0;19,-1;
2:1,12,10,Courier,1,12,0,0,0;1,12,10,Courier,1,12,65535,0,0;
:[font = output; output; inactive; preserveAspect]
(200*Sin[Pi*x])/(E^(Pi*t)*Pi) + (200*Sin[3*Pi*x])/(3*E^(3*Pi*t)*Pi) + 
 
  (40*Sin[5*Pi*x])/(E^(5*Pi*t)*Pi) + 
 
  (200*Sin[7*Pi*x])/(7*E^(7*Pi*t)*Pi)
;[o]
200 Sin[Pi x]   200 Sin[3 Pi x]   40 Sin[5 Pi x]   200 Sin[7 Pi x]
------------- + --------------- + -------------- + ---------------
   Pi t             3 Pi t           5 Pi t            7 Pi t
  E     Pi       3 E       Pi       E       Pi      7 E       Pi
:[font = text; inactive; preserveAspect; fontSize = 14]
(We approximate infinity by the large number 7!)  Let us plot the solution in time steps of 0.2:
:[font = input; preserveAspect; endGroup; endGroup]
Table[Plot[%,{x,0,1}, PlotRange -> {0,60}], {t, 0, 1, .2}]
;[s]
3:0,0;11,1;12,0;59,-1;
2:2,12,10,Courier,1,12,0,0,0;1,12,10,Courier,1,12,65535,0,0;
:[font = subsection; inactive; Cclosed; preserveAspect; startGroup]
2.  k = 1, 0 < x < 1, and the rod is initially at temperature 100 x. 
:[font = text; inactive; preserveAspect; fontSize = 14]
Model Problem VIII.1.  Find the future temperature within a rod of length 1, k=1, insulated along its length, and with ends in thermal contact with ice water, if at time 0 the rod is at a uniform temperature u(0,x) = 50.
 
;[s]
2:0,1;20,0;223,-1;
2:1,16,12,Times,0,14,0,0,0;1,16,12,Times,1,14,0,0,0;
:[font = input; preserveAspect; startGroup]
HeatSol[100 x,t,x, 7]
;[s]
2:0,1;7,0;22,-1;
2:1,12,10,Courier,1,12,0,0,0;1,12,10,Courier,1,12,65535,0,0;
:[font = output; output; inactive; preserveAspect]
(200*Sin[Pi*x])/(E^(Pi*t)*Pi) - (100*Sin[2*Pi*x])/(E^(2*Pi*t)*Pi) + 
 
  (200*Sin[3*Pi*x])/(3*E^(3*Pi*t)*Pi) - 
 
  (50*Sin[4*Pi*x])/(E^(4*Pi*t)*Pi) + 
 
  (40*Sin[5*Pi*x])/(E^(5*Pi*t)*Pi) - 
 
  (100*Sin[6*Pi*x])/(3*E^(6*Pi*t)*Pi) + 
 
  (200*Sin[7*Pi*x])/(7*E^(7*Pi*t)*Pi)
;[o]
200 Sin[Pi x]   100 Sin[2 Pi x]   200 Sin[3 Pi x]   50 Sin[4 Pi x]
------------- - --------------- + --------------- - -------------- + 
   Pi t            2 Pi t             3 Pi t           4 Pi t
  E     Pi        E       Pi       3 E       Pi       E       Pi
 
  40 Sin[5 Pi x]   100 Sin[6 Pi x]   200 Sin[7 Pi x]
  -------------- - --------------- + ---------------
     5 Pi t            6 Pi t            7 Pi t
    E       Pi      3 E       Pi      7 E       Pi
:[font = text; inactive; preserveAspect; fontSize = 14]
(We approximate infinity by the large number 7!)  Let us plot the solution in time steps of 0.2:
:[font = input; preserveAspect; endGroup; endGroup]
Table[Plot[%,{x,0,1}, PlotRange -> {0,110}], {t, 0, 1, .2}]
;[s]
3:0,0;11,1;12,0;60,-1;
2:2,12,10,Courier,1,12,0,0,0;1,12,10,Courier,1,12,65535,0,0;
:[font = subsection; inactive; Cclosed; preserveAspect; startGroup]
3.  Observing the smoothing property, approach to equilibrium, and maximum principle.
:[font = text; inactive; preserveAspect; fontSize = 14]
In the graphs of the previous two examples we easily several general features:
       1.  u[t,x] is a smooth function
       2.  u[t,x] tends to the equilibrium value 0
       3.  The way in which u[t,x] approaches equilibrium is the temperature distribution rather rapidly rearranges itself to have the shape of the fundamental normal mode Sin[Pi x].
       4.  The maximum and minimum values of the temperature occur at the space-time boundary:  x = 0, x = 1, t = 0, or t = tmax (whatever the stopping time is).

The initial conditions for these problems were, however, particularly simple.  Here are some other solutions of the heat equation the graphs of which will reveal the same features.  We construct them as linear combinations of simple solutions of the heat equation, with various boundary conditions.

1.
;[s]
3:0,0;477,1;480,0;818,-1;
2:2,16,12,Times,0,14,0,0,0;1,24,16,Times,64,14,0,0,0;
:[font = input; preserveAspect; startGroup]
HeatOp[2t + x^2- 5] == 0
:[font = output; output; inactive; preserveAspect; endGroup]
True
;[o]
True
:[font = input; preserveAspect]
u1[t_,x_] = 2t + x^2- 5 + Sum[(-1)^(n+1) 10 /(1 + Sqrt[n]) Exp[- n^2 Pi^2 t] \
      Sin[n Pi x], {n, 1, 6}];
:[font = input; preserveAspect]
Table[Plot[u1[t,x],{x,0,1}, PlotRange -> {-15,15}], {t, 0, .2, .05}]
;[s]
3:0,0;11,1;13,0;69,-1;
2:2,12,10,Courier,1,12,0,0,0;1,12,10,Courier,1,12,65535,0,0;
:[font = text; inactive; preserveAspect]
2.
;[s]
1:0,1;3,-1;
2:0,13,9,Times,0,12,0,0,0;1,16,12,Times,0,14,0,0,0;
:[font = input; preserveAspect; startGroup]
HeatOp[Exp[- n^2 Pi^2 t] Sin[n Pi x + n]] == 0
:[font = output; output; inactive; preserveAspect; endGroup]
True
;[o]
True
:[font = input; preserveAspect]
u2[t_,x_] = Sum[(-1)^n n Exp[- n^2 Pi^2 t] \
      Sin[n Pi x + n], {n, 1, 6}];
:[font = input; preserveAspect; endGroup]
Table[Plot[u2[t,x],{x,0,1}, PlotRange -> {-15,15}], {t, 0, .1, .02}]
;[s]
3:0,0;11,1;13,0;69,-1;
2:2,12,10,Courier,1,12,0,0,0;1,12,10,Courier,1,12,65535,0,0;
:[font = section; inactive; Cclosed; preserveAspect; startGroup]
The temperature distribution with radiative loss of energy, using method of separation of variables
:[font = text; inactive; dontPreserveAspect; fontSize = 14]
If the insulation along the rod is imperfect, we may model the loss with a modified heat operator of the form:
:[font = input; preserveAspect]
ModHeatOp[u_] := D[u,t] - D[u,{x,2}] - .2 u
:[font = text; inactive; preserveAspect]
We shall impose Dirichlet conditions at 0 and 1.
;[s]
1:0,1;49,-1;
2:0,13,9,Times,0,12,0,0,0;1,16,12,Times,0,14,0,0,0;
:[font = input; preserveAspect; startGroup]
Simplify[ModHeatOp[T[t] X[x]]/(T[t] X[x])]
;[s]
3:0,0;9,1;18,0;43,-1;
2:2,12,10,Courier,1,12,0,0,0;1,12,10,Courier,1,12,65535,0,0;
:[font = output; output; inactive; dontPreserveAspect; endGroup]
-0.2 + Derivative[1][T][t]/T[t] - Derivative[2][X][x]/X[x]
;[o]
       T'[t]   X''[x]
-0.2 + ----- - ------
       T[t]     X[x]
:[font = text; inactive; preserveAspect; fontSize = 14]
The eigenvalue equation becomes:
:[font = input; preserveAspect; startGroup]
 DSolve[- X''[x] + (1/5) X[x] == mu X[x],X[x],x]
:[font = output; output; inactive; dontPreserveAspect; endGroup]
{{X[x] -> C[1]/E^(((1 - 5*mu)^(1/2)*x)/5^(1/2)) + 
 
     E^(((1 - 5*mu)^(1/2)*x)/5^(1/2))*C[2]}}
;[o]
                     C[1]
{{X[x] -> --------------------------- + 
           (Sqrt[1 - 5 mu] x)/Sqrt[5]
          E
 
      (Sqrt[1 - 5 mu] x)/Sqrt[5]
     E                           C[2]}}
:[font = text; inactive; preserveAspect; fontSize = 14]
In order that x[0] == 0, C[2] = - C[1], so the eigenfunctions are proportional to
:[font = input; preserveAspect]
X[x_,mu_] := Sin[Sqrt[5 mu - 1] x/Sqrt[5]]
:[font = text; inactive; preserveAspect; fontSize = 14]
In order for this to equal 0 when x=1, mu -> (5 n^2 Pi^2 - 1)/5, so the separated solutions become 
        Exp[- ((1 - 5 n^2 Pi^2)/5) t] Sin[n Pi x]
:[font = input; preserveAspect]
ModHeatSol[f_,t_,x_, infinity_] :=                               \
            Sum[B[f,n] Exp[- ((5 n^2 Pi^2 - 1)/5) t ] Sin[n Pi x],          \
            {n,1,infinity}]
:[font = input; preserveAspect]
ModHeatSol[4 x^2,t,x,8];
;[s]
2:0,1;10,0;25,-1;
2:1,12,10,Courier,1,12,0,0,0;1,12,10,Courier,1,12,65535,0,0;
:[font = input; preserveAspect; endGroup]
Table[Plot[%,{x,0,1}, PlotRange -> {0,5}], {t, 0, 1, .2}]
;[s]
3:0,0;11,1;12,0;58,-1;
2:2,12,10,Courier,1,12,0,0,0;1,12,10,Courier,1,12,65535,0,0;
:[font = section; inactive; Cclosed; preserveAspect; startGroup]
Exercises
:[font = subsection; inactive; preserveAspect]
1.  Use Mathematica to automate the solution of the heat equation with Neumann boundary conditions at a and b.  Use the result to solve some specificinitial-boundary-value problems of your own design.
;[s]
3:0,0;8,1;19,0;201,-1;
2:2,16,12,Times,1,14,0,0,0;1,16,12,Times,3,14,0,0,0;
:[font = subsection; inactive; preserveAspect]
2.  Modify the solutions in this notebook to solve an initial-boundary-value problem with nonhomogeneous boundary conditions, for example:
	              ut  == uxx
	              ux[t,0] = 0, ux[t,1] = 1
	              u[0,x] = 0. 
;[s]
9:0,0;155,1;156,0;162,1;164,0;181,1;182,0;194,1;195,0;233,-1;
2:5,16,12,Times,1,14,0,0,0;4,24,16,Times,65,14,0,0,0;
:[font = subsection; inactive; preserveAspect; endGroup; endGroup]
3.  We may model a rod in ice water which is in thermal contact at its ends, imperfectly insulated along its length, and with an internal source of heat (chemical or radioactive) by:
	              ut  == uxx= - (1/5) u + 1
	              ux[t,0] = 0, u[t,1] = 0
	              u[0,x] = 0. 
	Solve.
;[s]
7:0,0;199,1;200,0;206,1;208,0;240,1;241,0;299,-1;
2:4,16,12,Times,1,14,0,0,0;3,24,16,Times,65,14,0,0,0;
^*)