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.
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) ;
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. :)