12
$\begingroup$

Hey there and thanks for giving time to look at my question. This is a updated version of my question which I posted earlier in physics.stackexchange.com

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 where time in quantum mechanics is replaced by imaginary one $$ t = -i \tau $$ This substitution causes the high energy particles in my 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.

  • Note! The nonlinear Schrödinger equation includes the laplacian and some other differential terms in space. These are all solved using fast Fourier transform. In the end we have only a time ODE. *

My problem and question: The calculations go from $t_0$ to $t_f$. The ode45 is put in a for loop so it doesn't calculate a giant vector $[t_0,\dots,t_f]$ at the same time. The first round would start with ode45(odefun,$[t_0, t_0+\Delta/2, t_0 + \Delta],y,\dots$ ) and then go on next time from $t_0 + \Delta$. Here the time step $\Delta$ is my problem. Different choices in time steps gives me different ground state solutions and I have no idea how to determine which time step gives me the "most" correct ground state!

My attempt: I realize that in this scheme large time steps will cause large number of particles to decay before being re-normalized to the original number of particles whereas small time steps will cause smaller amount of particles to decay before getting re-normalized. My initial thought is that small time steps should give a more accurate solution but it seems to be the opposite.

I'm not a numerical expert so the choice of ode45 was simply arbitrary. ode113 gives me the same thing. :(

Does anyone have any thoughts on this matter. Let me know if any additional detail is needed.

Thank you.

Update 1: I've been researching the imaginary time method and the ODEs. It would seem that if the time step isn't small enough the whole thing becomes unstable. This makes me wonder if my nonlinear equations are stiff which makes things a lot more difficult from what I understand. I'll keep you updated.

Update 2: FIXED: The problem was indeed having the normalization outside of the ODE. If the normalization is kept inside odefun then the ODE returns the same outcome for different choices of "outside" time steps. My colleague showed me older codes and I simply added one line in my odefun.

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

    ...
    [ Nonlinear equations evaluated ]  
    ...


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

The last line calculates the difference in current number of particles (Ntemp) and the number of particles which the system should hold (N0). It adds a portion of the particles back to the output and thus creates a total particle number stability in the system instead of having them all decay away.

I will pose also a new question regarding the dimensionality of the problem and some differences in working with either picoseconds or nanoseconds as time steps in the ODE.

Thank you all. :)

Improve this question
$\endgroup$
5
  • 3
    $\begingroup$ The fundamental problem is that you're forcibly using an adaptive method like ode45() to take equispaced steps. Why, precisely, are you avoiding the generation of the "giant vector"? If you absolutely need equispaced points, have ode45() proceed as usual, and then use interpolation. $\endgroup$
    – J. M.
    Commented May 24, 2013 at 13:05
  • $\begingroup$ Hmm...this could be the problem. The origin of these fixed steps was that somewhere I needed to re-normalize the number of particles before they would all decay away. But maybe I can do this by putting the normalization in odefun and use the "giant time vector". Also, the input $y$ into the ode45 is 4*129*129 numbers. I was afraid if I didn't use time steps I wouldn't have enough memory. $\endgroup$
    – user4388
    Commented May 24, 2013 at 13:29
  • $\begingroup$ If memory serves, there should be an option in ode45() that would allow you to retain steps bigger than a certain threshold; you might want to look into that. $\endgroup$
    – J. M.
    Commented May 24, 2013 at 13:34
  • 1
    $\begingroup$ The answer is just to use a local error estimate. There is one built in to ODE45, so the easiest thing is to use that, but you could alternatively code up your own. $\endgroup$
    – David Ketcheson
    Commented May 24, 2013 at 15:57
  • 1
    $\begingroup$ In the answer to the follow-up question it turns out that $0.1$ is a dimensional quantity with dimensions $1/\text{time}$. Maybe more consistent results are obtained with $\frac{\alpha}{\Delta t} \cdot (N_t - N_0)$ where $\Delta t$ is the time step? $\endgroup$
    – Stefano M
    Commented Jun 20, 2013 at 22:56

3 Answers 3

Reset to default
4
$\begingroup$

Since you didn't post your MATLAB code, I'm not sure how you're calling ode45. I'm guessing you are changing the tspan vector (second argument) on each call to ode45. The first thing to understand is that the tspan vector has (almost) no effect on the time step used by ode45. The tspan vector simply allows you to pass to ode45 the time span of the integration and at what times you want output. The time step used by the Runga-Kutta algorithm in ode45 is adjusted internally to achieve a prescribed accuracy. The two parameters that control this accuracy are RelTol and AbsTol in the options structure passed to ode45. They have reasonable defaults and since you didn't mention these, I'm assuming you didn't change them.

I said "almost" no effect on the normal ode45 time step. If you are requesting output at time intervals very small relative to the time step ode45 would otherwise take, then it will have to reduce the time step to satisfy your output request. I believe this is what J.M. is assuming is going on. Why do you need the solution at so many output times? Usually it is sufficient to just request output at enough times to generate a smooth plot.

As to the change in solution you are seeing, maybe the default values of RelTol and AbsTol are not appropriate for your problem. I suggest replacing your loop on ode45 with a single call, request output at some reasonable number of times, and experiment with smaller values of RelTol and AbsTol until you get a converged solution.

Improve this answer
$\endgroup$
1
  • $\begingroup$ Thanks for the answer. The reason why I needed an solution at so many output times is because if the wave function isn't normalized regularly then everything decays and my system is empty. That's why I put the ode45 into a loop with small tspan vectors so I could re-normalize after each tspan vector. $\endgroup$
    – user4388
    Commented May 26, 2013 at 13:03
2
$\begingroup$

Since the nonlinear Schrödinger equation is, well, nonlinear, it can have a large number of stationary states, some of which may be stable. In the physical reality, starting from one particular state, the system will deterministically evolve into one final state. If the numerical scheme gives you different results for different discretizations (time steps), then this is fundamental flaw of your discretization. Check your code.

If you end up with a state $\psi_0$, then it is easy to verify if it is really a stationary state: If the time evolution is given by $$ \frac{\text{d}\psi}{\text{d}t} = F(\psi), $$ then $$ F(\psi_0) = 0. $$ Should you indeed end up with different stationary states, you can compare their Gibbs energies $$ G(\psi) = \int_\Omega E(\psi) $$ where $E(\cdot)$ is the energy density. When $F(\psi)=0$, $E(\psi)$ often looks quite simple, e.g., for Ginzburg-Landau equations, $E(\psi) = -|\psi|^4$.

Improve this answer
$\endgroup$
1
  • $\begingroup$ Yes. I plot the density profile of my output solution and when it doesn't change for a long time, basically stopped evolving, I assume I have reached a stationary state. But I'm not sure that looking at the energy density can help since the wave function is a spinor with (+2, +1, -1, -2) spin components. I don't think integrating each component will tell me the energy of the condensate but when I get to the ground state, the energy density should be stationary and thus a constant in time, this is a clue for a correct solution. $\endgroup$
    – user4388
    Commented May 27, 2013 at 9:08
1
$\begingroup$

Problem solved:

Normalization needs to be a part of the function evaluated in the ODE. Breaking up the ODE in many steps and normalizing between them causes seemingly numerical instability and produces different results depending on the time intervals the ODE is broken into. (See edit 2 in question for more details.)

Improve this answer
$\endgroup$

Your Answer

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