How to turn model of Schrödinger's Equation 2D?

[IsaiahvH]

Hi,

I am a student in the Netherlands, currently 17 years old and at the end of my 'middelbare school', meaning that next year I'll be a bachelor student at a university.

I am doing an extended essay/research thing that is custom you do in your last year, with a friend of mine.
We picked the topic of antimatter and we want to look at the energy values of electrons.
Our teacher suggested we would create a model based on Schrödinger's Equation in MATLAB.

I have created such a model, 1D, but I don't fully understand the equation, which seems rather abstract to me, and have no idea how to create a 2D version and eventually a 3D version. Our teacher's advice? Ask you guys (or girls, of course).

My first question: What does Schrödinger's equation mean?
This is the least important of my questions; I understand what it provides (particle probability) and how to construct a model out of it, but it feels wrong to work with something I don't fully understand. Examples found on the internet are difficult to comprehend for me, but nevertheless I'd hoped someone could shed some light on this.

The formula we use, is this one.

It's a somewhat rewritten version from a schoolbook.
This has been turned into the following code.

for x = 0:dx:max
    fu(z) = C * x.^2;
    curve = -1 * massfactor * (E-fu(z)) * wave;
    slope = slope + curv * dx;
    d_wave = slope * dx;
    wave = wave + d_wave;

    fy(z) = wave;
    z = z + 1;
end

Spoiler: Variable explanation

x is the distance in the set dimension (1D) of the situation.
C is not explained in our book
Curve refers to the current curve (our book is not clearer)
Massfactor is the 8*pi^2*m/h^2 part, bundled together
Slope is the current slope of the graph
Dx is the delta x
D_wave is the delta wave
Wave is the square of the probability function, also Ψ, if I'm correct

Of course, all variables are given a start value beforehand, but basically this is the piece that concerns the formula. The point is, this piece basically draws the wave function (or Ψ, if I'm correct) based on the given E, the energy value. The graph is also dependant on a given situation, described in fu g, which is in this example set to C (= 1) * x^2. The graph results in the wave function (being the square of the probability function?), which should, as fu gets incredibly high, get smaller and smaller and have an asymptote to zero. However, our E is not entirely precise, resulting in the graph only staying horizontal for a while, before having an asymptote either up or down. The more accurate the given E is (more decimals), the longer it 'stays' horizontal (before going up or down, but this can be ignored).

The perfect E is found on this idea;
The end of the graph holds information on where the nearest perfect E value is. If it's pointing up, the perfect E is actually above the given one right now. Based on this principle, our code determines the best decimal and then moves on to the next decimals, determining the E up until 15 decimals, and indeed, the more accurate decimals, the longer our graph stays flat.

This is my second question: Can I use this idea that the tail of the graph always points to the nearest energy value? It definitely seems to work, and I can calculate seemingly okay E-values, and the graph gets flatter on the end using this idea. I also got an alternative version where the graph is given points on how long it's flat, and based on that principle it finds the perfect E-value, but I'd prefer this one.

Another piece surrounding these basically leads to the next E-value, as in Schrödinger's Equation you got these multiple E-values corresponding with the fixed energy states the particle (electron) can be in.

All this together results in the following graph

Spoiler: Using the following start values

massfactor = 5;
C = 1;
wave = 1;
slope = 0;

The blue line represents the situation, x^2.
The orange one is the first energy state, the yellow one the second and the purple one the third (there are more), with the corresponding energy values of
E-1 = 0.4471758673234
E-2 = 4.0247949984073
E-3 = 7.6024556510729

I haven't got the slightest idea if this is correct, but hell, it seems okay?
(That's what my teacher said as well.)

This all leads to my third and main question: How do I expand this to a 2D (and eventually) a 3D environment to fully create a model of Hydrogen?
Is my current code usable (if correct at all), and expandable to a 2D environment? What needs to be done? Is it not usable at all, is it completely incorrect?
____________________________________________________________________

I would massively appreciate any help, even if not directly an answer to one of my questions.
I hope there's someone who knows about this and would be willing to help.

I did my best to explain my current idea, questions and model as well as possible, but I can easily imagine it's insufficient. I have had no experience beforehand with either this formula, or MATLAB.
Any further explanation I'd be happy to provide!

Also, if necessary, my code can be found https://drive.google.com/open?id=0B5_0h9k46DFuQXlGUUM3WmxhNEU.
Basically, if the Supporting Functions folder is included, quantumrun.m can be run without problems (supposedly). The code I described, which calculates the graph, can be found in calculate_fy.m

Thanks to anyone in advance!

Kind regards,

Isaiah van Hunen

 

[vanhees71]

First of all, you did a great job on yourself already. This is all very non-trivial stuff. What you discovered here is an algorithm to find the eigenfunctions for the Hamilton operator, known in the "computational physics" literature as the Fox-Goodwin or Numerov algorithm.

The math behind it is pretty tough to explain at the high-school level, but there is a huge mathematical theory behind it, which was developed in the early 20th century by Hilbert and many other mathematicians (like Fredholm) who thought about such eigenvalue problems of linear differential operators or integral operators on funcions a lot. Schrödinger just had to pick up these methods from the famous textbook on mathematical methods in physics by Courant and Hilbert when he developed his version of quantum theory "as a eigenvalue problem" ("wave mechanics"). The yellow one is that of $E_2$ and the purple one is that of $E_4$. So there should be energy eigenvalues $E_1$ and $E_3$ between these.

You can tune your algorithm to find these eigenvalues and eigenfunctions by counting the "nodes" (i.e., the zeroes) of the wave function. If you look for $E_1$ you choose a trial value and check whether your wave function has exactly one node. If it has no node, the value of $E$ must be larger if it has more than one it must be smaller. In this way you can find these states and eigenvalues. Of course, the exact value of $E_n$ must be such that in addition to the correct number of $n$ nodes the wave function must also go to 0 at infinity as you've figured already out.

For the harmonic oscillator the Hamiltonian has only discrete eigenvalues given by $E_n=\hbar \omega (n+1/2)$ with $n \in \mathbb{N}_0=\{0,1,\ldots \}$. As you have realized with your numerics the eigenvalues are those values for $E$ where the wave function goes to zero at infinity, and it must do so, because the wave function must be square integrable, i.e., the integral
$$\int_{-\infty}^{\infty} \mathrm{d} x |\psi(x)|^2$$
must exist, having a finite value. If you normalize the wave function such that this integral is 1, $|\psi(x)|^2$ is the probability density for finding the particle at position $x$.

Further the mathematics tells you that the eigenfunction for the eigenvalue $E_n$ has $n$ zeroes. So from your plot it's clear that indeed the orange curve is (a good approximation of) the ground state for $E_0=\hbar \omega/2$.

 

[IsaiahvH]

First of all, you did a great job on yourself already. This is all very non-trivial stuff. What you discovered here is an algorithm to find the eigenfunctions for the Hamilton operator, known in the "computational physics" literature as the Fox-Goodwin or Numerov algorithm.

The math behind it is pretty tough to explain at the high-school level, but there is a huge mathematical theory behind it, which was developed in the early 20th century by Hilbert and many other mathematicians (like Fredholm) who thought about such eigenvalue problems of linear differential operators or integral operators on funcions a lot. Schrödinger just had to pick up these methods from the famous textbook on mathematical methods in physics by Courant and Hilbert when he developed his version of quantum theory "as a eigenvalue problem" ("wave mechanics"). The yellow one is that of $E_2$ and the purple one is that of $E_4$. So there should be energy eigenvalues $E_1$ and $E_3$ between these.

You can tune your algorithm to find these eigenvalues and eigenfunctions by counting the "nodes" (i.e., the zeroes) of the wave function. If you look for $E_1$ you choose a trial value and check whether your wave function has exactly one node. If it has no node, the value of $E$ must be larger if it has more than one it must be smaller. In this way you can find these states and eigenvalues. Of course, the exact value of $E_n$ must be such that in addition to the correct number of $n$ nodes the wave function must also go to 0 at infinity as you've figured already out.

For the harmonic oscillator the Hamiltonian has only discrete eigenvalues given by $E_n=\hbar \omega (n+1/2)$ with $n \in \mathbb{N}_0=\{0,1,\ldots \}$. As you have realized with your numerics the eigenvalues are those values for $E$ where the wave function goes to zero at infinity, and it must do so, because the wave function must be square integrable, i.e., the integral
$$\int_{-\infty}^{\infty} \mathrm{d} x |\psi(x)|^2$$
must exist, having a finite value. If you normalize the wave function such that this integral is 1, $|\psi(x)|^2$ is the probability density for finding the particle at position $x$.

Further the mathematics tells you that the eigenfunction for the eigenvalue $E_n$ has $n$ zeroes. So from your plot it's clear that indeed the orange curve is (a good approximation of) the ground state for $E_0=\hbar \omega/2$.

First of all, thank you immensely for your reply!

So basically the idea of using the tail direction to determine the nearest $E_n$ value is correct?
Also, I made the mistake of labelling $E_0$ as $E_1$, and $E_1$ as $E_2$ and so forth.
(By the way, I also forgot to square $|\psi(x)|$, I changed that now)

Anyway, I changed my code to detect when a node is 'added' (my previous method is now unnecessary, I believe).

The current graph looks like this;

With
$E_0 = 0.4471758673234$
$E_1 = 2.2359736473542$
$E_2 = 4.0247949984073$
$E_3 = 5.8136234137805$

It seems okay, right?

If anyone has any idea on how to expand my code to turn this into a 2D model, that'd also be greatly appreciated.

For now, many thanks @vanhees17!

 

[vanhees71]

Now, I wonder why all curves seem to have $|\psi(0)|^2=1$ for all your wave functions. Since the potential is symmetric under reflections, i.e, $V(x)=V(-x)$ the wave functions should be even/odd for the even and odd $n$.

 

[A. Neumaier]

For hydrogen, you need to decompose your solutions into representations of the rotation group (partial wave analysis). This leaves for each value of angular momentum a scalar radial differential equation, to which your technique can probably be applied.

 

[IsaiahvH]

Now, I wonder why all curves seem to have $|\psi(0)|^2=1$ for all your wave functions. Since the potential is symmetric under reflections, i.e, $V(x)=V(-x)$ the wave functions should be even/odd for the even and odd $n$.

I understand why you ask that.
It correlates to the start value of $wave$ (see my earliest code), but not directly.
The previous screenshot is with a start value for $wave$ set to $1$.
For example, if the start value of $wave$ is set to $0.5$, it results in this graph, starting at exactly $0.25$. (https://s11.postimg.org/spvrrsamr/Wave_function.png)

Seemingly it only scales it up and down?
I am not entirely sure what to set for the start value of $wave$.

I would like to ask you, (1) what is the usual measurement of the probability density?
0% to 100%?

Also, (2) how come that for the higher energy values, the maximum probability density is higher than the maximum probability density for $E_0$?

Also, I wonder what you mean with

Since the potential is symmetric under reflections, i.e, $V(x)=V(-x)$ the wave functions should be even/odd for the even and odd $n$.

I agree that the graph is horizontally symmetric; and thus our teacher instructed us only to create the right side of the graph (and then you could mirror it). But what do you mean with 'the wave function being even/odd'? How is a function even or odd?

Kind regards and thank you immensely for your help so far,

Isaiah van Hunen

 

[vanhees71]

The Hamiltonian is invariant under spatial reflections, i.e., under the transformation $x \rightarrow -x$. That's why the eigenfunctions of energy are either even or odd functions, i.e., $\psi(x)=\psi(-x)$ (even) or $\psi(x)=-\psi(-x)$ (odd). It turns out that $\psi_n(x)=(-1)^n \psi_n(-x)$ for the harmonic oscillator. For the odd function you have to choose $\psi(0)=0$ and then you can choose any value for $\psi(h)$ (where $h$ is the stepsize of your numerics). For the even function you must have $\psi(0)\neq 0$ arbitrary. Since further $\psi(-h)=\psi(h)$ you can resolve the recursion of the Fox-Goodwin-Numerov algorithm to choose $\psi(h)$. See, e.g.,

http://www.iiserpune.ac.in/~pghosh/phy453/numerov.pdf

for details.

To give probabilities the wave function should be normalized to 1, i.e.,
$$\int_{-\infty}^{\infty} |\psi(x)|^2=1.$$
Of course, the probability to be in some finite interval is between 0 and 1, but the density can take any value!