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

[rchretien]

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.