Chciałbym uruchomić kilka prostych symulacji rozpraszania pakietów fal z prostych potencjałów w jednym wymiarze.

Czy istnieją proste sposoby numerycznego rozwiązania jednowymiarowego TDSE dla pojedynczej cząstki? Wiem, że ogólnie rzecz biorąc, próba zastosowania naiwnych podejść do zintegrowania równań różniczkowych cząstkowych może szybko zakończyć się katastrofą. Dlatego szukam algorytmów, które

• są stabilne numerycznie,
• są łatwe do wdrożenia lub mają łatwo dostępne implementacje biblioteki kodów,
• biegnij dość szybko i miejmy nadzieję
• są stosunkowo proste do zrozumienia.

Chciałbym również omijać metody spektralne, a zwłaszcza metody, które są niczym więcej niż rozwiązywaniem równania Schrödingera niezależnego od czasu jak zwykle. Byłbym jednak zainteresowany metodami pseudo-spektralnymi, które wykorzystują splajny B lub tym podobne. Jeśli metoda może wymagać potencjału zależnego od czasu, jest to zdecydowanie bonus.

Oczywiście, każda taka metoda zawsze będzie miała szereg wad, więc chciałbym o nich usłyszeć. Kiedy to nie działa? Jakie są typowe pułapki? Które sposoby mogą być popchnięte, a które nie?

Równanie Schroedingera jest faktycznie równaniem dyfuzji reakcyjnej (wszystkie stałe to 1). Jeśli chodzi o równanie różniczkowe cząstkowe, istnieją dwa sposoby ich rozwiązania:

$$i\frac{\partial\psi}{\partial t}=-\nabla^2\psi+V\psi\tag{1}$$

1. Metoda niejawna (adv: duże kroki czasowe i bezwarunkowo stabilna, disadv: wymaga rozwiązania macierzy, które może dawać złe dane)
2. Metoda jawna (adv: łatwa do wdrożenia, disadv: wymaga małych kroków czasowych dla stabilności)

W przypadku równań parabolicznych (liniowych i 2 rzędu w x ) metoda niejawna jest często lepszym wyborem. Powodem jest warunek stabilności dla metody jawnej wymaga d t ∝ d x 2 , która będzie bardzo mała. Możesz uniknąć tego problemu, stosując metodę niejawną, która nie ma takich ograniczeń co do kroku czasowego (chociaż w praktyce zwykle nie czynisz go niesamowicie dużym, ponieważ możesz stracić część fizyki). Następnie opiszę metodę Crank-Nicolsona , wspólny niejawny schemat dokładności drugiego rzędu (czasoprzestrzeni).$t$$x$$dt\propto dx^2$

Przystawki

Aby rozwiązać obliczeniowo PDE, musisz go zdyskretyzować (dopasować zmienne do siatki). Najbardziej prosta jest prostokątna, kartezjańska siatka. Tutaj oznacza indeks czasu (i zawsze jest superskryptem), a j indeks pozycji (zawsze jest indeksem dolnym). Dzięki zastosowaniu rozszerzenia Taylora dla zmiennej zależnej od pozycji równanie (1) staje się i ψ n + 1 j - ψ $n$$j$ Gdzie przyjęliśmy, żeV=V(x). Następnie następuje grupowanie podobnych wskaźników przestrzennych i czasowych (możesz dwukrotnie sprawdzić matematykę):

$$i\frac{\psi^{n+1}_j-\psi_j}{dt} = -\frac12\left(\frac{\psi_{j+1}^{n+1}-2\psi_j^{n+1}+\psi_{j-1}^{n+1}}{dx^2}+\frac{\psi_{j+1}^{n}-2\psi_j^{n}+\psi_{j-1}^{n}}{dx^2}\right) \\ +\frac12\left(V_j\psi_j^{n+1} +V_j\psi_j^n\right)$$ $V=V(x)$ To równanie ma postać (A0A-00⋯A+A0A-0⋯0A+A0A-⋯⋮⋮⋮⋱⋮)(ψ n + 1 0 ψ n + 1 1 ⋮ψ n + 1 J - 1 )= $$\frac12\frac{dt}{dx^2}\psi_{j+1}^{n+1}+\left(i-\frac{dt}{dx^2}-\frac12V_j\right)\psi_j^{n+1}+\frac12\frac{dt}{dx^2}\psi_{j-1}^{n+1}=\\ i\psi_j^n-\frac12\frac{dt}{dx^2}\left(\psi_{j+1}^n-2\psi_j^n+\psi_{j-1}^n\right)+\frac12V_j\psi_j^n\tag{2}$$ co jest zwane matryca trój przekątną i maznanego rozwiązania! (oraz przykładów roboczych, w tym piśmie przeze mnie). Sposób ten wyraźne zadrapania spośród całej lewej strony (a raczej górną granicę?) W równaniu (2), z wyjątkiemı* F n + 1, j okresie. $$\left(\begin{array}{ccccc}A_0 & A_- & 0 & 0 &\cdots\\ A_+ & A_0 & A_- & 0 &\cdots\\ 0 & A_+ & A_0 & A_- &\cdots \\ \vdots & \vdots & \vdots & \ddots & \vdots\end{array}\right)\left(\begin{array}{c}\psi_0^{n+1} \\ \psi_1^{n+1}\\ \vdots \\ \psi_{J-1}^{n+1}\end{array}\right)=\left(\begin{array}{c}\psi_0^{n} \\ \psi_1^{n}\\ \vdots \\ \psi_{J-1}^{n}\end{array}\right)$$ $i\psi_j^{n+1}$

Zagadnienia

Największym problemem, jaki znalazłem przy zastosowaniu metod niejawnych, jest to, że są one silnie zależne od warunków brzegowych. Jeśli masz źle zdefiniowane / zaimplementowane warunki brzegowe, możesz uzyskać fałszywe oscylacje w swoich komórkach, które mogą prowadzić do złych wyników (zobacz mój post SciComp na podobny temat). Prowadzi to do uzyskania dokładności pierwszego rzędu w przestrzeni, a nie drugiej, którą powinien dać twój plan .

Metody niejawne są również rzekomo trudne do zrównoleglenia, ale użyłem ich tylko do równań ciepła 1D i nie potrzebowałem wsparcia równoległego, więc nie mogę zweryfikować ani odrzucić twierdzenia.

Nie jestem również pewien, w jaki sposób złożony charakter funkcji falowej wpłynie na obliczenia. Praca, którą wykonałem, wykorzystuje równania dynamiczne płynu Eulera , a zatem są całkowicie rzeczywiste o wielkościach nieujemnych.

Potencjał zależny od czasu

Jeśli masz analityczny potencjał zależny od czasu (np ), to by po prostu użyć aktualnego czasu, t , dla V j na RHS (2) oraz czasu przyszłego, t + d t , na LHS. Nie wierzę, że spowodowałoby to jakiekolwiek problemy, ale nie przetestowałem tego, więc nie mogę również zweryfikować ani zaprzeczyć temu aspektowi.$V\propto \cos(\omega t)$$t$$V_j$$t+dt$

Alternatywy

Istnieje również kilka interesujących alternatyw dla metody Crank-Nicolson. Pierwszym krokiem jest tak zwana metoda „przeskakiwania w czasie”. W tej jawnej metodzie bierzesz krok czasowy ( ) i używasz pierwiastków wielomianów Czebyszewa, aby uzyskać zoptymalizowany zestaw kroków czasowych, które szybko sumują się do d t szybciej niż robienie d t / N kroków N razy (faktycznie dostajesz Δ T = N 2 d t, tak że każdy krok N posuwa cię naprzód N d $dt\propto dx^2$$dt$$dt/N$$N$$\Delta T=N^2dt$$N$$Ndt$o czasie). (Używam tej metody w swoich badaniach, ponieważ masz dobrze zdefiniowany „strumień” z jednej komórki do drugiej, który jest używany do scalania danych z jednego procesora do drugiego, używając schematu Crank-Nicolson, którego nie byłem w stanie zrobić).

EDYCJA Należy zwrócić uwagę na to, że metoda ta jest dokładna z dokładnością do pierwszego rzędu w czasie, ale jeśli użyjesz metody Runge-Kutta 2 w połączeniu, da ci to dokładny schemat drugiego rzędu w czasie.

Drugi nazywany jest jawnym kierunkiem przemiennym . Ta metoda wymaga znanych i dobrze zdefiniowanych warunków brzegowych. Następnie przystępuje do rozwiązania równania przy użyciu granicy bezpośrednio w obliczeniach (nie trzeba jej stosować po każdym kroku). W tej metodzie rozwiązujesz PDE dwa razy, raz w górę i raz w dół. Przeciąganie w górę używa podczas gdy wobulacji w dół używa ∂2

$$\frac{\partial^2\psi}{\partial x^2}\approx\frac{\psi_{j-1}^{n+1}-\psi_j^{n+1}-\psi_j^n+\psi_{j+1}^n}{dx^2}$$ dla równania dyfuzji, podczas gdy pozostałe warunki pozostaną takie same. Krok czasowyn+1jest następnie rozwiązany przez uśrednienie dwóch kierunkowych przebiegów. $$\frac{\partial^2\psi}{\partial x^2}\approx\frac{\psi_{j+1}^{n+1}-\psi_j^{n+1}-\psi_j^n+\psi_{j-1}^n}{dx^2}$$ $n+1$

Na początku lat 90. szukaliśmy metody rozwiązania TDSE wystarczająco szybko, aby wykonywać animacje w czasie rzeczywistym na PC i natknęliśmy się na zaskakująco prostą, stabilną, wyraźną metodę opisaną przez PB Visschera w Computers in Physics : „ Szybki, wyraźny algorytm dla zależnego od czasu równania Schrödingera ”. Visscher zauważa, że ​​jeśli podzielisz funkcję fali na rzeczywistą i wymyśloną część, , SE stanie się systemem:$\psi=R+iI$

$$\begin{eqnarray}\frac{dR}{dt}&=&HI \\ \frac{dI}{dt}&=&-HR \\ H&=&-\frac{1}{2m}\nabla^2+V\end{eqnarray}$$

Jeśli następnie obliczyć i że czasami przesuniętych ( R na 0 , hemibursztynianu t , 2 hemibursztynianu t , . . . I I w 0,5 hemibursztynianu t , 1,5 Δ t $R$$I$$R$$0,\Delta t,2\Delta t,...$$I$ , Można uzyskać dyskretyzacji:$0.5\Delta t, 1.5\Delta t,...)$

$$R(t+\frac{1}{2} \Delta t)=R(t-\frac{1}{2} \Delta t)+\Delta t HI(t)$$

$$I(t+\frac{1}{2} \Delta t)=I(t-\frac{1}{2} \Delta t)-\Delta t HR(t)$$

$$\nabla^2\psi(r,t)=\frac{\psi(r+\Delta r,t)-2\psi(r,t)+\psi(r-\Delta r,t)}{\Delta r^2}$$

$\Delta t$

Definiowanie gęstości prawdopodobieństwa jako

$$P(x,t)=R^2(x,t)+I(x,t+\frac{1}{2} \Delta t)I(x,t-\frac{1}{2} \Delta t)$$

$$P(x,t)=R(x,t+\frac{1}{2} \Delta t)R(x,t-\frac{1}{2} \Delta t)+I^2(x,t)$$

sprawia, że ​​algorytm jest jednolity, co pozwala zachować prawdopodobieństwo.

Przy wystarczającej optymalizacji kodu byliśmy w stanie uzyskać bardzo ładne animacje obliczane w czasie rzeczywistym na maszynach 80486. Uczniowie mogli „narysować” dowolny potencjał, wybrać całkowitą energię i obserwować ewolucję w czasie pakietu gaussowskiego.

$x$$t$

To, co będziesz robić, to ukryta wersja Metody propagacji wiązki do optycznej propagacji przez falowód o zmiennym przekroju (analogicznie do potencjałów zmieniających się w czasie), więc warto byłoby to również sprawdzić.

Sposób, w jaki patrzę na SSFM / BPM jest następujący. Jego podstawą jest formuła produktu Trottera z teorii Liego:

$$\tag{1}\lim\limits_{m\to\infty}\left(\exp\left(\mathcal{D}\,\frac{t}{m}\right)\,\exp\left(\mathcal{V}\,\frac{t}{m}\right)\right)^m = \exp((\mathcal{D+V}) t)$$

$x-y$$x-y-z$$\psi(x,y,z)$$t$$N$$\Psi$ (for a $1024\times1024$ grid we have $N=1024^2=1\,048\,576$) and then your Schrödinger equation is of the form:

$$\tag{2}\mathrm{d}_t \Psi = K\Psi = (\mathcal{D+V}(t)) \Psi$$

where $K = \mathcal{D+V}$ is an $N\times N$ skew-Hermitian matrix, an element of $\mathfrak{u}(N)$, and $\Psi$ is going to be mapped with increasing time by an element of the one parameter group $\exp(K\,t)$. (I've sucked the $i\hbar$ factor into the $K = \mathcal{D+V}$ on the RHS so I can more readily talk in Lie theoretic terms). Given the size of $N$, the operators' natural habitat $\mathfrak{U}(N)$ is a thoroughly colossal Lie group so PHEW! yes I am still talking in wholly theoretical terms!. Now, what does $\mathcal{D+V}$ look like? Still imagining for now, it could be thought of as a finite difference version of $i\,\hbar\,\nabla^2/(2\,m) - i\hbar^{-1}V_0 + i\hbar^{-1}(V_0-V(x,y,z,t_0))$, where $V_0$ is some convenient "mean" potential for the problem at hand.

We let:

$$\tag{3}\begin{array}{lcl}\mathcal{D} &=& i\frac{\hbar}{2\,m} \nabla^2 - i\hbar^{-1}V_0\\ \mathcal{V}&=&i\hbar^{-1}(V_0-V(x,y,z,t))\end{array}$$

Why I have split them up like this will become clear below.

The point about $\mathcal{D}$ is that it can be worked out analytically for a plane wave: it is a simple multiplication operator in momentum co-ordinates. So, to work out $\Psi\mapsto\exp(\Delta t\,\mathcal{D}) \Psi$, here are the first three steps of a SSFM/BPM cycle:

1. Impart FFT to dataset $\Psi$ to transform it into a set $\tilde{\Psi}$ of superposition weights of plane waves: now the grid co-ordinates have been changed from $x,\,y,\,z$ to $k_x,\,k_y,\,k_z$;
2. Impart $\tilde{\Psi}\mapsto\exp(\Delta t\,\mathcal{D}) \tilde{\Psi}$ by simply multiplying each point on the grid by $\exp(i\,\Delta t (V_0-k_x^2+k_y^2+k_z^2)/\hbar)$;

3. Impart inverse FFT to map our grid back to $\exp(\Delta t\,\mathcal{D}) \Psi$

   .Now we're back in position domain. This is the better domain to impart the operator $\mathcal{V}$ of course: here $\mathcal{V}$ is a simple multiplication operator. So here is your last step of your algorithmic cycle:

4. Impart the operator $\Psi\mapsto\exp(\Delta t\,\mathcal{V}) \Psi$ by simply multiplying each point on the grid by the phase factor $\exp(i\,\Delta t\,(V_0-V(x,y,z,t))/\hbar)$

....and then you begin your next $\Delta t$ step and cycle over and over. Clearly it is very easy to put time-varying potentials $V(x,y,z,t)$ into the code.

So you see you simply choose $\Delta t$ small enough that the Trotter formula (1) kicks in: you're simply approximating the action of the operator $\exp(\mathcal{D+V}\,\Delta t)\approx\exp(\mathcal{D}\,\Delta t)\,\exp(\mathcal{V}\,\Delta t)$ and you flit back and forth with your FFT between position and momentum co-ordinates, i.e. the domains where $\mathcal{V}$ and $\mathcal{D}$ are simple multiplication operators.

Notice that you are only ever imparting, even in the discretised world, unitary operators: FFTs and pure phase factors.

One point you do need to be careful of is that as your $\Delta t$ becomes small, you must make sure that the spatial grid spacing shrinks as well. Otherwise, suppose the spatial grid spacing is $\Delta x$. Then the physical meaning of the one discrete step is that the diffraction effects are travelling at a velocity $\Delta x/\Delta t$; when simulating Maxwell's equations and waveguides, you need to make sure that this velocity is much smaller than $c$. I daresay like limits apply to the Schrödinger equation: I don't have direct experience here but it does sound fun and maybe you could post your results sometime!

A second "experience" point with this kind of thing - I'd be almost willing to bet this is how you'll wind up following your ideas. We often have ideas that we want to do simple and quick and dirty simulations but it never quite works out that way! I'd begin with the SSFM as I've described above as it is very easy to get running and you'll quickly see whether or not its results are physical. Later on you can use your, say Mathematica SSFM code check the results of more sophisticated code you might end up building, say, a Crank Nicolson code along the lines of Kyle Kanos's answer.

---

Error Bounds

The Dynkin formula realisation of the Baker-Campbell-Hausdorff Theorem:

$$\exp(\mathcal{D}\Delta t)\,\exp(\mathcal{V})\Delta t) = \exp\left((\mathcal{D}+\mathcal{V})\Delta t + \frac{1}{2} [\mathcal{D},\,\mathcal{V}]\,\Delta t^2 + \cdots\right)$$ converging for some $\Delta t>0$ shows that the method is accurate to second order and can show that:

$$\exp(\mathcal{D}\Delta t)\,\exp(\mathcal{V})\Delta t)\,\exp\left(-\frac{1}{2} [\mathcal{D},\,\mathcal{V}]\,\Delta t^2\right) = \exp\left((\mathcal{D}+\mathcal{V})\Delta t + \mathcal{O}(\Delta t^3)\right)$$

You can, in theory, therefore use the term $\exp(\mathcal{V})\Delta t)\,\exp\left(-\frac{1}{2} [\mathcal{D},\,\mathcal{V}]\,\Delta t^2\right)$ to estimate the error and set your $\Delta t$ accordingly. This is not as easy as it looks and in practice bounds end up being instead rough estimates of the error. The problem is that:

$$\frac{\Delta t^2}{2}[\mathcal{D},\,\mathcal{V}] = -\frac{i\,\Delta t^2}{2\,m}\,\left(\partial_x^2 V(x,\,t) + 2 \partial_x V(x,\,t)\,\partial_x\right)$$

and there are no readily transformed to co-ordinates wherein $[\mathcal{D},\,\mathcal{V}]$ is a simple multiplication operator. So you have to be content with $\exp\left(-\frac{1}{2} [\mathcal{D},\,\mathcal{V}]\,\Delta t^2\right) \approx e^{-i\,\varphi\,\Delta t^2}\left(\mathrm{id} -\left(\frac{1}{2} [\mathcal{D},\,\mathcal{V}]\,-i\,\varphi(t)\right)\,\Delta t^2\right)$ and use this to estimate your error, by working out $\left(\mathrm{id} -\left(\frac{1}{2} [\mathcal{D},\,\mathcal{V}]\,-i\,\varphi(t)\right)\,\Delta t^2\right) \,\psi$ for your currently evolving solution $\psi(x,\,t)$ and using this to set your $\Delta t$ on-the-fly after each cycle of the algorithm. You can of course make these ideas the basis for an adaptive stepsize controller for your simulation. Here $\varphi$ is a global phase pulled out of the dataset to minimise the norm of $\left(\frac{1}{2} [\mathcal{D},\,\mathcal{V}]\,-i\,\varphi(t)\right)\,\Delta t^2$; you can of course often throw such a global phase out: depending on what you're doing with the simulation results often we're not bothered by a constant phase global $\exp\left(\int \varphi\,\mathrm{d}t\right)$.

A relevant paper about errors in the SSFM/BPM is:

Lars Thylén. "The Beam Propagation Method: An Analysis of its Applicability", Optical and Quantum Electronics 15 (1983) pp433-439.

Lars Thylén thinks about the errors in non-Lie theoretic terms (Lie groups are my bent, so I like to look for interpretations of them) but his ideas are essentially the same as the above.

I can recommend using the finite-difference time-domain (FDTD) method. I even wrote a tutorial some time back that should answer most of your questions:

J. R. Nagel, "A review and application of the finite-difference time-domain algorithm applied to the Schrödinger equation," ACES Journal, Vol. 24, No. 1, February 2009

I have some Matlab codes that run nicely for 1D systems. If you have experience with FDTD doing electromagnetics, it works great for quantum mechanics as well. I can post my codes if you're interested.

Basically, it just operates on the wavefunctions directly by splitting the derivatives up into finite differences. It is kind of similar to the Crank-Nicholson scheme, but not exactly. If you are familiar with FDTD from electromagnetic wave theory, then FDTD will be very intuitive when solving the Schrodinger equation.

The most straightforward finite difference method is fast and easy to understand, but is not unitary in time - so probability is not conserved. Crank-Nicholson-Crout averages the forward and backward finite difference methods to produce a hybrid implicit/explicit method that is still pretty easy to understand and to implement and is unitary in time. This site explains the method well, provides pseudocode, and gives the relevant properties:

http://www.physics.utah.edu/~detar/phycs6730/handouts/crank_nicholson/crank_nicholson/ Note: There is a - sign missing from the LHS of equation one of this link, which propagates throughout the page.

Where does the nonunitarity come from?

In a nut shell, solving the TDSE comes down to figuring out how to deal with

$| \psi(x,t)\rangle = e^{-iHt}|\psi(x,0)\rangle$

which contains a differential operator in an exponential.

Applying a forward finite difference turns the differential operator into a tridiagonal matrix (converting the Reals to a grid) and the exponential into the first two terms of its Taylor series

$e^{-iHt}\approx1-iHt$

This discretization and linearization is what gives rise to the nonunitarity. (You can show that the tridiagonal matrix is not unitary by direct computation.) Combining the forward finite difference with the backward finite difference produces the approximation

$e^{-iHt}\approx \frac {1-\frac{1}{2} iHt} {1+\frac{1}{2} iHt}$

which, kindly, happens to be unitary (again you can show it by direct computation).

A few answers and comments here conflate confusingly the TDSE with a wave equation; perhaps a semantics issue, to some extent. The TDSE is the quantized version of the classical non-relativistic hamiltonian

$$H=\frac{p^2}{2m} + V(x)= E.$$ With the rules $$p \rightarrow i\hbar\partial_x,\ \ E\rightarrow i\hbar\partial_t, \ \ x\rightarrow x,$$ (as discussed in chapter 1 of d'Espagnat, Conceptual foundations of quantum mechanics, https://philpapers.org/rec/ESPCFO), it therefore reads $$\left[-\frac{\hbar^2}{2m}\partial_{xx} + V(x)\right]\psi = i\hbar\partial_t\psi,$$ so it is clearly a diffusion-like equation. If one used the relativistic energy, which contains a E$^2$ term, then a wave-like equation such as $$\partial_{xx}\psi = \partial_{tt}\psi +\dots$$ would obtain (for V=0 for simplicity), such as the Pauli or Klein-Gordon equations. But that is, of course, a completely different matter.

Now, back to the TDSE, the obvious method is Crank-Nicolson as has been mentioned, because it is a small-time expansion that conserves the unitarity of the evolution (FTCS, e.g., does not). For the 1-space-D case, it can be treated as a matrix iteration, reading

$$\vec{\psi}^{n+1}=({\cal I} + \frac{i\tau}{2\hbar}\tilde{H})^{-1} ({\cal I} - \frac{i\tau}{2\hbar}\tilde{H})\vec{\psi}^{n}$$ with ${\cal I}$ the identity matrix and $$H_{jk}=(\tilde{H})_{jk}=\frac{-\hbar^{2}}{2m}\left[\frac{\delta_{j+1,k}+\delta_{j-1,k}-2\delta_{jk}}{h^{2}}\right]+V_{j}\delta_{jk}.$$ (Details e.g. in Numerical methods for physics, http://algarcia.org/nummeth/nummeth.html, by A. L. Garcia). As most clearly seen in periodic boundary conditions, a space-localized $\psi$ spreads out in time: this is expected, because the initial localized $\psi$ is not an eigenstate of the stationary Schroedinger equation, but a superposition thereof. The (classical massive free particle) eigenstate with fixed momentum (for the non-relativistic kinetic operator) is simply $\psi_s=e^{ikx}/\sqrt{L}$, i.e. fully delocalized as per Heisenberg principle, with constant probability density 1/L everywhere (note that I am avoiding normalization issues with continuum states by having my particle live on a finite, periodically repeated line). Using C-N, the norm $$\int |\psi|^2 dx$$ is conserved, thanks to unitarity (this is not the case in other schemes, such as FTCS, e.g.). Incidentally, notice that starting from an energy espression such as $$cp=E$$ with $c$ fixed, you'd get $$ic\hbar\partial_x=i\hbar\partial_t$$ i.e. the advection equation, which has no dispersion (if integrated properly with Lax-Wendroff methods), and your wavepacket will not spread in time in that case. The quantum analog is the massless-particle Dirac equation.