{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 }{CSTYLE "" -1 258 "" 0 1 0 0 
0 0 0 0 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 259 "" 0 1 0 0 0 0 0 0 1 0 0 0 
0 0 0 0 }{CSTYLE "" -1 260 "" 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 }
{CSTYLE "" -1 261 "" 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 
262 "" 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 263 "" 0 1 0 0 
0 0 0 0 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 264 "" 0 1 0 0 0 0 0 0 1 0 0 0 
0 0 0 0 }{CSTYLE "" -1 265 "" 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 14 "Fourier series" }}{PARA 
0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 17 "In the workshee
t " }{TEXT 19 11 "HeatEqn.mws" }{TEXT -1 553 " we have seen how a part
ial differential equation (the heat equation) in one space and one tim
e dimension admitted solutions in product form (separation of variable
s technique), and how the solution was approximated by expanding the s
patial part in a Fourier sine series. In this worksheet we generalize \+
the idea of representing functions defined on a finite interval (part \+
of the real axis) by means of an expansion in basis functions. The con
cepts can then be extended to basis functions defined as the eigenfunc
tions of some Sturm-Liouville operator." }}{PARA 0 "" 0 "" {TEXT -1 0 
"" }}{PARA 0 "" 0 "" {TEXT -1 296 "We begin by defining our basic inte
rval on which the function is to be defined. This would normally be th
e interval on which solutions to a differential equation are seeked, i
.e., usually the unknown function (solution to the ODE) has specified \+
values at the edges of the interval. The choice of " }{TEXT 19 1 "a" }
{TEXT -1 5 " and " }{TEXT 19 1 "b" }{TEXT -1 155 " below is arbitrary,
 and depends on the problem posed. In quantum mechanics we could be as
king for particles to be located in a box bounded by these edges." }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 21 "restart; a:=-2; b:=3;" }}}
{EXCHG {PARA 0 "" 0 "" {TEXT -1 128 "We define an inner product for re
al-valued functions (without weight function) assuming that the indepe
ndent variable is called " }{TEXT 19 1 "x" }{TEXT -1 1 ":" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 35 "IP:=(f,g)->int(expand(f*g),x=a..b);
" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 132 "This inner (dot) product wil
l be used to calculate the norm of functions, and the projection of fu
nctions onto basis functions. The " }{TEXT 19 6 "expand" }{TEXT -1 
104 " function has been included to help Maple calculate integrals by \+
breaking them into many simpler pieces." }}{PARA 0 "" 0 "" {TEXT -1 0 
"" }}{PARA 0 "" 0 "" {TEXT -1 94 "The basis will be defined as the eig
enfunctions of the following simple differential operator:" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "DO:=f->-diff(f,x$2);" }}}{EXCHG 
{PARA 0 "" 0 "" {TEXT -1 209 "In order to determine the eigenfunctions
 and eigenvalues of this differential operator on the interval [a,b] w
e have to specify some boundary conditions (BC). We anticipate (based \+
on elementary Calculus) that " }{TEXT 19 8 "sin(k*x)" }{TEXT -1 5 " an
d " }{TEXT 19 8 "cos(k*x)" }{TEXT -1 49 " are candidates for the eigen
functions, and that " }{TEXT 19 3 "k^2" }{TEXT -1 83 " will appear as \+
an eigenvalue. The boundary conditions will impose which values of " }
{TEXT 19 1 "k" }{TEXT -1 15 " are permitted." }}{PARA 0 "" 0 "" {TEXT 
-1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 30 "In the example carried out in \+
" }{TEXT 19 11 "HeatEqn.mws" }{TEXT -1 108 " we encountered two types \+
of BCs: we can fix the value of the basis function, or we can fix its \+
derivative. " }}{PARA 0 "" 0 "" {TEXT -1 22 "The eigenvalue problem" }
}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT 19 34 "DO(f) = \+
lambda*f     ,  lambda=k^2" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 
"" 0 "" {TEXT -1 37 "is homogeneous, which means that if  " }{TEXT 19 
4 "f(x)" }{TEXT -1 68 "  is a solution, then any multiple of it will b
e a solution as well." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 
"" {TEXT -1 59 "If we do not impose boundary conditions, then the valu
e of " }{TEXT 19 1 "k" }{TEXT -1 248 " will be unrestricted, and the i
nterval [a,b] becomes somehow irrelevant. This means that the function
s are defined on the entire real axis (something which becomes importa
nt in quantum mechanics when discussing plane-wave free-particle solut
ions)." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 86 
"Let us start with the simple case of imposing that the basis function
s vanish at both " }{TEXT 19 1 "a" }{TEXT -1 5 " and " }{TEXT 19 1 "b
" }{TEXT -1 35 ". These BCs are called homogeneous." }}}{EXCHG {PARA 
0 "> " 0 "" {MPLTEXT 1 0 12 "BC1:=f(a)=0;" }}}{EXCHG {PARA 0 "> " 0 "
" {MPLTEXT 1 0 22 "SL:=DO(f(x))=k^2*f(x);" }}}{EXCHG {PARA 0 "" 0 "" 
{TEXT -1 160 "How do we find the basis functions that satisfy this equ
ation? We can't use the second boundary condition, because that would \+
force the boring trivial solution " }{TEXT 19 6 "f(x)=0" }{TEXT -1 42 
". We have to use the choice of eigenvalue " }{TEXT 19 6 "lambda" }
{TEXT -1 31 " to find non-trivial solutions." }}}{EXCHG {PARA 0 "> " 
0 "" {MPLTEXT 1 0 27 "sol:=dsolve(\{SL,BC1\},f(x));" }}}{EXCHG {PARA 
0 "" 0 "" {TEXT -1 35 "The choice of integration constant " }{TEXT 19 
3 "_C1" }{TEXT -1 105 " is arbitrary and defines the normalization of \+
the solution. We impose now the second boundary condition:" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "fB:=subs(_C1=1,rhs(sol));" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "BC2:=subs(x=b,fB)=0;" }}}
{EXCHG {PARA 0 "" 0 "" {TEXT -1 198 "This is a nonlinear transcendenta
l equation which admits infinitely many roots. To find the roots we ha
ve to employ a numerical solver. First, let us graph the expression wh
ose roots we are seeking." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 
55 "plot(lhs(BC2),k=-5..5,view=[-5..5,-4..4],discont=true);" }}}
{EXCHG {PARA 0 "" 0 "" {TEXT -1 29 "Clearly, we need to consider " }
{TEXT 19 6 "k >= 0" }{TEXT -1 48 " only. The functions are labelled by
 eigenvalue " }{TEXT 19 3 "k^2" }{TEXT -1 57 ", which cannot distingui
sh positive from negative values." }}{PARA 0 "" 0 "" {TEXT -1 168 "The
 poles that appear in-between the roots can make life difficult for a \+
systematic root search, but at least it is clear where they occur, as \+
they are associated with " }{TEXT 19 8 "tan(2*k)" }{TEXT -1 23 ", i.e.
, they occur for " }{TEXT 19 10 "2*k = m*Pi" }{TEXT -1 8 ", where " }
{TEXT 19 9 "m=1,2,..." }}{PARA 0 "" 0 "" {TEXT -1 9 "The root " }
{TEXT 19 5 "k = 0" }{TEXT -1 109 " does not lead to an interesting sol
ution (due to the homogeneous BCs imposed we don't have a constant ter
m)." }}{PARA 0 "" 0 "" {TEXT -1 30 "The expression is periodic in " }
{TEXT 19 1 "k" }{TEXT -1 13 " with period " }{TEXT 19 4 "2*Pi" }{TEXT 
-1 1 "." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 38 "k1:=fsolve(BC2,k
=0.0001..Pi/4-0.0001);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 36 "k
2:=fsolve(BC2,k=Pi/4+0.0001..Pi/2);" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 38 "k3:=fsolve(BC2,k=Pi/2..3*Pi/4-0.0001);" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 23 "k4:=fsolve(BC2,k=2..3);" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 33 "k5:=fsolve(BC2,k=3*Pi/4..5*Pi/4);" 
}}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 36 "k6:=fsolve(BC2,k=Pi+0.000
1..5*Pi/4);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 16 "evalf(k6-Pi)
,k1;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 63 "Maple's solve command fin
ds analytic expressions for the roots:" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 18 "sol:=solve(BC2,k);" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 17 "simplify(sol[3]);" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 11 "evalf(sol);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 154 "
We will be content with using the first five basis states, but the per
iodicity allows one to construct all other eigenvalues easily from the
 five numbers " }{TEXT 19 6 "k1..k5" }{TEXT -1 1 "." }}{PARA 0 "" 0 "
" {TEXT -1 30 "Let us graph the basis states:" }}}{EXCHG {PARA 0 "> " 
0 "" {MPLTEXT 1 0 78 "plot([seq(subs(k=k||j,fB),j=1..5)],x=a..b,color=
[red,blue,green,black,brown]);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 
116 "We recognize that there are sine-type and cosine-type solutions p
resent by looking at the behaviour with respect to " }{TEXT 19 3 "x=0
" }{TEXT -1 1 "." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 23 "Bf:=i->
subs(k=k||i,fB);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 16 "IP(Bf(1
),Bf(1));" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 16 "IP(Bf(2),Bf(2)
);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 37 "The functions are not norma
lized yet." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 16 "IP(Bf(1),Bf(2
));" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 16 "IP(Bf(1),Bf(3));" }}
}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 16 "IP(Bf(2),Bf(3));" }}}
{EXCHG {PARA 0 "" 0 "" {TEXT -1 121 "They are orthogonal within numeri
cal accuracy (numerics come into play due to the approximate nature of
 the eigenvalues)." }}{PARA 0 "" 0 "" {TEXT -1 52 "We can generate a l
ist of normalized eigenfunctions:" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 47 "Bfn:=[seq(Bf(i)/sqrt(IP(Bf(i),Bf(i))),i=1..5)];" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 52 "plot(Bfn,x=a..b,color=[red,b
lue,green,black,brown]);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 138 "Now \+
let us expand some known functions in this basis. We begin with a func
tion that satisfies the boundary conditions obeyed by the basis." }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "f1:=(x-a)^3*(x-b)^2;" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 16 "plot(f1,x=a..b);" }}}{EXCHG 
{PARA 0 "" 0 "" {TEXT -1 54 "The expansion coefficients are obtained b
y projection:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 33 "c_n:=[seq(
IP(f1,Bfn[j]),j=1..5)];" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 33 "
f1_ap:=add(c_n[n]*Bfn[n],n=1..5);" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 41 "plot([f1,f1_ap],x=a..b,color=[red,blue]);" }}}{EXCHG 
{PARA 0 "" 0 "" {TEXT -1 70 "Five basis functions are sufficient to re
ach almost 1%-level accuracy!" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}
{PARA 0 "" 0 "" {TEXT 258 11 "Exercise 1:" }}{PARA 0 "" 0 "" {TEXT -1 
120 "Change the power for the (x-a) and (x-b) factors and repeat the c
alculation with five basis states. What do you observe?" }}{PARA 0 "" 
0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT 257 11 "Exercise 2:" }}
{PARA 0 "" 0 "" {TEXT -1 170 "Increase the order of the expansion by f
inding the next five eigenvalues (you can use the periodicity), and ei
genfunctions. Does the convergence behaviour become evident?" }}{PARA 
0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 513 "This was an ex
ample of applying linear algebra techniques to a function space. We us
ed an eigenvalue problem to define a set of eigenvectors (basis functi
ons), and then we expanded some arbitrary vector (function) in this ba
sis. The economy is obvious: we can have approximate knowledge of a fu
nction by knowing a set of numbers (the expansion coefficients). Often
 one can adapt a basis to a problem at hand, and, thus, limit the numb
er of expansion coefficients required for a reasonably accurate repres
entation." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 
72 "We demonstrate now what happens when we look outside the interval \+
[a,b]:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 83 "plot([f1,f1_ap],x
=-10..10,color=[red,blue],view=[-10..10,-200..200],numpoints=500);" }}
}{EXCHG {PARA 0 "" 0 "" {TEXT -1 124 "We observe that the answer obtai
ned from the Fourier series is periodic, i.e., it repeats itself outsi
de the interval [a,b]." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 
0 "" {TEXT 259 11 "Exercise 3:" }}{PARA 0 "" 0 "" {TEXT -1 177 "Change
 the interval [a,b] to cover a wider range (e.g., [-3,4]), repeat the \+
calculation for the Fourier series, and make your observation. Investi
gate the convergence behaviour." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}
{PARA 0 "" 0 "" {TEXT -1 395 "For functions that do not vanish on the \+
edge of the interval some interesting pheonomena can be observed. The \+
textbooks on the subject develop the calculation of Fourier coefficien
ts for simple periodic functions (square wave train, sawtooth function
, etc.), for which a formula for the expansion coefficients can be der
ived. This permits to study convergence behaviour in a straightforward
 way." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 90 "
Let us define the next wavenumbers for the Fourier basis. We found k5 \+
to be related to k1." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 53 "plo
t(lhs(BC2),k=0..9,view=[0..9,-4..4],discont=true);" }}}{EXCHG {PARA 0 
"> " 0 "" {MPLTEXT 1 0 3 "k6;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 
0 27 "k7:=fsolve(BC2,k=4.1..4.8);" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 27 "k8:=fsolve(BC2,k=4.8..5.5);" }}}{EXCHG {PARA 0 "> " 
0 "" {MPLTEXT 1 0 25 "k9:=fsolve(BC2,k=5.5..6);" }}}{EXCHG {PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 26 "k10:=fsolve(BC2,k=6..6.5);" }}}{EXCHG {PARA 0 
"> " 0 "" {MPLTEXT 1 0 20 "k11:=k1+evalf(2*Pi);" }}}{EXCHG {PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 20 "k12:=k2+evalf(2*Pi);" }}}{EXCHG {PARA 0 "> " 
0 "" {MPLTEXT 1 0 20 "k13:=k3+evalf(2*Pi);" }}}{EXCHG {PARA 0 "> " 0 "
" {MPLTEXT 1 0 20 "k14:=k4+evalf(2*Pi);" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 20 "k15:=k5+evalf(2*Pi);" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 20 "k16:=k6+evalf(2*Pi);" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 42 "plot([seq([n,k||n],n=1..16)],style=point);" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 33 "k1,k2-k1,k3-k2,k4-k3,k5-k4,k
6-k5;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 78 "We can actually infer th
e behaviour of the eigenvalues as a function of order:" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 15 "evalf(2*Pi/10);" }}}{EXCHG {PARA 0 
"" 0 "" {TEXT -1 118 "Now we are able to carry out a limited convergen
ce study. Let us pick an arbitrary function, and expand it over [a,b].
" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 22 "f2:=exp(-x/15)*sin(x);
" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 48 "Bfn:=[seq(Bf(i)/sqrt(IP
(Bf(i),Bf(i))),i=1..15)]:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 
34 "c_n:=[seq(IP(f2,Bfn[j]),j=1..15)];" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 34 "f2_ap5:=add(c_n[n]*Bfn[n],n=1..5):" }}}{EXCHG {PARA 
0 "> " 0 "" {MPLTEXT 1 0 36 "f2_ap10:=add(c_n[n]*Bfn[n],n=1..10):" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 36 "f2_ap15:=add(c_n[n]*Bfn[n],n
=1..15):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 70 "plot([f2_ap5,f2
_ap10,f2_ap15,f2],x=a..b,color=[red,blue,green,black]);" }}}{EXCHG 
{PARA 0 "" 0 "" {TEXT -1 110 "What happens at the left-hand side is kn
own as the Gibbs phenomenon. Given that all basis functions vanish at \+
" }{TEXT 19 3 "x=a" }{TEXT -1 8 " and at " }{TEXT 19 3 "x=b" }{TEXT 
-1 436 ", the approximate Fourier series representation cannot be accu
rate there. The Fourier expansion tries to approximate a step function
 behaviour: at the endpoints the function is assumed to rise immediate
ly to the non-zero value with infinite slope. Fourier expansions with \+
a large number of terms are able to produce a large slope, but the pri
ce is paid by the wiggling around the true function even some distance
 away from the boundary." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "
" 0 "" {TEXT 261 11 "Exercise 4:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}
{PARA 0 "" 0 "" {TEXT -1 52 "Use the analytic expression found for the
 values of " }{TEXT 19 4 "k||n" }{TEXT -1 44 " to push the expansion t
o high order (e.g., " }{TEXT 19 8 "n_max=50" }{TEXT -1 35 ") and obser
ve the Gibbs phenomenon." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG 
{PARA 0 "" 0 "" {TEXT 260 35 "Solution of a differential equation" }}
{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 287 "So far w
e have not calculated anything particularly useful given that the func
tions to be represented could be easily expressed as elementary expres
sions. We will show now that the Fourier series expansion allows to co
nvert differential equation problems into linear algebra assignments.
" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 142 "Let \+
us define an inhomogeneous linear second-order differential equation w
ith boundary conditions specified at the ends of the interval [a,b]." 
}}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 51 "DE:=diff(u(x),x$2)+1/4*x*
diff(u(x),x)+x^2*u(x)=x^4;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 
19 "BCd:=u(a)=0,u(b)=0;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 29 "
solDE:=dsolve(\{DE,BCd\},u(x));" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 
106 "Maple can't solve this one. We can generate a numerical solution \+
to this boundary value problem using the " }{TEXT 19 15 "dsolve[numeri
c]" }{TEXT -1 12 " integrator." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 
1 0 58 "solDE:=dsolve(\{DE,BCd\},u(x),numeric,output=listprocedure):" 
}}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "U:=subs(solDE,u(x)):" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 18 "plot(U(x),x=a..b);" }}}
{EXCHG {PARA 0 "" 0 "" {TEXT -1 217 "Now we produce a basis-expansion \+
solution. We generate a matrix from the left-hand side of the equation
 by acting with the differential operator on basis function Bfn[j], wh
ile multiplying with Bfn[i] and integrating." }}}{EXCHG {PARA 0 "> " 
0 "" {MPLTEXT 1 0 8 "lhs(DE);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 
0 6 "Nb:=9;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "with(LinearA
lgebra):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "A:=Matrix(Nb,Nb
):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 99 "for i from 1 to Nb do
: for j from 1 to Nb do: A[i,j]:=IP(Bfn[i],subs(u(x)=Bfn[j],lhs(DE))):
 od: od:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 9 "evalm(A);" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 21 "Ai:=MatrixInverse(A);" }}}
{EXCHG {PARA 0 "" 0 "" {TEXT -1 100 "Now we form a right-hand-side vec
tor using a projection of the rhs of the ODE onto the basis states:" }
}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "bvec:=Vector(Nb):" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 55 "for i from 1 to Nb do: bvec[
i]:=IP(rhs(DE),Bfn[i]): od:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 
12 "evalm(bvec);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "cvec:= \+
Ai . bvec;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 35 "solAP:=add(cv
ec[i]*Bfn[i],i=1..Nb):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 43 "p
lot([solAP,U(x)],x=a..b,color=[red,blue]);" }}}{EXCHG {PARA 0 "" 0 "" 
{TEXT 262 11 "Exercise 5:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "
" 0 "" {TEXT -1 289 "Pick a linear differential equation of your choic
e with appropriate boundary conditions and compare numerical solutions
 as provided by Maple with the Fourier-series solution. Observe the ef
ficiency of the Fourier method: the solution is specified using a smal
l number of real coefficients." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}
{PARA 0 "" 0 "" {TEXT 263 56 "Solution of a differential equation - ei
genvalue problem" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" 
{TEXT -1 344 "Here we apply the Fourier-series expansion technique to \+
a Schroedinger equation. This will result in a matrix eigenvalue probl
em. We will continue to use our asymmetric interval [a,b], but apply i
t to a situation where a symmetric interval would be more appropriate.
 We cosnsider the bound states in a potential that has an inverted bel
l shape:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 23 "Vpot:=-2/(1+(2*
x)^2)^2;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 18 "plot(Vpot,x=a..
b);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 87 "We will ask the following \+
question: Given that the kinetic energy can be represented as" }}
{PARA 0 "" 0 "" {TEXT 19 29 "Tkin = - 1/2 diff(psi(x),x$2)" }}{PARA 0 
"" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 44 "what are the allo
wed energy eigenvalues for " }{TEXT 19 9 "Tkin+Vpot" }{TEXT -1 67 "? [
for physicists: we have chosen units in which Planck's constant " }
{TEXT 19 6 "h=2*Pi" }{TEXT -1 179 ", and the mass of the particle equa
ls unity, and a third constant was chosen so that the strength of the \+
potential energy is dimensonless]. The kinetic energy operator represe
nts " }{TEXT 19 9 "p^2/(2*m)" }{TEXT -1 43 " using the association (pr
escription) that " }{TEXT 19 21 "p = - I*h/(2*Pi)*d/dx" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 49 "SE:=-1/2*diff(psi(x),x$2)+Vpot*psi(
x) = E*psi(x);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 42 "We are looking \+
for bound solutions, i.e., " }{TEXT 19 5 "E < 0" }{TEXT -1 174 ". The \+
ODE problem is homogeneous (eigenvalue problems always are). This prob
lem cannot be solved easily by dsolve[numeric], because the energy eig
envalue is unknown a priori." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 
0 "" 0 "" {TEXT -1 46 "We generate a matrix representation as before:
" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 6 "Nb:=5;" }}}{EXCHG {PARA 
0 "> " 0 "" {MPLTEXT 1 0 17 "A:=Matrix(Nb,Nb):" }}}{EXCHG {PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 121 "for i from 1 to Nb do: for j from 1 to Nb do:
 A[i,j]:=IP(Bfn[i],subs(psi(x)=Bfn[j],lhs(SE))): print(i,j,A[i,j]);  o
d: od:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 9 "evalm(A);" }}}
{EXCHG {PARA 0 "" 0 "" {TEXT -1 18 "What do we notice?" }}{PARA 0 "" 
0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 70 "1) the imaginary par
t in the matrix elements looks like spurious noise" }}{PARA 0 "" 0 "" 
{TEXT -1 125 "2) the matrix is almost symmetric [the small deviations \+
are a result of using non-symmetric boundary conditions in the basis]
" }}{PARA 0 "" 0 "" {TEXT -1 75 "Before we proceed with finding the ei
genvalues we correct the two problems:" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 13 "A:=map(Re,A);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 
1 0 26 "A:=0.5*(A + Transpose(A));" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 15 "Eigenvalues(A);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 
143 "We have found one bound state. To determine how reliable the ener
gy estimate is we have to increase the matrix size and repeat the calc
ulation." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT 19 
12 "E0 = -.52511" }{TEXT -1 14 " from 4-by-4 (" }{TEXT 19 4 "Nb=4" }
{TEXT -1 1 ")" }}{PARA 0 "" 0 "" {TEXT 19 12 "E0 = -.54981" }{TEXT -1 
14 " from 9-by-9 (" }{TEXT 19 4 "Nb=9" }{TEXT -1 1 ")" }}{PARA 0 "" 0 
"" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT 264 11 "Exercise 6:" }}{PARA 
0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 143 "Carry out a co
nvergence analysis. Concentrate on the lowest-lying eigenvalue (the gr
ound state). The higher-lying states are artificial states." }}{PARA 
0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT 265 11 "Exercise 7:" }
}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 173 "Deepen \+
the potential well by increasing the constant that pre-multiplies the \+
potential. Re-compute the ground state. Are there additional bound sta
tes in the potential well?" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 
"" 0 "" {TEXT -1 29 "To compute the eigenfunction " }{TEXT 19 8 "psi_0
(x)" }{TEXT -1 33 " associated with the eigenenergy " }{TEXT 19 2 "E0
" }{TEXT -1 102 " we have to look at the eigenvector and assemble the \+
superposition. The syntax is somewhat cumbersome." }}}{EXCHG {PARA 0 "
> " 0 "" {MPLTEXT 1 0 20 "EV:=Eigenvectors(A);" }}}{EXCHG {PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 6 "n0:=1;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 
0 10 "EV[1][n0];" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 29 "cf:=seq
(EV[2][j,n0],j=1..Nb);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 36 "p
si0:=add(Re(cf[j]*Bfn[j]),j=1..Nb):" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 20 "plot(psi0^2,x=a..b);" }}}{EXCHG {PARA 0 "" 0 "" 
{TEXT -1 313 "This is an approximation to the square of the ground-sta
te wavefunction which (when normalized properly) carries the probabili
ty information of where to find particles in position space (if they a
re in the ground state). Observe carefully how the wings of the distri
bution are in the classically forbidden regime:" }}}{EXCHG {PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 70 "plot([Vpot,psi0^2+EV[1][n0],EV[1][n0]],x=a..b,
color=[red,blue,green]);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 165 "The \+
baseline for the wavefunction graph intersects the potential energy at
 the classical turning point (all kinetic energy is converted into pot
ential energy there)." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 
"" {TEXT -1 59 "The asymmetry is caused by the basis. For larger value
s of " }{TEXT 19 2 "Nb" }{TEXT -1 92 " the ground state is calculated \+
to be symmetric (as it should be for a symmetric potential)." }}{PARA 
0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 42 "The state is ac
tually properly normalized:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 
14 "IP(psi0,psi0);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 247 "For suffic
iently deep (and wide) potentials additional bound states appear. The \+
next bound state (first excited state) in a symmetric potential is ant
i-symmetric (has odd parity). The probability density (square of the w
avefunction) has a node at " }{TEXT 19 3 "x=0" }{TEXT -1 64 " in this \+
case. Verify this feature on an example of your choice." }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}{MARK "0 0 0" 0 }{VIEWOPTS 1 
1 0 1 1 1803 1 1 1 1 }{PAGENUMBERS 0 1 2 33 1 1 }