WICHTIG: Der Betrieb von goMatlab.de wird privat finanziert fortgesetzt. - Mehr Infos...

Mein MATLAB Forum - goMatlab.de

Mein MATLAB Forum

 
Gast > Registrieren       Autologin?   

Partner:




Forum
      Option
[Erweitert]
  • Diese Seite per Mail weiterempfehlen
     


Gehe zu:  
Neues Thema eröffnen Neue Antwort erstellen

Calculation issue with octave (data unexpectedly go to zero

 

rchretien
Forum-Newbie

Forum-Newbie


Beiträge: 7
Anmeldedatum: 22.02.14
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 22.02.2014, 18:28     Titel: Calculation issue with octave (data unexpectedly go to zero
  Antworten mit Zitat      
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 ?

I’m looking forward to reading your opinions and suggestions,

Renaud
Private Nachricht senden Benutzer-Profile anzeigen
Verschoben: 22.02.2014, 20:13 Uhr von denny
Von Programmierung nach Octave-Forum


Jan S
Moderator

Moderator


Beiträge: 11.057
Anmeldedatum: 08.07.10
Wohnort: Heidelberg
Version: 2009a, 2016b
     Beitrag Verfasst am: 22.02.2014, 22:54     Titel: Re: Calculation issue with octave (data unexpectedly go to z
  Antworten mit Zitat      
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.

Could you post any details about your code?

Kind regards, Jan
Private Nachricht senden Benutzer-Profile anzeigen
 
rchretien
Themenstarter

Forum-Newbie

Forum-Newbie


Beiträge: 7
Anmeldedatum: 22.02.14
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 23.02.2014, 10:49     Titel:
  Antworten mit Zitat      
Dear Jan,

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

Thank you for your help,

Renaud
Private Nachricht senden Benutzer-Profile anzeigen
 
rchretien
Themenstarter

Forum-Newbie

Forum-Newbie


Beiträge: 7
Anmeldedatum: 22.02.14
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 23.02.2014, 11:05     Titel:
  Antworten mit Zitat      
Here is an illustration of what I said.

The output under matlab :


expectedPopulationOfExcitedState =

0.0000 0 0 0
0.0064 0 0 0
0.0225 0 0 0
0.0451 0 0 0
0.0709 0 0 0
0.0974 0 0 0
0.1227 0 0 0
0.1451 0 0 0
0.1633 0 0 0
0.1767 0 0 0
0.1845 0 0 0
0.1867 0 0 0
0.1836 0 0 0
0.1756 0 0 0
0.1638 0 0 0
0.1493 0 0 0
0.1334 0 0 0
0.1173 0 0 0
0.1022 0 0 0
0.0891 0 0 0
0.0786 0 0 0
0.0710 0 0 0
0.0663 0 0 0
0.0644 0 0 0
0.0649 0 0 0
0.0674 0 0 0
0.0713 0 0 0
0.0760 0 0 0
0.0811 0 0 0
0.0862 0 0 0
0.0908 0 0 0
0.0947 0 0 0
0.0977 0 0 0
0.0998 0 0 0
0.1008 0 0 0
0.1008 0 0 0
0.1001 0 0 0
0.0986 0 0 0
0.0967 0 0 0
0.0945 0 0 0
0.0921 0 0 0
0.0899 0 0 0
0.0878 0 0 0
0.0860 0 0 0
0.0846 0 0 0
0.0836 0 0 0
0.0831 0 0 0
0.0828 0 0 0
0.0829 0 0 0
0.0832 0 0 0
0.0837 0 0 0

The output under octave :



expectedPopulationOfExcitedState =

0.0000 0 0 0
0.0451 0 0 0
0.0974 0 0 0
0.1227 0 0 0
0.1767 0 0 0
0.1845 0 0 0
0.1836 0 0 0
0.1756 0 0 0
0.1493 0 0 0
0.1022 0 0 0
0.0786 0 0 0
0.0644 0 0 0
0.0649 0 0 0
0.0674 0 0 0
0.0713 0 0 0
0.0908 0 0 0
0.0947 0 0 0
0.1001 0 0 0
0.0986 0 0 0
0.0945 0 0 0
0.0921 0 0 0
0.0831 0 0 0
0.0828 0 0 0
0.0829 0 0 0
0.0832 0 0 0
0.0837 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0

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.
Private Nachricht senden Benutzer-Profile anzeigen
 
Harald
Forum-Meister

Forum-Meister


Beiträge: 24.434
Anmeldedatum: 26.03.09
Wohnort: Nähe München
Version: ab 2017b
     Beitrag Verfasst am: 23.02.2014, 11:16     Titel:
  Antworten mit Zitat      
Hi,

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.

Best wishes,
Harald
Private Nachricht senden Benutzer-Profile anzeigen
 
rchretien
Themenstarter

Forum-Newbie

Forum-Newbie


Beiträge: 7
Anmeldedatum: 22.02.14
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 27.02.2014, 20:35     Titel:
  Antworten mit Zitat      
Hi Harald,

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
 


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.
Private Nachricht senden Benutzer-Profile anzeigen
 
Harald
Forum-Meister

Forum-Meister


Beiträge: 24.434
Anmeldedatum: 26.03.09
Wohnort: Nähe München
Version: ab 2017b
     Beitrag Verfasst am: 27.02.2014, 21:30     Titel:
  Antworten mit Zitat      
Hi,

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


The statement
Code:
if mod(t, recordingTime) == 0

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
Code:
if abs(mod(t, recordingTime)) < tol

with some suitable tolerance tol.

Ultimately. I'd say: if MATLAB gives you the correct results and Octave doesn't, why not just use MATLAB? :)

Cheers,
Harald
Private Nachricht senden Benutzer-Profile anzeigen
 
rchretien
Themenstarter

Forum-Newbie

Forum-Newbie


Beiträge: 7
Anmeldedatum: 22.02.14
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 27.02.2014, 23:24     Titel:
  Antworten mit Zitat      
Hi Harald,

I am using octave (whilst I am using matlab at home) because the clusters on which I have to work have no matlab license.

I am running out of time for tonight (but I am going to test what you told about tomorrow, no mistakes) however, this might help.

I have posted a much (much !) simpler version of my code here :

https://www.dropbox.com/sh/hiuv8u4jxqk2hs3/tI_GfBHJ3Y

Simply uncompress the archive, composed of a minimalistic version of the code. To launch it, use the function called main that way :

[expectedPopulationOfGroundState] = main(N)

where N is a parameter denoting the number of realisations of the method (I suggest to put N=1, since the problem does not depend on N).

Then, look at what happens when launched under Octave :

Zitat:

expectedPopulationOfGroundState =

0.00000
0.00000
1.00000
1.00000
1.00000
1.00000
1.00000
1.00000
1.00000
1.00000
1.00000
1.00000
1.00000
1.00000
1.00000
1.00000
1.00000
1.00000
1.00000
1.00000
1.00000
1.00000
1.00000
1.00000
1.00000
1.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000


And what happens when launched under matlab

Zitat:

expectedPopulationOfGroundState =

0.0000
0
0
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000


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

Thank you very much for your help.

Best regards,

Renaud
Private Nachricht senden Benutzer-Profile anzeigen
 
Harald
Forum-Meister

Forum-Meister


Beiträge: 24.434
Anmeldedatum: 26.03.09
Wohnort: Nähe München
Version: ab 2017b
     Beitrag Verfasst am: 28.02.2014, 08:56     Titel:
  Antworten mit Zitat      
Hi,

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? Wink I'm sure the folks from MathWorks could help you there Smile

Harald
Private Nachricht senden Benutzer-Profile anzeigen
 
rchretien
Themenstarter

Forum-Newbie

Forum-Newbie


Beiträge: 7
Anmeldedatum: 22.02.14
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 28.02.2014, 19:51     Titel:
  Antworten mit Zitat      
Hi Harald,

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

Anyway, I parallely asked the same question here : http://octave.1599824.n4.nabble.com.....d4662177i20.html#a4662512

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                                  # ??????????????

octave:23> mod(0.8, 0.2)
ans = 0


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.

Code:
function m = mod_ml(x,y)
  if fix(y) != y
   err      = abs (x./y - round(x./y)) < sqrt (eps);
   m       = mod (x,y);
   m(err) = 0;
  end
end


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.

Anyway, I thank you for your help and amability.

Renaud
Private Nachricht senden Benutzer-Profile anzeigen
 
Harald
Forum-Meister

Forum-Meister


Beiträge: 24.434
Anmeldedatum: 26.03.09
Wohnort: Nähe München
Version: ab 2017b
     Beitrag Verfasst am: 28.02.2014, 20:26     Titel:
  Antworten mit Zitat      
Hi,
Zitat:

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:
if mod(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):
Code:
0.6 - 0.2 - 0.2 - 0.2
ans =
  -5.5511e-17

The MathWorks people apparently were smart in implementing mod in a way that it would not be affected by this issue - at least not in your case.

Harald
Private Nachricht senden Benutzer-Profile anzeigen
 
rchretien
Themenstarter

Forum-Newbie

Forum-Newbie


Beiträge: 7
Anmeldedatum: 22.02.14
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 01.03.2014, 22:43     Titel:
  Antworten mit Zitat      
Hi,

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
Private Nachricht senden Benutzer-Profile anzeigen
 
Neues Thema eröffnen Neue Antwort erstellen



Einstellungen und Berechtigungen
Beiträge der letzten Zeit anzeigen:

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
.





 Impressum  | Nutzungsbedingungen  | Datenschutz | FAQ | goMatlab RSS Button RSS

Hosted by:


Copyright © 2007 - 2024 goMatlab.de | Dies ist keine offizielle Website der Firma The Mathworks

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.