{VERSION 5 0 "IBM INTEL NT" "5.0" }
{USTYLETAB {CSTYLE "Maple Input" -1 0 "Courier" 0 1 255 0 0 1 0 1 0 0 
1 0 0 0 0 1 }{CSTYLE "2D Math" -1 2 "Times" 0 1 0 0 0 0 0 0 2 0 0 0 0 
0 0 1 }{CSTYLE "2D Input" 2 19 "" 0 1 255 0 0 1 0 0 0 0 0 0 0 0 0 1 }
{CSTYLE "" -1 256 "" 1 16 0 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 
257 "" 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 }{PSTYLE "Normal" -1 0 1 
{CSTYLE "" -1 -1 "Times" 1 12 0 0 0 1 2 2 2 2 2 2 1 1 1 1 }1 1 0 0 0 
0 1 0 1 0 2 2 0 1 }}
{SECT 0 {EXCHG {PARA 0 "" 0 "" {TEXT 256 13 "Heat Equation" }}{PARA 0 
"" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 187 "A derivation is \+
given in the notes (jpg-files in the folder). We follow David Betounes
: Partial Differential Equations for Computational Science, Springer V
erlag 1997, ISBN 0-387-98300-7." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}
{PARA 0 "" 0 "" {TEXT -1 176 "We begin with a one-dimensional example.
(Example 2.2 in the book). A thin rod lies on the x-axis from 0 to 1 (
appropriate units are chosen to make the equation dimensionless). " }}
{PARA 0 "" 0 "" {TEXT -1 65 "At time zero the rod is at zero temperatu
re (inital condition).  " }{TEXT 19 10 "u(x,0) = 0" }{TEXT -1 1 "." }}
{PARA 0 "" 0 "" {TEXT -1 158 "On the left boundary the rod is kept at \+
zero temperature by a heat bath (make it zero Celsius to avoid trouble
 with physics, zero Kelvin are not reachable!). " }}{PARA 0 "" 0 "" 
{TEXT 19 10 "u(0,t) = 0" }{TEXT -1 6 "  for " }{TEXT 19 5 "t > 0" }}
{PARA 0 "" 0 "" {TEXT -1 101 "On the right boundary the rod is insulat
ed, i.e., no flow of heat across this boundary is permitted. " }{TEXT 
19 28 "diff(u(x,t),x)[@x=1,t>0] = 0" }{TEXT -1 28 " (this is not Maple
 syntax!)" }}{PARA 0 "" 0 "" {TEXT -1 172 "The heat equation contains \+
a source term, i.e., it is not homogeneous. The inhomogeneity is chose
n as a constant, i.e., heat is being generated in all locations of the
 rod " }{TEXT 19 9 "0 < x < 1" }{TEXT -1 9 " (not at " }{TEXT 19 5 "x \+
= 0" }{TEXT -1 2 ")." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "
" {TEXT -1 96 "The heat equation has a steady-state solution, i.e., a \+
solution that is being reached for large " }{TEXT 19 1 "t" }{TEXT -1 
38 ". For this to be possible, a solution " }{TEXT 19 6 "u(x,t)" }
{TEXT -1 127 " has to be found with vanishing time derivative in the l
imit of large time. This function will depend on the position variable
 " }{TEXT 19 1 "x" }{TEXT -1 6 " only." }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 8 "restart;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 43 
"HEq:=diff(u(x,t),t) = diff(u(x,t),x$2) + 1;" }}}{EXCHG {PARA 0 "" 0 "
" {TEXT -1 46 "In order to express the boundary condition at " }{TEXT 
19 3 "x=1" }{TEXT -1 31 " we have to understand how the " }{TEXT 19 1 
"D" }{TEXT -1 50 " operator works on functions of several variables:" 
}}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 28 "BC:=u(0,t)=0,D[1](u)(1,t)
=0;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 13 "IC:=u(x,0)=0;" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 33 "sol:=pdsolve(\{HEq,BC,IC\},u
(x,t));" }{TEXT -1 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 138 "It see
ms as if Maple has not found a unique solution. Four constants appear \+
in the answer. Is it possible that many solutions truly exist?" }}
{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 141 "Let us e
xtract the steady state solution by solving the ordinary differential \+
equation that follows from setting the time derivative to zero." }}
{PARA 0 "" 0 "" {TEXT -1 38 "Let us call the steady-state solution " }
{TEXT 19 4 "w(x)" }{TEXT -1 1 ":" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 23 "DE:=diff(w(x),x$2)+1=0;" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 24 "BCode:=w(0)=0,D(w)(1)=0;" }}}{EXCHG {PARA 0 "> " 0 "
" {MPLTEXT 1 0 26 "solDE:=dsolve(\{DE,BCode\});" }}}{EXCHG {PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 24 "plot(rhs(solDE),x=0..1);" }}}{EXCHG {PARA 0 "
" 0 "" {TEXT -1 61 "The heat equation is supposed to describe the evol
ution from " }{TEXT 19 8 "u(x,0)=0" }{TEXT -1 9 " towards " }{TEXT 19 
18 "u(x,infinity)=w(x)" }{TEXT -1 1 "." }}{PARA 0 "" 0 "" {TEXT -1 
177 "Note how the boundary conditions are satisified (prescribed tempe
rature value on the left, and zero derivative on the right). These con
ditions have to be satisfied at all times." }}{PARA 0 "" 0 "" {TEXT 
-1 201 "The equation is usually solved by the separation of variables \+
method. Before we do that let us look at the solution generated by map
le. Without the boundary conditions the following answer is obtained:
" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 18 "sol:=pdsolve(HEq);" }}}
{EXCHG {PARA 0 "" 0 "" {TEXT -1 76 "Now it becomes clearer what the  c
onstants mean: pdsolve is telling us that:" }}{PARA 0 "" 0 "" {TEXT 
-1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 152 "1) there is a time-independent
, additive part to the solution, which is characterized by three const
ants. From our asymptotic solution it is clear that " }}}{EXCHG {PARA 
0 "> " 0 "" {MPLTEXT 1 0 47 "solHEq:=subs(_C1=1,_C2=1,_C3=0,rhs(op(1,s
ol)));" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 128 "The heat equation admi
ts a solution of the form where the transient part is expressed as a p
roduct of one-dimensional functions " }{TEXT 19 11 "F1(x)*F2(t)" }
{TEXT -1 77 ". The two functions satisfy simple ordinary differential \+
equations listed in " }{TEXT 19 3 "sol" }{TEXT -1 1 "." }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 28 "eq1:=diff(F2(t),t)=c1*F2(t);" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 30 "eq2:=diff(F1(x),x$2)=c1*F1(x
);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 46 "sol_x:=dsolve(\{eq2,F
1(0)=0,D(F1)(1)=0\},F1(x));" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 68 "We
 are not interested in the trivial solution, so what do we do now?" }}
}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 35 "sol_x:=dsolve(\{eq2,F1(0)=0
\},F1(x));" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 131 "Somehow we will ha
ve to find acceptable solutions (derivative to vanish on the right-han
d boundary) by choosing c1 rather than _C2." }}}{EXCHG {PARA 0 "> " 0 
"" {MPLTEXT 1 0 31 "sol_x1:=subs(_C2=1,rhs(sol_x));" }}}{EXCHG {PARA 
0 "> " 0 "" {MPLTEXT 1 0 40 "plot(subs(x=1,diff(sol_x1,x)),c1=0..50);
" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 35 "This shows that positive valu
es of " }{TEXT 19 2 "c1" }{TEXT -1 88 " won't do the trick. So we try \+
negative values, which makes the solution complex-valued." }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 72 "[Re(subs(x=1,c1=-5,diff(sol_x1,x)))
,Im(subs(x=1,c1=-5,diff(sol_x1,x)))];" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 119 "plot([Re(subs(x=1,diff(sol_x1,x))),Im(subs(x=1,diff(
sol_x1,x)))],c1=-220..0,color=[red,blue],axes=boxed,numpoints=500);" }
}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 185 "What are we learning here? We t
ook the condition for the right-hand boundary, and we showed that it c
an be satisfied by the solution to the spatial ODE provided the separa
tion constant " }{TEXT 19 2 "c1" }{TEXT -1 132 " is chosen to be negat
ive. The real part of the condition is always zero (boring, but import
ant). The derivative of the solution at " }{TEXT 19 3 "x=1" }{TEXT -1 
105 " does vanish when the imaginary part of the expression crosses ze
ro. We can use the integration constant " }{TEXT 19 3 "_C2" }{TEXT -1 
59 " to make the solution real-valued (by making it imaginary)." }}
{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 117 "Therefor
e, the product part of the solution, i.e., the transient part, is not \+
unique. There are different choices of " }{TEXT 19 2 "c1" }{TEXT -1 
152 " that work, and we have to take them all into account in order to
 find a unique solution that matches our boundary conditions and the i
nitial condition." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" 
{TEXT -1 32 "Each spatial solution (function " }{TEXT 19 5 "F1(x)" }
{TEXT -1 86 " that satisfies the ODE and the boundary conditions) has \+
an appropriate time solution " }{TEXT 19 5 "F2(t)" }{TEXT -1 112 " tha
t goes with it. Before we find it, let us look at the spatial modes. F
irst we determine a few of the lowest " }{TEXT 19 2 "c1" }{TEXT -1 19 
" values accurately." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 82 "sol
_c1:='sol_c1': sol_c1[1]:=fsolve(Im(subs(x=1,diff(sol_x1,x))),c1=-10..
-0.0001);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 60 "sol_c1[2]:=fso
lve(Im(subs(x=1,diff(sol_x1,x))),c1=-50..-10);" }}}{EXCHG {PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 61 "sol_c1[3]:=fsolve(Im(subs(x=1,diff(sol_x1,x)))
,c1=-100..-50);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 62 "sol_c1[4
]:=fsolve(Im(subs(x=1,diff(sol_x1,x))),c1=-150..-100);" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 62 "sol_c1[5]:=fsolve(Im(subs(x=1,diff(
sol_x1,x))),c1=-250..-150);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 203 "O
f course, there are infinitely many solutions, and it is not unreasona
ble to ask at this point whether we will get anywhere at all... Nevert
heless, we'll proceed by looking at some of the spatial modes:" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 90 "plot([seq(subs(c1=sol_c1[i],
I*sol_x1),i=1..5)],x=0..1,color=[red,blue,green,black,brown]);" }}}
{EXCHG {PARA 0 "" 0 "" {TEXT -1 67 "These functions look familiar. So \+
let us understand why that is so." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}
{PARA 0 "" 0 "" {TEXT -1 233 "The differential equation for the spatia
l modes is asking: which functions are proportional to their second de
rivative. This selects the class of exponentials, including complex ex
ponentials. Then we asked the functions to vanish at " }{TEXT 19 3 "x=
0" }{TEXT -1 165 ", and we wanted real-valued functions. This selected
 from sines and cosines just sine functions. Next we asked only for th
ose functions which had zero derivative at " }{TEXT 19 3 "x=1" }{TEXT 
-1 25 ". This selected from all " }{TEXT 19 15 "sin(sqrt(c1)*x)" }
{TEXT -1 54 " those with an appropriate wavenumber, i.e., value of " }
{TEXT 19 2 "c1" }{TEXT -1 91 ". Of those we plotted the lowest four. F
unctions which do not have too much variation with " }{TEXT 19 1 "x" }
{TEXT -1 216 " should be representable with a few low-lying modes to a
 reasonable degree of accuracy. This is where our hope comes from: we \+
proceed even though we know that in principle there are infinitely man
y product solutions." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "
" {TEXT -1 51 "Let us now understand how the modes evolve in time." }}
{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 
0 27 "sol_t:=dsolve(\{eq1\},F2(t));" }}}{EXCHG {PARA 0 "" 0 "" {TEXT 
-1 43 "We observe that we want negative values of " }{TEXT 19 2 "c1" }
{TEXT -1 185 ", because this generates dying solutions. Each of the sp
atial modes dies with time with a time constant of its own! The higher
 modes (many oscillations) die faster than the lower modes." }}{PARA 
0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 113 "Now one might \+
ask, why do we need to worry about dying modes, when the solution to t
he heat equation starts from " }{TEXT 19 8 "u(x,0)=0" }{TEXT -1 286 " \+
? The answer is that we have the additive asymptotic piece. To be cons
istent with it, and with the stated initial condition we have to super
impose the modes in such a way that they cancel the asymptotic piece. \+
This we do by a projection method (formally we construct a Fourier ser
ies)." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 6 "solDE;" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 54 "basis:=seq(evalc(subs(c1=sol_c1[i],
I*sol_x1)),i=1..5);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 130 "We claim \+
that the inner product (dot product) of two basis functions is just th
e integral of the product of functions from 0 to 1." }}}{EXCHG {PARA 
0 "> " 0 "" {MPLTEXT 1 0 30 "int(basis[1]*basis[2],x=0..1);" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 30 "int(basis[1]*basis[3],x=0..1
);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 81 "To numerical accuracy the f
unctions are orthogonal! What about the normalization?" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 30 "int(basis[1]*basis[1],x=0..1);" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 30 "int(basis[2]*basis[2],x=0..1
);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 353 "OK, so the functions are n
ormalized in a funny way, but that shouldn't deter us. Verify the orth
ogonality and normalization for some of the other basis functions (or \+
basis vectors in the Hilbert space). Let us expand by calculating proj
ections. To make up for the wrong normalization we re-define the basis
 functions by dividing by the square root of 2." }}}{EXCHG {PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 46 "b[1]:=int(rhs(solDE)*basis[1]/sqrt(2),x=0..1);
" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 46 "b[2]:=int(rhs(solDE)*ba
sis[2]/sqrt(2),x=0..1);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 46 "
b[3]:=int(rhs(solDE)*basis[3]/sqrt(2),x=0..1);" }}}{EXCHG {PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 46 "b[4]:=int(rhs(solDE)*basis[4]/sqrt(2),x=0..1);
" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 46 "b[5]:=int(rhs(solDE)*ba
sis[5]/sqrt(2),x=0..1);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 168 "Good \+
news: the projection of the function to be expanded onto higher basis \+
states are getting smaller. We are ready to compare the expansion with
 the original function:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 38 "
FE:=add(b[i]*basis[i]/sqrt(2),i=1..5);" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 46 "plot([FE,rhs(solDE)],x=0..1,color=[red,blue]);" }}}
{EXCHG {PARA 0 "" 0 "" {TEXT -1 202 "We see that the expansion is bett
er than 1% accurate (otherwise we would notice a difference in the gra
ph). Repeat the graph with less than five basis states and observe how
 the convergence comes about!" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}
{PARA 0 "" 0 "" {TEXT -1 218 "Now we are ready to put everything toget
her. The total solution is obtained by multiplying the basis functions
 with the appropriate time-factors, and by adding the transient and st
eady-state parts together as required:" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 6 "sol_t;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 67 "F
Et:=rhs(solDE)-add(b[i]*basis[i]/sqrt(2)*exp(sol_c1[i]*t),i=1..5);" }}
}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 81 "plot([seq(subs(t=j*0.25,FEt
),j=1..5)],x=0..1,color=[red,blue,green,black,brown]);" }}}{EXCHG 
{PARA 0 "" 0 "" {TEXT -1 154 "We observe how the solution rapidly appr
oaches the steady state. Explore the time behavior in detail (by choos
ing different time intervals and more steps " }{TEXT 19 1 "j" }{TEXT 
-1 2 ")." }}{PARA 0 "" 0 "" {TEXT -1 120 "Repeat the observation for a
 more approximate solution by selecting a smaller expansion in the Fou
rier series, i.e., in " }{TEXT 19 3 "FEt" }{TEXT -1 1 "." }}{PARA 0 "
" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 165 "Here is a graph i
n three dimensions. It shows, how well the Fourier expansion approxima
tes the initial condition, and how quickly the asymptotic values are a
ttained." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 37 "plot3d(FEt,x=0.
.1,t=0..2,axes=boxed);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 69 "Rotate \+
the graph to look at the solution from different perspectives." }}
{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 106 "To remin
d ourselves that we do not have an exact solution we show the deviatio
n from the intial condition:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 
0 27 "plot(subs(t=0,FEt),x=0..1);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT 
257 9 "Exercise:" }}{PARA 0 "" 0 "" {TEXT -1 103 "Increase the expansi
on to 10 basis functions and observe how well the initial condition ca
n be matched." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}
{MARK "66" 0 }{VIEWOPTS 1 1 0 1 1 1803 1 1 1 1 }{PAGENUMBERS 0 1 2 33 
1 1 }