{VERSION 6 0 "IBM INTEL NT" "6.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 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 263 "" 0 1 0 0 0 0 0 1 0 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 } {CSTYLE "" -1 266 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 267 "" 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 268 "" 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 269 "" 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 270 "" 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 20 "Multidimensional FFT" }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 171 "A test o f Maple's new discrete Fourier transform package (available as of Mapl e9). First, we test it in one dimension and demonstrate the effect of \+ two mesh distributions." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 507 "First case: a one-dimensional mesh that includes th e origin, and is similar in both the coordinate and momentum domain. W e restrict ourselves to meshes with size in powers of two (as in the o ld FFT program), even though the algorithm is capable of handling othe r meshes as well. We will use the interpretation of the Fourier transf orm as familiar from quantum mechanics, i.e., connecting the coordinat e and momentum representations. The alternative would be to look at ti me signals in the frequency domain. " }}{PARA 0 "" 0 "" {TEXT -1 160 " We do not overemphasize the QM aspect, i.e., we work in units where Pl anck's constant (hbar) equals unity. Our domains of integration will b e symmetric both in " }{TEXT 19 1 "x" }{TEXT -1 8 " and in " }{TEXT 19 1 "p" }{TEXT -1 1 "." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT 263 22 "1. One-dimensional FFT" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "restart; Digits:=15; " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "with(DiscreteTransforms );" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 6 "N:=64;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "dx:=0.3;" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 24 "dp:=evalhf(2*Pi/(N*dx));" }}}{EXCHG {PARA 0 "" 0 " " {TEXT -1 120 "Our choice of parameters resulted in a similar spacing in coordinate and momentum. We pick a function to be transformed." }} }{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 16 "fx:=1/(1+x^2)^2;" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 35 "Z:=Array(1..N,datatype=compl ex[8]);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 327 "In the first 1d imple mentation we introduce a coordinate mesh which includes the origin. Th e mesh is arranged in such a way that one starts with zero, goes to th e largest represented positive value, then continues with the largest \+ negative value (which is one unit further than the largest positive) a nd returns to the origin. " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 23 "for i from 1 to N/2 do:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 162 "pv[ i]:=evalhf((i-1)*dp); pv[N-i+1]:=evalhf(-i*dp); xv[i]:=evalhf((i-1)*dx ); Z[i]:=eval(fx,x=xv[i]); xv[N-i+1]:=evalhf(-i*dx); Z[N-i+1]:=eval(f x,x=xv[N-i+1]): od:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 56 "plot ([seq([xv[i],Z[i]] ,i=1..N)],x=-10..10,style=point);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 75 "Note that our mesh spends a lot of points where the function is very small." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 90 "Now we carry out the Fourier transform whereby \+ the result is stored in the original array:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 33 "FourierTransform(Z,inplace=true);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 15 "Z[1],Z[2],Z[N];" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 149 "Note the negligible imaginary parts. The discret e FT recognized the symmetry of the function. We graph both real and i maginary parts to confirm this:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 50 "plot([seq([pv[i],Re(Z[i])] ,i=1..N)],style=point);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 50 "plot([seq([pv[i],Im(Z[i])] ,i=1..N) ],style=point);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 50 "How well does \+ the DFT agree with the known answer?" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 61 "fp:=int(fx*exp(I*p*x),x=-infinity..infinity) assuming p,real;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 21 "evalf(limit(fp, p=0));" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 19 "plot(fp,p=-10..10 );" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 5 "Z[1];" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 140 "Before we can compare in detail we have \+ to resolve some normalization issues. Neither our integral, nor the Ma ple DFT contain the customary " }{TEXT 19 12 "1/sqrt(2*Pi)" }{TEXT -1 32 " factor. Maple's DFT contains a " }{TEXT 19 9 "1/sqrt(N)" }{TEXT -1 26 " factor, but is missing a " }{TEXT 19 2 "dx" }{TEXT -1 103 " fa ctor to form the Riemann sum (which connects the result to the integra l). Thus, we correct for both:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 37 "f_DFT:=evalf(sqrt(N)*dx): f_DFT*Z[1];" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 60 "P1:=plot([seq([pv[i],Re(f_DFT*Z[i])] ,i=1..N)],s tyle=point):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 57 "P2:=plot(fp ,p=-10..10,color=blue): plots[display](P1,P2);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 27 "The agreement is very good." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT 257 13 "Exercise 1.1:" }}{PARA 0 "" 0 "" {TEXT -1 121 "Explore the question of agreement/disagreement in detail: calculate the percentage error and display it for all mome nta (" }{TEXT 19 1 "p" }{TEXT -1 9 "-values)." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT 258 13 "Exercise 1.2:" }}{PARA 0 "" 0 "" {TEXT -1 48 "Explore the behaviour using different values of " }{TEXT 19 6 "N, dx:" }{TEXT -1 46 " can one obtain acceptable resul ts with small " }{TEXT 19 1 "N" }{TEXT -1 1 "?" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT 259 13 "Exercise 1.3:" }}{PARA 0 "" 0 "" {TEXT -1 25 "Modify the function pair " }{TEXT 19 8 "\{fx, f p\}" }{TEXT -1 22 " (you can replace the " }{TEXT 19 1 "1" }{TEXT -1 141 " in the denominator by a constant), or change to a different pair (different power of the entire denominator?) and explore Exercises 1 \+ and 2." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT 260 13 "Mini-project:" }{TEXT -1 1 " " }}{PARA 0 "" 0 "" {TEXT -1 129 "Exp lore further issues with the FFT: calculate the second derivate of the function by going to transform space, multiplying with " }{TEXT 19 3 "p^2" }{TEXT -1 30 ", and transforming back using " }{TEXT 19 23 "Inve rseFourierTransform" }{TEXT -1 41 ", and compare with the analytical r esult." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 70 "Now we try a modification of the mesh that excl udes the origin in the " }{TEXT 19 1 "x" }{TEXT -1 80 "-domain. This i s useful for cases when the original function is undefined there." }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 20 "restart; Digits:=15:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "with(DiscreteTransforms);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 6 "N:=64;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 35 "Z :=Array(1..N,datatype=complex[8]);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "dx:=0.3;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 14 "xM:=Vector(N):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 65 "for i fr om 1 to N/2 do: xM[i]:=(i-1/2)*dx; xM[N-i+1]:=-xM[i]; od:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 30 "xM[1],xM[N/2],xM[N/2+1],xM[N];" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 24 "dpx:=evalf(2*Pi/(N*dx));" }} }{EXCHG {PARA 0 "" 0 "" {TEXT -1 71 "It will turn out that the transfo rmed variable domain does contain the " }{TEXT 19 3 "p=0" }{TEXT -1 7 " point:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 44 "for i from 1 to N/2+1 do: pxM[i]:=(i-1)*dpx:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 41 "if i>1 then pxM[N-i+2]:=-pxM[i]: fi: od: " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 41 "pxM[1],pxM[2],pxM[N/2],pxM[N/2+1],pxM[N];" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 5 "a:=1;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 18 "fx :=1/(a^2+x^2)^2;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 61 "fp:=int (fx*exp(I*p*x),x=-infinity..infinity) assuming p,real;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 31 "for i from 1 to N do: x:=xM[i]:" }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 21 "Z[i]:=fx: od: x:='x':" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 33 "FourierTransform(Z,inplace=true);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 10 "Z[1],Z[2];" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 27 "Why is this complex-valued?" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 73 "Arrange the density in an array with p ra nging from negative to positive." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 50 "plot([seq([pxM[i],Re(Z[i])],i=1..N)],style=point);" } }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 50 "plot([seq([pxM[i],Im(Z[i]) ],i=1..N)],style=point);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "f_DFT:=evalf(sqrt(N)*dx);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 60 "P1:=plot([seq([pxM[i],Re(f_DFT*Z[i])],i=1..N)],style=point):" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 57 "P2:=plot(fp,p=-10..10,color= blue): plots[display](P1,P2);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 219 "As could have been guessed, the real part alone does not agree with t he answer. If we take the absolute value, matters are corrected. Thus, for this mesh arrangement we obtain a DFT which contains a wrong comp lex phase." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 61 "P1:=plot([seq ([pxM[i],abs(f_DFT*Z[i])],i=1..N)],style=point):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 57 "P2:=plot(fp,p=-10..10,color=blue): plots[disp lay](P1,P2);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 184 "Therefore, it mi ght be better to work with the mesh that contains the origin, and to u se a regulator in those instances where the function to be transformed has problems at the origin." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT 261 13 "Exercise 1.4:" }}{PARA 0 "" 0 "" {TEXT -1 50 " Explore the effect of the phase error by reducing " }{TEXT 19 2 "dx" } {TEXT -1 19 ". You can increase " }{TEXT 19 1 "N" }{TEXT -1 78 " or ke ep it constant. What happens to the real and imaginary parts of the DF T?" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT 265 13 "Ex ercise 1.5:" }}{PARA 0 "" 0 "" {TEXT -1 132 "An interesting example to Fourier analyze is the Gaussian wavepacket which one can control thro ugh various parameters such as width " }{TEXT 19 2 "Dx" }{TEXT -1 99 " (coordinate space parameter corresponds to inverse momentum width par ameter), and center momentum " }{TEXT 19 2 "p0" }{TEXT -1 163 ". It is relatively straightforward to find parameter choices which exceed the mesh capabilities, either by exceeding the momentum resolution (very \+ broad packet in " }{TEXT 19 1 "x" }{TEXT -1 29 "-space = narrow struct ure in " }{TEXT 19 1 "p" }{TEXT -1 23 "-space), or by putting " } {TEXT 19 2 "p0" }{TEXT -1 70 " at the edge of the covered momentum ran ge. Explore this on your own. " }{TEXT 19 35 "GWP=Norm*exp(-(x/Dx)^2)* exp(I*p0*x)" }{TEXT -1 66 ". Find the analytic Fourier image and compa re with the FFT result." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 " " }}}{EXCHG {PARA 0 "" 0 "" {TEXT 262 22 "2. Two-dimensional FFT" }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 192 "Now we p roceed with the two-dimensional case. Comparison with analytical resul ts will be more complicated. In contrast to the old FFT routines multi -dimensional FFTs are calculated within the " }{TEXT 19 18 "DiscreteTr ansforms" }{TEXT -1 9 " package." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 46 "restart; Digits:=15: with(DiscreteTransforms);" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 7 "N:=128;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 40 "Z:=Array(1..N,1..N,datatype=complex[8]);" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 16 "dx:=0.3; dy:=dx:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 62 "The simplest function to choose would be \+ one which factorizes." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 5 "a:= 1;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 18 "fx:=1/(a^2+x^2)^3;" } }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 64 "fpx:=int(fx*exp(I*px*x),x= -infinity..infinity) assuming px,real;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 5 "b:=2;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 18 "fy :=1/(b^2+y^2)^3;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 64 "fpy:=in t(fy*exp(I*py*y),y=-infinity..infinity) assuming py,real;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 29 "xM:=Vector(N): yM:=Vector(N):" }}}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 104 "for i from 1 to N/2 do: xM[i]:=(i-1/2)*dx; x M[N-i+1]:=-xM[i]; yM[i]:=(i-1/2)*dy; yM[N-i+1]:=-yM[i]; od:" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 42 "xM[1],xM[N/2],xM[N/2+1],xM[N ],yM[1],yM[N];" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 125 "We chose a cub e, but nothing prevented us from using a parallelepiped as a computati onal domain. Simply use different dy, dz." }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 63 "for i from 1 to N do: x:=xM[i]: for j from 1 to N d o: y:=yM[j]:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 40 "r:=sqrt(x^2+y^2): Z [i,j]:=fx*fy: od: od:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 44 "DA :=Array(1..N,1..N,1..3,datatype=float[8]):" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 119 "#for i from 1 to N do: x:=xM[i]: for j from 1 to N do: y:=yM[j]: DA[i,j,1]:=x: DA[i,j,2]:=y: DA[i,j,3]:=Z[i,j]; od: od: " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 128 "for i from 1 to N do: \+ x:=xM[i]: for j from 1 to N do: y:=yM[j]: DA[i,j,1]:=x: DA[i,j,2]:=y: \+ DA[i,j,3]:=log(abs(Z[i,j])); od: od:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 4 "The " }{TEXT 19 8 "surfdata" }{TEXT -1 85 " procedure is terribl e at plotting the original surface. That's why we graph the log." }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 63 "plots[surfdata](DA,axes=boxe d,style=patchcontour,shading=zhue);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 33 "FourierTransform(Z,inplace=true);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 7 "Z[1,2];" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 123 "Why is this complex-valued? The one-dimensional results imply tha t the mesh distribution without an origin in the original " }{TEXT 19 1 "r" }{TEXT -1 22 "-domain are to blame." }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 24 "dpx:=evalf(2*Pi/(N*dx));" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 24 "dpy:=evalf(2*Pi/(N*dy)):" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 63 "for i from 1 to N/2+1 do: pxM[i]:=(i-1)*dpx: pyM[i] :=(i-1)*dpy:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 62 "if i>1 then pxM[N-i +2]:=-pxM[i]: pyM[N-i+2]:=-pyM[i]: fi: od: " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 41 "pxM[1],pxM[2],pxM[N/2],pxM[N/2+1],pxM[N];" }}} {EXCHG {PARA 0 "" 0 "" {TEXT -1 190 "At one point we were so dismayed \+ by the bad surface plot that we exported the data, and used Surfer 7 t o graph them. The grid generation and contour plot in that program is \+ much superior to " }{TEXT 19 8 "surfdata" }{TEXT -1 76 ". We left the \+ commented commands to show how to export data via ASCII files." }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 51 "# fd := fopen(\"C:/marko/map le9/testsrf.dat\",WRITE);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 44 "DA:=Array(1..N,1..N,1..3,datatype=float[8]):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 129 "for i from 1 to N do: px:=pxM[i]: for j from 1 to N \+ do: py:=pyM[j]: DA[i,j,1]:=px: DA[i,j,2]:=py: DA[i,j,3]:=abs(Z[i,j]); \+ od: od:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 19 "px:='px': py:='p y':" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 44 "# fprintf(fd,\"%g %g %g\\n\",px,py,DA[i,j,3]); " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 123 "for i from 1 to N do: for j from 1 to N do: if DA[i,j,3]>0 then D A[i,j,3]:=log10(DA[i,j,3]): else DA[i,j,3]:=0: fi: od: od:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 181 "We cheat here: the zero values that occu r at the largest momentum value are transferred to the origin, where w e enter the known value. Didn't know how else to get rid off the value ." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 125 "for i from 1 to N do: for j from 1 to N do: if DA[i,j,3]=0 then DA[i,j,1]:=0: DA[i,j,2]:=0; DA[i,j,3]:=DA[1,1,3]; fi: od: od:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 13 "# fclose(fd);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 63 "plots[surfdata](DA,axes=boxed,style=patchcontour,shading=zhue) ;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 81 "Why the surface looks differ ent viewed from the top versus the bottom is for the " }{TEXT 19 8 "su rfdata" }{TEXT -1 25 " programmers to explain. " }}{PARA 0 "" 0 "" {TEXT -1 102 "There are cases when the graph looks bad from the top, b ut has OK contours when looked at from below. " }{TEXT 19 11 "FT(exp(- r))" }{TEXT -1 15 " is an example." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT -1 139 "How do we know the calculated 2dimension al FT is correct? We would have to use our known FT pair and plug in t he correct conversion factor." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 37 "f_DFT2:=evalf(sqrt(N)*dx*sqrt(N)*dy);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 14 "Z[1,1]*f_DFT2;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 31 "evalf(eval(fpx*fpy,\{x=0,y=0\}));" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 100 "plot3d(log10(abs(fpx*fpy/f_DFT2)),px=-10.. 10,py=-10..10,axes=boxed,style=patchcontour,shading=zhue);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 38 "This looks OK what about the details?." } }{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "n_x:=15; n_y:=10;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 " DA[n_x,n_y];" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 66 "evalf(eval( log10(abs(fpx*fpy/f_DFT2)),\{px=pxM[n_x],py=pyM[n_y]\}));" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 59 "The agreement at intermediate momentum va lues is excellent!" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT 264 13 "Exercise 2.1:" }}{PARA 0 "" 0 "" {TEXT -1 95 "Find other FT function pairs (ideally non-factorized ones) and observe the accur acy of the DFT." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}} {EXCHG {PARA 0 "" 0 "" {TEXT 266 24 "3. Three-dimensional FFT" }{TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 46 "restart; Digits:= 15: with(DiscreteTransforms);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 6 "N:=64;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 45 "Z:=Array(1.. N,1..N,1..N,datatype=complex[8]);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 24 "dx:=0.3; dy:=dx: dz:=dy:" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 44 "xM:=Vector(N): yM:=Vector(N): zM:=Vector(N):" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 143 "for i from 1 to N/2 do: xM[ i]:=(i-1/2)*dx; xM[N-i+1]:=-xM[i]; yM[i]:=(i-1/2)*dy; yM[N-i+1]:=-yM[ i]; zM[i]:=(i-1/2)*dz; zM[N-i+1]:=-zM[i]; od:" }}}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 36 "xM[1],xM[N],yM[1],yM[N],zM[1],zM[N];" }}} {EXCHG {PARA 0 "" 0 "" {TEXT -1 118 "We chose a cube, but nothing prev ented us from using a parallelepiped as a computational domain. Simply use different " }{TEXT 19 6 "dy, dz" }{TEXT -1 1 "." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 5 "a:=1;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 95 "for i from 1 to N do: x:=xM[i]: for j from 1 to N do: y:=yM[j]: for k from 1 to N do: z:=zM[k]:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 54 "r:=sqrt(x^2+y^2+z^2): Z[i,j,k]:=exp(-a*r): od: od: od :" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 33 "FourierTransform(Z,inp lace=true);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 9 "Z[1,1,2];" }} }{EXCHG {PARA 0 "" 0 "" {TEXT -1 27 "Why is this complex-valued?" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 24 "dpx:=evalf(2*Pi/(N*dx));" }} }{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 49 "dpy:=evalf(2*Pi/(N*dy)): dp z:=evalf(2*Pi/(N*dz)):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 82 "f or i from 1 to N/2+1 do: pxM[i]:=(i-1)*dpx: pyM[i]:=(i-1)*dpy: pzM[i]: =(i-1)*dpz:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 83 "if i>1 then pxM[N-i+ 2]:=-pxM[i]: pyM[N-i+2]:=-pyM[i]: pzM[N-i+2]:=-pzM[i]: fi: od: " }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 41 "pxM[1],pxM[2],pxM[N/2],pxM[N /2+1],pxM[N];" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 92 "Now look at a cu t. Arrange the density in an array with p ranging from negative to pos itive." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 44 "DA:=Array(1..N,1. .N,1..3,datatype=float[8]):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 198 "for i from 1 to N do: px:=pxM[i]: for j from 1 to N do: py:=pyM[j]: DA[i, j,1]:=px: DA[i,j,2]:=py: DA[i,j,3]:=0: for k from 1 to N do: pz:=pzM[ k]: DA[i,j,3]:=DA[i,j,3]+abs(Z[i,j,k])*dpz; od: od: od:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 123 "for i from 1 to N do: for j from 1 to N do: if DA[i,j,3]>0 then DA[i,j,3]:=log10(DA[i,j,3]): else DA[i,j ,3]:=0: fi: od: od:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 141 "for i from 1 to N do: for j from 1 to N do: if DA[i,j,3]=0 then DA[i,j,3] :=DA[1,1,3]: DA[i,j,1]:=DA[1,1,1]: DA[i,j,2]:=DA[1,1,2]: fi: od: od:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 63 "plots[surfdata](DA,axes=b oxed,style=patchcontour,shading=zhue);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 65 "This is a case that looks better from the bottom than fro m above." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 48 "f_DFT3:=evalf(s qrt(N)*dx*sqrt(N)*dy*sqrt(N)*dz);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 170 "How can we verify the validity? We need to calculate the analytic FT for our function. This is done in spherical polar coordinates. The azimuthal angular integral yields " }{TEXT 19 4 "2*Pi" }{TEXT -1 87 " . The polar angle integration is straightforward. The radial integral \+ is done in Maple:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 82 "r:='r' : fp:=simplify(4*Pi/p*int(exp(-a*r)*sin(p*r)*r,r=0..infinity)) assumin g p>0;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 126 "We have a spherically \+ symmetric answer. We should be able to verify the FT result by choosin g a line along a beam starting at " }{TEXT 19 3 "p=0" }{TEXT -1 1 "." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 33 "ResN:=Vector(N): ResA:=Ve ctor(N):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 151 "for i from 1 t o N/2 do: px:=pxM[i]: py:=pyM[i]: pz:=pzM[i]: p:=sqrt(px^2+py^2+pz^2): ResA[i]:=[p,p*evalf(fp)]: ResN[i]:=[p,p*abs(f_DFT3*Z[i,i,i])]: od:" } }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 16 "ResA[2],ResN[2];" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 46 "P1:=plot([seq(ResN[i],i=1..N /2)],style=point):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 68 "P2:=p lot([seq(ResA[i],i=1..N/2)],color=blue): plots[display](P1,P2);" }}} {EXCHG {PARA 0 "" 0 "" {TEXT -1 81 "The FFT works remarkably well give n that we have not resolved the maximum in the " }{TEXT 19 6 "p*f(p)" }{TEXT -1 64 " distribution, i.e., we are working on a relatively coar se mesh." }}{PARA 0 "" 0 "" {TEXT -1 127 "Keep in mind that we used th e absolute value of the transform. A mesh including the origin should \+ not suffer from this problem." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT 267 13 "Exercise 2.2:" }}{PARA 0 "" 0 "" {TEXT -1 21 "Change the parameter " }{TEXT 19 1 "a" }{TEXT -1 213 " in the w avefunction and observe the change in the momentum distribution. Expla in what it means physically. Find other function pairs (separable in C artesians, and non-separable) and verify how well the FFT works." }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT 269 13 "Exercise 2.3:" }}{PARA 0 "" 0 "" {TEXT -1 149 "For QM-literate students: consi der the calculation of the kinetic energy expectation value in momentu m space. You can calculate it by multiplying in " }{TEXT 19 1 "p" } {TEXT -1 40 "-space with the multiplicative operator " }{TEXT 19 5 "p^ 2/2" }{TEXT -1 15 " (in our units " }{TEXT 19 8 "m=hbar=1" }{TEXT -1 81 ") and forming the average (and normalizing properly). This will em phasize larger " }{TEXT 19 1 "p" }{TEXT -1 87 "-values. To check how a ccurate this calculation will be you can probe the agreement of " } {TEXT 19 8 "p^2*f(p)" }{TEXT -1 77 " in the above plot. What do you fi nd? How can you improve upon the situation?" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT 270 13 "Exercise 2.4:" }}{PARA 0 "" 0 "" {TEXT -1 88 "Change the mesh in the above code to include the origi n also in the configuration space " }{TEXT 19 7 "(x,y,z)" }{TEXT -1 95 " domain. Demonstrate that the problems with the incorrect complex \+ phase disappear in this case." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT 268 13 "Mini-project:" }}{PARA 0 "" 0 "" {TEXT -1 118 "Demonstrate the usefulness of the 3d FFT by making other cuts \+ for comparison between numerical and analytical results." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}{MARK "146" 0 }{VIEWOPTS 1 1 0 1 1 1803 1 1 1 1 }{PAGENUMBERS 0 1 2 33 1 1 }