weird initial state probability at t=0, but otherwise fine for the rest of the timesteps?

[dondondooooon]

Im getting a huge spike at the start of my probability, and I can't figure out why... This is by the way for the first time I switched on the double slit wall in the box. (Problem 7 part 2)

However, I did get a seemingly correct graph for when the slit was switched off... (Problem 7 part 1)

[dondondooooon]

I do have several speculations, but I don't know which one it is, and I'm afraid to change too much stuff permanently before consulting for help... My biggest speculation as to why im getting this issue is that:

In my c++ program, I use a cube matrix to store the different u_ij for a chosen n. When I initialize this cube matrix (i.e. for the first slice) I add in the potential (including the one from the middle wall in the box where the slit is, if the slit is switched on), before I actually initialize and calculate the inner point matrix U using the gaussian wavepacket function. So I was thinking doesn't this mean that the potential is taken into account here as well? since the high potential slit wall is within the inner points? Thus resulting in perhaps a huge probability at the initial state of the system.

(p.s. unsure if I normalize the function properly in the wavepoint function as well, but im pretty sure I did?)...

Anyways, am I suppose to not add the potential to the cube matrix at all? And just have U in the inner points w/o the potential, whilst my A and B matrix will take the potential into account when updating the points in U?

[dondondooooon]

Nvm, i think i found out the problem. for a given time slice, I was plotting the probability of the whole matrix instead of only the inner points of the matrix :))) thus giving a weirdly high spike in probability.

Though would still like to know if we were meant to add the potential to the cube matrix that we are working with...

[anderkve]

Hi @dondondooooon!

Here are some quick comments:

• If you want to store the (constant) potential V as well, e.g. for plotting, I'd recommend storing that in a separate matrix from U, and write it to a separate file. Keep in mind that the potential V and the wavefunction U are fundamentally different types of quantities, and while you can of course store them in memory however you like, these quantities should never be added mathematically. (If we hadn't scaled away all units these quantities would have different units, so adding them would be meaningless.) So I'm a bit worried about lines like these

  grid_tid.slice(0) = V;
  U = grid_tid.slice(0)(span(1,L),span(1,L));

  since the second line here effectively sets U equal to (the inner part of) V, which can't make sense.

• This probably explains the plotting issue you had, as your code seems to have been saving V along the "border" of your matrix (at every time slice), while replacing the inner points with U. I don't think that's a sensible way of saving the data. I'd suggest you use your cube to only save the wavefunction U at every time step. Whether you just want to save the inner points of U or the complete U (including the "border" of boundary points) is up to you -- though I'd recommend saving the complete U since you know exactly the value of U along the border at every time step.

• In your function quantum::func you construct a Gaussian wavepacket based on using integers i and j for the x and y coordinates, in the factors (i - xc) and (j - yc). But that can't be right. What you should have are factors like (x-xc) and (y-yc) where both x and xc are doubles in the range [0,1] (and similarly for y and yc). Keep in mind not to confuse the value of the coordinate x with the value of the corresponding index i. (See the Notation section of the project description.)

• If the y axis of the last plot is showing 1 - [total probability at time t], the variation looks much too large to be correct. (But perhaps this was just part of the plotting issue you solved.) Also, keep in mind that the y-axis label can't be p_{ij}^n if what you are plotting is the total probability across the entire xy plane at timestep n.

[dondondooooon]

Yeah, I suspected that V and U are suppose to be seperate from one another, but I made a mistake while coding thinking of how the plotting would look like at the end, and thought of adding in V in the main matrix. Thank you for the help and clarification. :))

I dont quite get what you meant by saving the inner points of U vs. the complete U (including the "border" of boundary points). Is the border in this case just 0s according to the boundary conditions?

Did you also mean to have xc and yc as a double in the range [0,1] (i.e. in a vector) or is it only x and y that needs to be in a vector [0,1]? I thought xc and yc were suppose to be a constant double which represents the coordinates of the centre of the initial wave packet?

[anderkve]

I dont quite get what you meant by saving the inner points of U vs. the complete U (including the "border" of boundary points). Is the border in this case just 0s according to the boundary conditions?

Yep! :) So it doesn't matter too much, but personally I find it easier to not get confused with indices etc (especially in making the final plots) if my programs always store the complete state, regardless of what the boundary conditions happen to be in the given problem.

Did you also mean to have xc and yc as a double in the range [0,1] (i.e. in a vector) or is it only x and y that needs to be in a vector [0,1]? I thought xc and yc were suppose to be a constant double which represents the coordinates of the centre of the initial wave packet?

My point was not that anything had to be in a vector, but that the numbers x, xc, y and yc are all floating-point numbers that live on the domain from 0 to 1 (since this is the coordinate system we decided to work with). You are right that xc and yc are simply constants. The main point was that the bracket (i - xc) necessarily would give the wrong results since i would be some large integer while xc would be e.g. xc = 0.25, while mathematically this bracket should have the form (x - xc) where x is some number between 0 and 1 (and similarly for the y coordinate).

[dondondooooon]

My point was not that anything had to be in a vector, but that the numbers x, xc, y and yc are all floating-point numbers that live on the domain from 0 to 1 (since this is the coordinate system we decided to work with). You are right that xc and yc are simply constants. The main point was that the bracket (i - xc) necessarily would give the wrong results since i would be some large integer while xc would be e.g. xc = 0.25, while mathematically this bracket should have the form (x - xc) where x is some number between 0 and 1 (and similarly for the y coordinate).

Am I constructing the wavepacket correctly now?

where xvec and yvec is defined in the constructor:

[anderkve]

I think so -- at least I can't spot anything wrong right away :)

[dondondooooon]

Since we're working with dimensionless variables, should I just have "t" on the x axis without specifying what t is (i.e. it is time). Also is the labelling on the y-axis alright considering I had to scale them. Is there a cleaner way, if this isn't acceptable?

sim1_prob.pdf

sim1_dev.pdf

[anderkve]

(I'll write p^n_total instead of the sum over i,j, but the sum notation is perfectly fine.)

When using bracket notation, a bracket "[X]" is typically read as "in units X". So in your sim1_dev.pdf plot I would have read your current y axis as "1 - p^n_total, in units 10^15", which is not what you wanted to convey. I suspect what you wanted so say was either

• "(1 - p^n_total) x 10^15"; or
• "1 - p^n_total [10^{-15}]".

Either of those would be correct. (Though, if you just set something like plt.xlim([-8e-15, 8e-15]) I think you don't have to worry about any manual scaling, since Matplotlib will show the factor 10^{-15} automatically either directly in the tick numbers or on top of the axis.)

Similarly, the axis label "t [time]" would be read as "t in units time", which doesn't make sense. Here I'd simply go with "t" or "time", and if you want to stress that it is dimensionless you could put something like "t [1]" or "t [dimensionless]" or something similar. But I'd say it's perfectly fine to just say "t" and point out elsewhere in the report that "t" is a dimensionless time variable.

[dondondooooon]

Last thing... is setting up the a & b vector values for the diagonal values of A and B matrices used in the Crank Nicolson method correct here?

if i vectorize (translate my V matrix into v vector) using V.as_col like this

instead of using my function that specifically does this in accordance to Problem 2

[anderkve]

The easiest way to double-check code like this is simply to print the matrix to screen (you can print Armadillo matrices with cout) and see if the values are what you expect. (To make the output fit in the terminal, either print a subset of the matrix or set a very large step size, which gives small matrices, and make the terminal font size sufficiently small.) Or in this case, you could simply print the diagonal as a vector. If you at the same time print the values of the variables k, r, dt_half and v(k) it's easy to check if the resulting diagonal looks like it should.

Looking at the code above, I'm worried about this line in quantum:af (and the corresponding one in quantum::bf`):

cx_double send( 1.0, 4*imag(r) + imag(dt_half) * v(k) ); 

In particular, I wouldn't expect the variable dt_half to have an imaginary component -- the time step should simply be a real number. So in that case imag(dt_half) would return 0 and thus make the diagonal wrong.

Personally I find it easier to use the complex number "i" directly in my code, rather than using the cx_double (in reality the std::complex<double>) constructor with the real and imaginary parts as separate arguments. At least I find that I make fewer mistakes when I write code that looks more similar to the math. So I would do something like this

cx_double send = 1.0 + 4 * r + 1.0i * dt_half * v(k); 

where r is of type cx_double while dt_half is of type double.

If you encounter compilation problems when using the "i" notation, take a look at this page: https://stackoverflow.com/questions/17925904/how-to-use-complex-number-i-in-c

Concerning the ordering in your vectorized potential: Look at the documentation of the as_col() function: http://arma.sourceforge.net/docs.html#as_col_row There it says that the vector is created "by concatenating all columns". Whether or not that corresponds to your function for translating from (i.j) --> (k) is something you need to figure out, so I won't give a straight yes/no answer. :)

(And if you're unsure, remember that you can always just test it by setting up a small example matrix, e.g. a random 3x3 matrix, and checking the resulting .as_col() vector against the element ordering you get with your index translation function.)

[dondondooooon]

In particular, I wouldn't expect the variable dt_half to have an imaginary component -- the time step should simply be a real number. So in that case imag(dt_half) would return 0 and thus make the diagonal wrong.

I initialized dt_half as a complex double with only imaginary...

So when I do imag(dt_half) it will return dt/2. In hindsight it's a bit roundabout way of doing it, but I started with this and just continued accordingly.