(*^
::[	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; 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; Cclosed; 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.
;[s]
13: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;721,-1;
3:7,13,9,Times,0,12,0,0,0;5,13,9,Times,2,12,0,0,0;1,13,9,Times,0,12,65535,0,0;
:[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.
:[font = section; inactive; Cclosed; preserveAspect; startGroup]
Introduction.  The heat equation and the heat operator.
:[font = text; inactive; dontPreserveAspect; fontSize = 14]
In Mathematica  notation, the wave equation can be written as
					
					D[u, T] = k D[u, {x,2}] 
					
We classify this PDE as second order  and linear.  
;[s]
5:0,0;3,1;14,0;72,2;96,0;156,-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 = text; inactive; preserveAspect; fontSize = 14]
"Second order" means that derivatives up to second order occur in the equation, but no higher-order derivatives.  This will be the case for the great majority of PDE's in this class.  We shall solve the heat equation with a the method known as separation of variables, which is often particularly good for linear equations.  "Linear" refers to the superposition property, i.e., if two functions u1 and u2 are solutions, then so is any linear combination, a1 u1 + a2 u2. 

In fact, (HE) is a statement that the function u is in the kernel of a linear operator. 

A convenient way to set up a linear differential operator in Mathematica is as follows - the operator acts on functions, so identify a function f with an understroke to show that it is the independent variable in the 
wave operator:
;[s]
15:0,0;396,1;397,0;403,1;404,0;455,2;456,1;457,0;459,1;460,0;463,2;464,1;465,0;467,1;468,0;795,-1;
3:7,16,12,Times,0,14,0,0,0;6,24,16,Times,64,14,0,0,0;2,16,12,Symbol,0,14,0,0,0;
:[font = input; preserveAspect]
HeatOp[f_,k_] := (1/k) D[f,t] -  D[f, {x,2}]
:[font = text; inactive; preserveAspect; fontSize = 14]
We'll often set k = 1 for convenience; this is not really a restriction, since we are free to choose a clock with scaled time variable t' = k t.  To make this simplification convenient, we define:
:[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;
^*)