Verfasst am: 22.02.2014, 18:28
Titel: Calculation issue with octave (data unexpectedly go to zero
Hi everybody,
I contact you because of my poor knowledge of octave. For my master thesis, I have to use Octave installed on a remote cluster whilst I am running matlab at home. For the same code, matlab and octave do not produce the same results. Or rather, should I precise, they produce the same results during a certain amount of time, until octave produces an unexpected collapse of the data to zero.
Here is an example of what happens :
Under matlab (I repeat that I have used the same version of the code), the output is basically the same, except for abscissae greater than ~975 (i.e. the abscissa from where octave collapses) where no crash is observed (contrarily to what can be seen here above). Roughly speaking, the code behaves as expected when run under matlab and not when run under octave.
My question is, since this is manifestedly not a coding problem (in which case, the output from matlab would also collapse), is this problem linked to an option, or some kind of precision, or to a version problem, in octave ? Has this problem already been observed ? If so, how can I fix it ?
Im looking forward to reading your opinions and suggestions,
Renaud
Verschoben: 22.02.2014, 20:13 Uhr von denny Von Programmierung nach Octave-Forum
Verfasst am: 22.02.2014, 22:54
Titel: Re: Calculation issue with octave (data unexpectedly go to z
Dear rchretien,
Without seeing any detail of the code, it is impossible to guess the reason for the observed behaviour. Is the applied algorithm is numerically instable, the breakdown could be correct in a sense of computational theory.
There is no problem for me posting the code, apart from the fact that it is a relatively massive one (approximately 20 functions, some of which being not straightforward) so that I am afraid that no one reads it (which is, by the way, totally logical).
However, I do not think that the algorithm is numerically unstable. Let me explain why. When I set an arbitrary duration of simulation, the collapse seems to occur at a given time, let us denote it t_0. Now, if I set a duration of simulation, say 100 times, greater than the previous one, then the collapse will NOT occur at time t = t_0 but rather at time t = 100*t_0 (although I am not sure of the proportionality factor, this is a visual observation). Therefore, if the algorithm was unstable (of course, I have checked this out a number of times and found nothing) then I can safely assume that the collapse would be insensitive to the duration of simulation and could occur at the same time, t_0. For example, I have launched the code for a tiny duration of simulation (1/100 of the one in the figure in my previous message), and the same behaviour was observed : collapse at ~1/2 of the duration of simulation.
Another argument in favour of this point of view is that the same launch of the code with matlab (that is, exactly the same version of the code with exactly the same simulation parameters) procudes an expected result, that is, the same result as octave wi NO collapse.
Then, because I know octave poorly, I wonder if this problem is not linked to an option to set in octave (f. ex. yesterday, I thought that the "save" command of octave produced same output as the one of matlab, but the .mat files procudes were unreadable under matlab and that because an option was to be set in octave).
The fact that the values are not rigorously equal before the collapse comes from the pseudo-random feature of the method : it is a Monte-Carlo simulation and the code must be launched a huge amount of times and the results are to averaged.
PS : compared to the first image I posted, only the duration of simulation has changed.
I doubt that we will be able to help you unless you can identify the source of the discrepancy.
If you have 20 functions involved, try testing them individually to see which one is causing the problem.
Zitat:
When I set an arbitrary duration of simulation, the collapse seems to occur at a given time, let us denote it t_0. Now, if I set a duration of simulation, say 100 times, greater than the previous one, then the collapse will NOT occur at time t = t_0 but rather at time t = 100*t_0
If you set the duration to 100 times the previous one, will your algorithm execute the same way up to the original duration? It might just be that you are using a certain number of simulation steps and that the simulation is just getting rescaled this way. This does not eliminate the possibility of a numerical problem.
Just guessing, but that's all we can do unless you give us something to work with.
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)
% 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
% Beginning ot the temporal loop
phi_f = phi_i;
for t = 0:dt:duration
ifmod(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);
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.
are you absolutely sure you are calling the function with the same input for duration and dt both times?
Using the Debugger may help you track down the issue further.
I personally do not like for-loops of that style, I'd prefer
Code:
allt = 0:dt:duration;
for n=1:numel(allt)
t = allt(n) ... end
may not a good idea (in particular if recordingTime is not an integer) as minor inaccuracies in t may result in things not being equal when they should be. I would instead use
Since the problem treated here is the spontaneous emission of a photon by an atom, the expected population of the ground state is expected to be one when the time is large (this explanation is given just in case someone tries to understand the code).
I do not have Octave installed, so I cannot test that.
Try debugging on Octave. I'd expect it to go through all iterations of the for-loop but the if-statement will likely evaluate to false.
Zitat:
I am using octave (whilst I am using matlab at home) because the clusters on which I have to work have no matlab license.
That could be changed, couldn't it? I'm sure the folks from MathWorks could help you there
As you suspected, that cannot be changed, I have to use it I did not asked the question to the MathWorks guys since the problem was concerning Octave (and I do not know whether they would have helped me).
As you correctly suspected, the problem comes from the modulo dealing with non-integer values. The function may exhibit an unexpected behaviour in such a case, as the following example demonstrates :
Code:
octave:20> mod(0.2, 0.2) ans = 0
octave:21> mod(0.4, 0.2) ans = 0
octave:22> mod(0.6, 0.2) ans = 0.20000 # ??????????????
The guys helping me in the topic which I gave the link herebefore have reported a possible bug for mod. One of them proposed me to use this function instead of mod.
Although this function seems to be doing the job (I have to perform further checkings just in case) I am still following the discussions. I shall post the outcome here (and change the name of the original topic) for completeness.
As you suspected, that cannot be changed, I have to use it Smile I did not asked the question to the MathWorks guys since the problem was concerning Octave (and I do not know whether they would have helped me).
Slight misunderstanding - my recommendation was to try and get MATLAB to run on these clusters, and MathWorks could have helped you with that.
My recommendation had a little flaw, it should have been more like
Code:
ifmod(t, recordingTime)) < tol || mod(t, recordingTime)) > recordingTime - tol
This should eliminate the need for a separate function.
The behavior of Octave is not as surprising as it may look like. Numbers like 0.2 do not have an exact representation in the binary system which introduces minimal numerical errors.
Even in MATLAB (and likely any other package that relies on numerical computations):
I was aware of the fact that decimal numbers do not have an exact representation in the binary system but I thought that was managed by the software. Indeed, it was the case under matlab where the mod did not yield any unexpected behaviour but it was not the case under octave. The good point is that I shall be aware of that from now on, the bad point is that it is surprising for someone using matlab.
Anyway, I think that I can mark the topic as resolved. Any suggestion for renaming it with a more appropriate name ?
Best regards,
Renaud
Einstellungen und Berechtigungen
Du kannst Beiträge in dieses Forum schreiben. Du kannst auf Beiträge in diesem Forum antworten. Du kannst deine Beiträge in diesem Forum nicht bearbeiten. Du kannst deine Beiträge in diesem Forum nicht löschen. Du kannst an Umfragen in diesem Forum nicht mitmachen. Du kannst Dateien in diesem Forum posten Du kannst Dateien in diesem Forum herunterladen
MATLAB, Simulink, Stateflow, Handle Graphics, Real-Time Workshop, SimBiology, SimHydraulics, SimEvents, and xPC TargetBox are registered trademarks and The MathWorks, the L-shaped membrane logo, and Embedded MATLAB are trademarks of The MathWorks, Inc.
I'm sure the folks from MathWorks could help you there 
RSS