Re: Calculation issue with octave (data unexpectedly go to zero at some point)

[Juan Pablo Carbajal]

On Thu, Feb 27, 2014 at 8:39 PM, rchretien <address@hidden> wrote:
> Hi Michael,
>
> I can easily imagine that my question sounds like "I have a problem, how can
> I fix it", it is anyway what I was thinking when writing it.
>
> However, unless someone particularly insists for me posting the code, I did
> not think of publishing it (because I seriously doubt someone wants to read
> it, which is completely normal).
>
> Anyway, although I did not have much time these last days, I think I am a
> bit closer of finding what goes wrong. Indeed, the zeros in my data come
> from the instanciation of the vector (which is initially full of zeros and
> supposed to be gradually filled with results during the simulation). The for
> loop which is supposed to compute the results stops before the end, which
> amounts to let the rest of the vector filled with zeros. The problem is that
> I have absolutely no idea of what mechanism could stop it, except one break
> (if it were to happen, a message would be displayed at the screen,
> indicating that a problem occured). Here is the associated code (which is
> short, reason why I do not hesitate to post it)
>
> [code]
> function [phi_f, energy, squareOfEnergy, expectedPopulationOfExcitedState,
> expectedPopulationOfGroundState, expectedp, expectedpSquare] =
> temporalEvolution(phi_i, hbar, duration, dt, recordingTime, HS, H, Cm, P_e,
> P_g, ge, gg, Je, Jg, dpmax, k_L)
>
> % The purpose of this function is to compute the evolution for the
> % wavefunction between t and t + dt by determining if a quantum jump occurs
> % or not and calling the proper functions with respect to the result of the
> % jump/no-jump event.
> % The output is the NORMALISED state vector, expressed in the {|Je, me, p>,
> |Jg, mg, p>} basis.
> % The duration is assumed to be given in terms of integer multiple of
> 1/Gamma
>
> % Definition of two vectors containing the values for the expected energy
> % and its square
> energy = zeros(1 + duration/recordingTime, 1);
> squareOfEnergy = energy;
> expectedPopulationOfExcitedState = zeros(1 + duration/recordingTime, 4);
> expectedPopulationOfGroundState = zeros(1 + duration/recordingTime, 2);
> expectedp = zeros(1 + duration/recordingTime, 1);
> expectedpSquare = zeros(1 + duration/recordingTime, 1);
> [p, pSquare] = buildMomentumAndItsSquare(hbar, ge, gg, Je, Jg, dpmax, k_L);
>
> n = 1; % Counter for the records of mean energy, recorded all integer
> multiples of recordingTime
>
> size = 2*dpmax + 1;
>
> P_e1 = zeros(length(P_e));
> P_e2 = zeros(length(P_e));
> P_e3 = zeros(length(P_e));
> P_e4 = zeros(length(P_e));
> P_g1 = zeros(length(P_e));
> P_g2 = zeros(length(P_e));
>
> P_e1(1:size,1:size) = P_e(1:size,1:size);
> P_e2(size+1:2*size,size+1:2*size) = P_e(size+1:2*size,size+1:2*size);
> P_e3(2*size+1:3*size,2*size+1:3*size) =
> P_e(2*size+1:3*size,2*size+1:3*size);
> P_e4(3*size+1:4*size,3*size+1:4*size) =
> P_e(3*size+1:4*size,3*size+1:4*size);
> P_g1(4*size+1:5*size,4*size+1:5*size) =
> P_g(4*size+1:5*size,4*size+1:5*size);
> P_g2(5*size+1:6*size,5*size+1:6*size) =
> P_g(5*size+1:6*size,5*size+1:6*size);
>
> % Beginning ot the temporal loop
> phi_f = phi_i;
> for t = 0:dt:duration
>     if mod(t, recordingTime) == 0
>        energy(n,1) = real(expectedValue(phi_i, HS));
>        squareOfEnergy(n,1) = real(expectedValue(phi_i, HS))^2;
>        expectedPopulationOfExcitedState(n,1) = real(expectedValue(phi_i,
> P_e1));
>        expectedPopulationOfExcitedState(n,2) = real(expectedValue(phi_i,
> P_e2));
>        expectedPopulationOfExcitedState(n,3) = real(expectedValue(phi_i,
> P_e3));
>        expectedPopulationOfExcitedState(n,4) = real(expectedValue(phi_i,
> P_e4));
>        expectedPopulationOfGroundState(n,1) = real(expectedValue(phi_i,
> P_g1));
>        expectedPopulationOfGroundState(n,2) = real(expectedValue(phi_i,
> P_g2));
>        expectedp(n,1) = real(expectedValue(phi_i, p));
>        expectedpSquare(n,1) = real(expectedValue(phi_i, pSquare));
>        n = n + 1 % This line echoes, which is how I knew that the loop does
> not go as far as expected
>     end
>     dp = real((dt*1i/hbar) * expectedValue(phi_i, H - H')); % Probability
> for a jump to occur
>     if dp > 0.1
>         disp('Caution, computation might be meaningless');
>         break;
>     end
>     if dp < rand(1) % In such a case, no jump occurs
>         phi_f = nonHermitianEvolution(phi_i, hbar, dt, H);
>     else
>         m = computeWhichJumpOccurs(phi_i, dt, dp, Cm);
>         phi_f = quantumJumpEvolution(phi_i, Cm, m);
>     end
>     phi_i = phi_f/norm(phi_f); % Normalisation
> end
> phi_f = phi_f/norm(phi_f);
>
> end
> [/code]
>
> For the parameters I have chosen, n is assumed to be equal to n = 51 at the
> end of the for loop (it is the case under matlab). Under Octave, n = 27 at
> the "end" and stops increasing after (which indicates that the loop is
> over...). As I said, since I see no mechanism responsible for unexpectedly
> stopping the loop (such as break f. ex., at the restriction I have mentioned
> herebefore), I do not understand what is going on.
>
> Best regards,
>
> Renaud
>
> PS : After having commented the "break" statement just in case of, it does
> not change anything.
>
>
>
> --
> View this message in context: 
> http://octave.1599824.n4.nabble.com/Calculation-issue-with-octave-data-unexpectedly-go-to-zero-at-some-point-tp4662177p4662436.html
> Sent from the Octave - General mailing list archive at Nabble.com.
> _______________________________________________
> Help-octave mailing list
> address@hidden
> https://mailman.cae.wisc.edu/listinfo/help-octave

Here in the list, if the code shows the problem we are ready to read
it. However, the code you provided does not work, so we cannot test to
see if we can reproduce your problem. Your code is filled with calls
to other functions you do not provide.
So when we ask for code we mean: "A minimalistic code that works and
shows the problem". You should forgive us if we weren't so clear about
that, you know, we do this all the time :D

We ask for a minimalistic code example, because in the work of writing
that code you most likely will find the problem yourself. In a way we
are helping you by asking for this!

So, can you isolate the problem you have an send a couple of lines
illustrating it?

Thanks