2
|
332 Views
|
2 Replies
|
5 Total Likes
View groups...
Share

# Help needed with 2D Schrödinger equation using discrete Fourier transform and Split-Operator method

Posted 1 month ago
 Hello, dear forum members,I have a question that I've been struggling with for some time, and despite extensive searching, I haven't been able to find an answer or identify the error.I’m working with the discrete Fourier transform to solve the one- and two-dimensional Schrödinger equation using the split-operator method. This involves splitting the Hamiltonian into kinetic and potential parts. As long as the time steps are sufficiently small, the split-operator method is generally a suitable approach for determining the discrete wave function. In the one-dimensional case, this method works very well compared to others (e.g., Crank-Nicolson, see attached file). However, I’m encountering issues in the two-dimensional case, where I'm not obtaining the expected approximate wave function. I believe the error lies in the transformation to momentum space or I may have overlooked something basic.In my example, a free wave packet with momentum k in the y-direction splits into two parts after a few time steps, which doesn't seem correct.Could someone please help me with this issue? The example is taken from: J. M. Feagin, Quantum Methods with Mathematica, New York, TELOS/Springer-Verlag, 1994.Thank you in advance! TDFourier2D[timesteps_, potheight_, ellip_] := Module[{NP, L, dx, dt, kx, ky, x0, y0, xmax, xmin, ymax, ymin, kxmax, kymax, \[Rho], TimeEnd, psi, K3, nt}, NP = 100; (* netpoints *); L = 1; dx = L/(NP - 1); dt = 1/20 dx^2; TimeEnd = timesteps; \[Rho] = 0.16; (* initial width *) x0 = 0.0 ; (* initial position *) y0 = 0.0; (* initial position *) kx = 0 ;(* wave number x *) ky = 50; (* wave number y *) xmax = 1; (* dimension of the box *) xmin = -1; (* dimension of the box *) ymax = 1; (* dimension of the box *) ymin = -1; (* dimension of the box *) kxmax = N[Pi/dx]; (* momentum lattice x-direction*); kymax = N[Pi/dx]; psi3[0] = Array[1/(\[Rho] Sqrt[Pi]) Exp[-I kx #1 - (#1 - x0)^2/(2 \[Rho]^2)]* Exp[-I ky #2 - (#2 - y0)^2/(2 \[Rho]^2)] &, {NP, NP}, {{xmin, xmax}, {ymin, ymax}}] // Quiet; K3 = Array[ 1/2 #1^2 + 1/2 #2^2 &, {NP, NP}, {{-kxmax, kxmax}, {-kymax, kymax}}]; V3 = Array[ potheight/2 #1^2 + ellip*potheight/2 #2^2 &, {NP, NP}, {{xmin, xmax}, {ymin, ymax}}]; UV[dt] = Exp[-I V3 dt]; UK[dt] = Exp[-I K3 dt]; U[dt_]@psi_ := UV[dt]*InverseFourier[UK[dt]*Fourier[psi]]; psi3[nt_, psi0_] := Nest[U[dt][#] &, psi0, nt]; Do[psi3[nt] = psi3[1, psi3[nt - 1]], {nt, 1, TimeEnd, 1}]; ] ellipticalform = 0.2; (*elliptical form from from > 0 and <= 1; 1=cylindrical*) TimeEnd = 160; (* timesteps 10 to 150 *) potheight = 0; (*potential height, free particle *) (TDFourier2D[TimeEnd, potheight, ellipticalform]) // AbsoluteTiming graphTable1 = GraphicsGrid[ Partition[ Table[Labeled[ ListDensityPlot[Abs[psi3[i]], Axes -> True, PlotRange -> All, ColorFunction -> "SunsetColors", PerformanceGoal -> "Quality"], {"x", "y"}, {Left, Top}, RotateLabel -> True], {i, 1, TimeEnd, TimeEnd/4}], 2], Frame -> All, ImageSize -> 400]  Attachments:
2 Replies
Sort By:
Posted 23 days ago
 Thank you very much for your prompt assistance. Reducing the grid steps in time and space does help improve the numerical calculation of the discrete Schrödinger equation. However, my main mistake was not shifting the kinetic energy to the center. The command "ResourceFunction["FourierShift"]" resolves this issue. Unfortunately, I only realized this some time after submitting the question. Using Fourier to get the numerical FT of a Gaussian Attachments:
Posted 24 days ago
 You're using the discrete Fourier transform and the whole split-operator method thing to solve the two-dimensional Schrödinger equation. In the first code block, the wave packet with momentum primarily in the y-direction, is splitting into two parts after a few time steps. Is this physically correct for a free particle? I think the error is not an intrinsic issue of splitting wave packets, but the potential issue lies with numerical resolution or accuracy. TDFourier2D[timesteps_, potheight_, ellip_] := Module[{NP, L, dx, dt, kx, ky, x0, y0, xmax, xmin, ymax, ymin, kxmax, kymax, \[Rho], TimeEnd, psi, K3, nt, UK, UV, kgridx, kgridy}, NP = 100; L = 1; dx = L/(NP - 1); dt = (1/20) dx^2; TimeEnd = timesteps; \[Rho] = 0.16; x0 = 0.0; y0 = 0.0; kx = 0; ky = 50; xmax = 1; xmin = -1; ymax = 1; ymin = -1; kgridx = Array[# &, {NP}, {-Pi/dx, Pi/dx}]; kgridy = Array[# &, {NP}, {-Pi/dx, Pi/dx}]; psi3[0] = Array[(1/(\[Rho] Sqrt[ Pi])) Exp[-I kx #1 - (#1 - x0)^2/(2 \[Rho]^2)]* Exp[-I ky #2 - (#2 - y0)^2/(2 \[Rho]^2)] &, {NP, NP}, {{xmin, xmax}, {ymin, ymax}}] // Quiet; K3 = 1/2 (kgridx^2 + kgridy^2); V3 = Array[(potheight/2) #1^2 + (ellip*potheight/2) #2^2 &, {NP, NP}, {{xmin, xmax}, {ymin, ymax}}]; UV[dt_] := Exp[-I V3 dt]; UK[dt_] := Exp[-I K3 dt]; U[dt_]@psi_ := UV[dt]*InverseFourier[UK[dt]*Fourier[psi]]; Do[psi3[nt] = U[dt][psi3[nt - 1]], {nt, 1, TimeEnd}]; psi3] ellipticalform = 0.2; TimeEnd = 160; potheight = 0; {timing, psi3Result} = AbsoluteTiming[TDFourier2D[TimeEnd, potheight, ellipticalform]]; graphTable1 = GraphicsGrid[ Partition[ Table[Labeled[ ListDensityPlot[Abs[psi3Result[i]], Axes -> True, PlotRange -> All, ColorFunction -> "SunsetColors", PerformanceGoal -> "Quality"], {"x", "y"}, {Left, Top}, RotateLabel -> True], {i, 1, TimeEnd, TimeEnd/4}], 2], Frame -> All, ImageSize -> 400] Now if I can just get my mouse working, I think that the importance of grid resolution and time step size is a critical way to accurately simulate the Schrödinger equation. It turns out that using the split-operator method, it's possible to increase the grid re-solution which therefore minimizes the effect, of aliasing in the Fourier transform. This might seem like something we could do by zooming in on the grid but no, we've got to clear out what we've got programmatically which will allow us to try a smaller time step. But you can try all the smaller time steps and when you do, you will find that the wave packet's evolution is accurately computed at each step, reducing the computational errors that could lead to the wave packet's artificial splitting. I imagine that in future versions of the Wolfram Language these Pandas-like artifacts such as splitting wherein a large dt might not capture the dynamics accurately, will never have a smaller dt that reduces these inaccuracies. In the following block of code you will discover that a smaller time step in fact does lead to more accurate time evolution, especially in this newfangled split-operator method, which relies on operator splitting over small intervals. TDFourier2D[timesteps_, potheight_, ellip_] := Module[{NP, L, dx, dt, kx, ky, x0, y0, xmax, xmin, ymax, ymin, kxmax, kymax, \[Rho], TimeEnd, psi, K3, nt, UK, UV, kgridx, kgridy}, NP = 256; L = 1; dx = L/(NP - 1); dt = (1/50) dx^2; TimeEnd = timesteps; \[Rho] = 0.16; x0 = 0.0; y0 = 0.0; kx = 0; ky = 50; xmax = 1; xmin = -1; ymax = 1; ymin = -1; kgridx = Array[# &, {NP}, {-Pi/dx, Pi/dx}]; kgridy = Array[# &, {NP}, {-Pi/dx, Pi/dx}]; psi3[0] = Array[(1/(\[Rho] Sqrt[ Pi])) Exp[-I kx #1 - (#1 - x0)^2/(2 \[Rho]^2)]* Exp[-I ky #2 - (#2 - y0)^2/(2 \[Rho]^2)] &, {NP, NP}, {{xmin, xmax}, {ymin, ymax}}] // Quiet; K3 = 1/2 (kgridx^2 + kgridy^2); V3 = Array[(potheight/2) #1^2 + (ellip*potheight/2) #2^2 &, {NP, NP}, {{xmin, xmax}, {ymin, ymax}}]; UV[dt_] := Exp[-I V3 dt]; UK[dt_] := Exp[-I K3 dt]; U[dt_]@psi_ := UV[dt]*InverseFourier[UK[dt]*Fourier[psi]]; Do[psi3[nt] = U[dt][psi3[nt - 1]], {nt, 1, TimeEnd}]; psi3] ellipticalform = 0.2; TimeEnd = 160; potheight = 0; {timing, psi3Result} = AbsoluteTiming[TDFourier2D[TimeEnd, potheight, ellipticalform]]; graphTable1 = GraphicsGrid[ Partition[ Table[Labeled[ ListDensityPlot[Abs[psi3Result[i]], Axes -> True, PlotRange -> All, ColorFunction -> "SunsetColors", PerformanceGoal -> "Quality"], {"x", "y"}, {Left, Top}, RotateLabel -> True], {i, 1, TimeEnd, TimeEnd/4}], 2], Frame -> All, ImageSize -> 400] That is how we make sure that in this third iteration we have increased the grid resolution "NP". A higher grid resolution, means we use more points to represent the wave function, which improves the accuracy of the discrete Fourier transform and so on and so forth. This helps in reducing numerical artifacts such as aliasing, which is an implicit unintended splitting of the wave packet. And that is why yes, the wave packet with momentum primarily in the y-direction is splitting into two parts after a few time steps, which isn't physically correct for a free particle. This begs the question of whether we can articulate how we achieve all our..wave packets with regard to sufficient resolution to minimize aliasing effects and gauging smaller time steps such that the wave packet's evolution is continuously and accurately computed at each step.