- #1
Loro
- 80
- 1
I'm trying to simulate wavepacket propagation in various potentials (so far in 1D). I'm writing in C++, and using allegro graphics library.
The TDSE is:
[itex] \frac{∂ψ}{∂t} = \frac{i \hbar}{2m} \frac{∂^2 ψ}{∂ x^2} - \frac{iV}{\hbar}ψ[/itex]
So what I've done was writing out the real and imaginary parts:
[itex] \frac{∂ψ_r}{∂t} = -\frac{i \hbar}{2m} \frac{∂^2 ψ_i}{∂ x^2} + \frac{iV}{\hbar}ψ_i[/itex]
[itex] \frac{∂ψ_i}{∂t} = +\frac{i \hbar}{2m} \frac{∂^2 ψ_r}{∂ x^2} - \frac{iV}{\hbar}ψ_r[/itex]
And then I changed the spatial derivatives into ratios of differences in a usual way.
To calculate the values of ψ's at the next time step, I use:
[itex] ψ_j(t + Δt) = ψ_j(t) + \frac{∂ψ_j (t) }{∂t}Δt + \frac{1}{2}\frac{∂^2ψ_j (t) }{∂t^2}Δt^2 [/itex]
(j is r or i)
Where for expressions for the time derivatives, I plug in the approximated right-hand sides of the TDSE.
So far I have no boundary conditions - I just don't process the boundary points. I've tried setting them to zero; also I thought of making them periodic but I'm not quite sure how to do it yet.
So the program simply calculates the ψ's as described above. Scales them to make sure normalization is ok (works perfectly). Plots the [itex] |ψ|^2 [/itex]. And starts over.
Problems:
There must be something wrong either with my equations, or/and with my numerical method, or/and with my boundary conditions, because:
1. The propagation without any potential works just fine. But if I put in constant potentials (V=const) - I get different spreading rates, depending on the value of V.
2. If I put in e.g. a harmonic oscillator potential - same thing happens regardless of whether I put in [itex]-x^2[/itex] or [itex]+x^2[/itex]. The wavepacket splits in 2 (with different ratios depending on the initial displacement) and they both travel outwards, more and more slowly, but never stop (changing the spring constant doesn't do much, unless we make it extremely small). Should the splitting even occur? In particular if we start in the equilibrium position, with no momentum - it splits into 2 identical pieces, moving apart - and I would've thought that it would rather stay there...?
3. The boundaries start talking to the wavepackets after some time, and I really don't know what to do about it. E.g. when I plug in a step potential (0 on the left, and constant on the right of the green y-axis) this happens:
The initial shape is Gaussian, propagating to the right:
https://dl.dropbox.com/u/94695102/fizyka/schr01.jpg
After hitting the step, it splits into reflected and transmitted parts, as expected, but we already start seeing something growing on the sides:
https://dl.dropbox.com/u/94695102/fizyka/schr02.jpg
A few seconds later, suddenly a bump grows out on the right, and kills the wavepackets:
https://dl.dropbox.com/u/94695102/fizyka/schr03.jpg
4. Sometimes (when V≠0) a propagating wavepacket leaves some junk in the place where it used to be, and after some time this junk suddenly grows, and kills the wavepacket.
If any of you have any suggestions about what I could try, or if you could refer me to a paper which adresses these issues, I would be very happy.
I've used Lax-Friedrich scheme already, and it solves some of the issues, but not most of them, and also introduces extra spreading.
I know I could try Crank-Nicholson method, but it would be a serious revolution in the code, and I suspect the error lies somewhere else.
The TDSE is:
[itex] \frac{∂ψ}{∂t} = \frac{i \hbar}{2m} \frac{∂^2 ψ}{∂ x^2} - \frac{iV}{\hbar}ψ[/itex]
So what I've done was writing out the real and imaginary parts:
[itex] \frac{∂ψ_r}{∂t} = -\frac{i \hbar}{2m} \frac{∂^2 ψ_i}{∂ x^2} + \frac{iV}{\hbar}ψ_i[/itex]
[itex] \frac{∂ψ_i}{∂t} = +\frac{i \hbar}{2m} \frac{∂^2 ψ_r}{∂ x^2} - \frac{iV}{\hbar}ψ_r[/itex]
And then I changed the spatial derivatives into ratios of differences in a usual way.
To calculate the values of ψ's at the next time step, I use:
[itex] ψ_j(t + Δt) = ψ_j(t) + \frac{∂ψ_j (t) }{∂t}Δt + \frac{1}{2}\frac{∂^2ψ_j (t) }{∂t^2}Δt^2 [/itex]
(j is r or i)
Where for expressions for the time derivatives, I plug in the approximated right-hand sides of the TDSE.
So far I have no boundary conditions - I just don't process the boundary points. I've tried setting them to zero; also I thought of making them periodic but I'm not quite sure how to do it yet.
So the program simply calculates the ψ's as described above. Scales them to make sure normalization is ok (works perfectly). Plots the [itex] |ψ|^2 [/itex]. And starts over.
Problems:
There must be something wrong either with my equations, or/and with my numerical method, or/and with my boundary conditions, because:
1. The propagation without any potential works just fine. But if I put in constant potentials (V=const) - I get different spreading rates, depending on the value of V.
2. If I put in e.g. a harmonic oscillator potential - same thing happens regardless of whether I put in [itex]-x^2[/itex] or [itex]+x^2[/itex]. The wavepacket splits in 2 (with different ratios depending on the initial displacement) and they both travel outwards, more and more slowly, but never stop (changing the spring constant doesn't do much, unless we make it extremely small). Should the splitting even occur? In particular if we start in the equilibrium position, with no momentum - it splits into 2 identical pieces, moving apart - and I would've thought that it would rather stay there...?
3. The boundaries start talking to the wavepackets after some time, and I really don't know what to do about it. E.g. when I plug in a step potential (0 on the left, and constant on the right of the green y-axis) this happens:
The initial shape is Gaussian, propagating to the right:
https://dl.dropbox.com/u/94695102/fizyka/schr01.jpg
After hitting the step, it splits into reflected and transmitted parts, as expected, but we already start seeing something growing on the sides:
https://dl.dropbox.com/u/94695102/fizyka/schr02.jpg
A few seconds later, suddenly a bump grows out on the right, and kills the wavepackets:
https://dl.dropbox.com/u/94695102/fizyka/schr03.jpg
4. Sometimes (when V≠0) a propagating wavepacket leaves some junk in the place where it used to be, and after some time this junk suddenly grows, and kills the wavepacket.
If any of you have any suggestions about what I could try, or if you could refer me to a paper which adresses these issues, I would be very happy.
I've used Lax-Friedrich scheme already, and it solves some of the issues, but not most of them, and also introduces extra spreading.
I know I could try Crank-Nicholson method, but it would be a serious revolution in the code, and I suspect the error lies somewhere else.
Last edited by a moderator: