3
$\begingroup$

Hello there and thanks for taking a look at this problem.

This problem is related to my previous question and I will therefore use a similar introduction from, Choice of step size using ODEs in matlab.

Introduction: I'm currently studying a 2D exciton spinor Bose-Einstein Condensate and am curious about the ground state of this system. The mathematical method of getting to the ground state is called the imaginary time method.

The method is very simple. Normal real time in quantum mechanics is replaced by imaginary one $$ t = -i \tau$$ This substitution causes the high energy particles in the system to decay faster than the low energy ones. Re-normalizing the number of particles in every step of the calculation we end up with a system of lowest energy particles, aka. the ground state.

The equation(s) in question is nonlinear, called the nonlinear Schrödinger equation, sometimes the Gross-Pitaevskii equation. To solve the problem I'm using Matlabs ode45 which evolves the system forward in time and eventually reaches the ground state.

The problem and the question: There are a lot of parameters controlling the physics here. It is very appropriate to work with constants and variables of comfortable dimensionality such as, micrometers, picoseconds, millielectronvolts and so on.

However, working with picoseconds and millielectronvolts seems to produce drastically different results then working with nanoseconds and microelectronvolts. $$ (\text{ns},\mu \text{eV}) \quad ^?\leftrightarrow \ ^? \quad (\text{ps}, \text{meV})$$

My question is: If all parameters are set so that these two representations correspond correctly, i.e. are effectively the same thing, why would the results be completely different for say...

...in the representation of ps and meV

[~,y_out] = ode113(odefun,[0:100:10000],y_in, ...variables_in_ps_and_meV)...)

and otherwise in the representation of ns and ueV?

[~,y_out] = ode113(odefun,[0:0.1:10],y_in, ...variables_in_ns_and_ueV...)

As you can see, the time step 100ps corresponds to 0.1ns.

Note: In the case of using nanoseconds, i.e. time steps 0.1, the ODE calculates much, much faster. This has me inclined to believe the method of using nanoseconds is somehow rushed and not accurate.

Thanks for any thoughts!

Update 1: Here are two pictures showing the completely different results

Update 2: I just found out that it doesn't matter at all whether I work with $\mu$eV or meV. The same results are displayed! This narrows the problem down to only the choice of appropriate time dimensionality for the ode.

Update 3: Major breakthrough. I decided to write a simple 1D Schrödinger equation in a harmonic potential and tried out nanoseconds vs. picoseconds. The results can be found here.

As you can see, there is no difference when using either nanoseconds or picoseconds (not very surprising). But now just remains to figure out what in the devil is causing this strange behavior in my 2D nonlinear Schrödinger equations.

Improve this question
$\endgroup$
4
  • $\begingroup$ Do you mean that you are not reaching the same ground state at the end ? Or that the path to the ground state is different ? $\endgroup$
    – nat chouf
    Jun 19, 2013 at 21:34
  • $\begingroup$ Well, at least the path is drastically different and I cannot see it will end up in the same ground state but I could be wrong. Here are pictures ( imgur.com/a/1XUCg ) showing the results at 1ns but one using picoseconds and another nanoseconds. I'll update it when I have calculated forward to 100ns. $\endgroup$
    – user4388
    Jun 20, 2013 at 16:07
  • 1
    $\begingroup$ If you are simply changing the units on the equation you are working with, an explicit ODE solver like ode45 should produce identical results. With an implicit solver you may run into some conditioning issues, but I imagine those would be minor. Also the time step you define to matlabs ode functions is not the actual time step it uses, it simply interpolates onto your grid for your output. $\endgroup$
    – Godric Seer
    Jun 20, 2013 at 18:42
  • $\begingroup$ This is true that the ODE uses it's own time steps to match the accuracy. However the way I go about calculating it from say...0 to 100ns is I create a loop where I let the ODE calculate for (0 -> 0.1ns) and then (0.1 -> 0.2ns) and so on. This method is also very accurate and allows me plot results and observe the progress of calculations. $\endgroup$
    – user4388
    Jun 20, 2013 at 19:38

2 Answers 2

Reset to default
4
$\begingroup$

You should rewrite the equation dimensionless with dimensional analysis. Proper step size should be evident then and can be compared more easily.

Improve this answer
$\endgroup$
3
  • $\begingroup$ By proper you mean steps around 1.0? Yes of course if I were to use meters, seconds and Joules then everything would go haywire since the quantities would be extremely small numbers. $\endgroup$
    – user4388
    Jun 20, 2013 at 19:40
  • 2
    $\begingroup$ +1 Correct suggestion: it was indeed a dimension error in the normalization, see HappyDog's own solution. $\endgroup$
    – Stefano M
    Jun 20, 2013 at 22:32
  • $\begingroup$ @HappyDog The point of rewriting the equation dimensionless is to be able to arbitrarily scale the system to a numerical stable region without bothering about Joule, seconds or nanoseconds (e.g. like saying "hbar = 1") because their relative relations are preserved. If relative system scales are in the correct order of magnitude there won't be any extrem small numbers. $\endgroup$
    – Bort
    Jun 21, 2013 at 6:43
2
$\begingroup$

BOOM SOLVED!!

Sorry, I'm just super happy about finding the one line who that the honor of making me question reality for the past five days.

Inside the odefun there is a line which was suggested by a friend of mine, who has great experience in these things (he's in India doing shenanigans, otherwise I would have contacted him right away).

This line restores a portion of the particles back to the system which is extremely important in the imaginary time method. Otherwise the particles decay all away leaving an empty system. This line was the solution to my previous problem regarding choice of time step size and normalization.

function y_out = odefun(t,y_in,...variables...) 

    ...
    [ Nonlinear equations evaluated ]  
    ...


    y_out = y_out + 0.1*y_in*(N0-Ntemp) ;
end

I knew that the choice of 0.1 in the last line was never perfectly arbitrary but I had no idea it was so explicitly connected to the type of time steps used in the ODE. The results for nanoseconds showed that a great portion of the particles seem to have decayed away, so I changed 0.1 to 100 since that is the order between picoseconds and nanoseconds and voila! Exactly same results as for the case where I used picoseconds.

Thanks for all who gave this problem a thought. I'm just happy it's solved!

Improve this answer
$\endgroup$
2
  • 3
    $\begingroup$ ehm... sorry for the late comment but dimensional analysis should have told you that 0.1 is not a pure number but has dimensions of $[T]^{-1}$ where $[T]$ is time. In fact y_out stands for $y' =\mathrm{d}y/\mathrm{d}t$ while y_in corresponds to $y$. Since N0-Ntemp is a pure number you cannot sum y_out and y_in without multiplying by a $[T]^{-1}$ constant. So very simply $10^{-1}\, \mathrm{ps}^{-1} = 10^2 \, \mathrm{ns}^{-1} = 10^{11}\,\mathrm{s}^{-1}$. $\endgroup$
    – Stefano M
    Jun 20, 2013 at 22:29
  • $\begingroup$ Thanks! I hadn't thought of it that way. This makes things much clearer. $\endgroup$
    – user4388
    Jun 21, 2013 at 0:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.