Problem with inf value

I am having inf values, while I am supposed to be having 1.0396e+01 value... (which is apparently not at all an infinite value)...

I have 4 loops in my code : and I am trying to get the ratios condGaz(i,j,k) / condGaz(i,j,g), where k=1:3 ; g=1:3.

Here is a part of my code :

----------- Before this, I have 3 loops i (for Temperature), j (for Pressure), k=1:3 and this loop is my final loop ------------

for g=1:3
    if k~=g
    Phigk(i,j,k,g) = condGaz(i,j,k) / condGaz (i,j,g);
    else
    Phigk(i,j,k,g) =0  ;
    end
end

end 
end
end

Here are the condGaz(i,j,k) values.

And here are the values that I obtain for my ratio Phigk (i,j,k,g) = condGaz(i,j,k) / condGaz (i,j,g);

I get the impression that for g=1, I get all values right...

As soon as it gets to g=2 I get 1 value right and the other is INF

And for g=3, I get all the values INF....

I am not supposed to be having INF values because I manually did the ratio of for example Phigk(1,1,2,3) = condGaz(1,1,2) / condGaz(1,1,3) = 1.23e+00. (which is INF value if you see the results with the loops for Phigk(:,:,2,3)...

I tried format shortE and stuff... And it did not work... So I am confused about the sense of this INF value.... It is supposed to give INF value when the result is too big... But here we are talking values of 1.23e+00.... Can you please help me?

[Atsushi Ueno]

Possible story is that, 1. the indexing is not working as expected, 2. the denominator becomes equal to zero, and 3. the result of the division by zero becomes inf. Division by zero does not become the error. You will need to output the denominator just before the division and debug your code.

[Plamen Bonev]

Thank you Atsushi for your answer.

Actually the indexing does work because It works well for g=1 (k=1, k=2, k=3), and it gives the right results (I checked them manually). Actually when g=1 with k=1, 2, 3 it has to do :

• condGaz (: , : , 2)/ condGaz (:, : ,1)
• condGaz (: , : , 3) / condGaz (:, :, 1)

So it actually devides values of the following order of magnitude :

• e-02 / e-01
• e-02 / e-01

So it yields a finite value.

But when it comes to calculate for g=2 (k=1, k=2, k=3) it has to do :

• condGaz (: , : , 1)/ condGaz (:, : ,2)
• condGaz (: , : , 3) / condGaz (:, :, 2)

So it actually devides values of the following order of magnitude:

• e-01/e-02
• e-01/e-02

So here the denominator is smaller than the numerator by 10... But this should not be a reason for having INF values because this division simply yields a value of approximatively 1.023e+00....

The denominator actually never gets to close to zero to have INF values.... Manually Matlab does calculate a finite value when I am doing myself the ratios...

So I still do not understand...

[Walter Roberson]

I recommend that you debug with

dbstop if naninf

[Plamen Bonev]

Thank you Walter,

This does not help because it stops my code and so nothing happens...

Actually, I think my code is showing INF value because 1e+01 is bigger than the other finite values which are 1e-02... But I still does not want INF values because I could not do further calculations...

I will appreciate any further help with this, please....

[Atsushi Ueno]

>Actually the indexing does work because It works wel, and it gives the right results (I checked them manually).

Why is that a reason for the indexing to be correct?

MATLAB does produce the same calculation results whether on the command prompt or on the script.

Why don't you try debugging like below?

    if k~=g
    disp(['i,j,k,g = ' num2str([i j k g]) 'condGaz (i,j,g) = ' num2str(condGaz(i,j,g))]);
    Phigk(i,j,k,g) = condGaz(i,j,k) / condGaz (i,j,g);
    else

[Plamen Bonev]

Atsushi, thank you for your answer.

I displayed all and actually the condGaz (i,j,g) are all arounded to 0, when in reality they are of the order of magnitude of 1e-2... condGaz(i,j,g) is not 0 only when condGaz(i,j,g) is bigger than 1e-01...

For example condGaz(1,1,1) = 2.5935e-01 and so when I "display" it shows 0.25.

But condGaz(1,1,2) = 2.4946e-02 and so when I "display" it shows 0....

So it rounds my results to 0 that is why I am getting the INF values...

How can I debug this please?

I tried short g, short e, long e and stuff... I still get the same...

I appreciate your help.

[Plamen Bonev]

Here the "display" results and the "actual calculation results":

The calculation results for condGaz(i,j,g) are :

And the displayed results are :

So you can see that the results of 1e-02 are all rounded to 0... So inf values come from there...

I don't understand why my results are rounded to 0....

[Walter Roberson]

This does not help because it stops my code and so nothing happens...

It stops your code allowing you to find the very first step the problem occurs at, and you can then examine the values of the variables.

You do not show us your outer loops. We can predict that something in your loops is setting the values to 0, either explicitly (indexing) or else that you initialize the variable to zeros inside a loop when you should have initialized it before the loop. For example,

for K = 1 : 2
   G = zeros(5,7,3);
   G(:,:,K+1) = rand(5,7);
end

compared to

G = zeros(5,7,3);
for K  = 1  : 2
  G(:,:,K+1) = rand(5,7);
end

[Atsushi Ueno]

Thank you. So the indexing is not the cause, but rounding is. You must know the reason why rounding happens, but I do not. Anyway, the cause of the Inf problem has been solved.

[Plamen Bonev]

Walter, thank you again, I appreciate it.

Actually It does not send them to 0, they are just rounded to 0.

For example you can see that when my values condGaz (i, j, g) are of the order of magnitude of e-01, it displays the right value. For example for my condGaz (1, 1, 1)=0.259e-01 when I "display" it shows 0.25.

However, when my values are of the order of magnitude of e-02, my "displayed" answer is 0... So it is rounded to 0...

For example for my condGaz (1, 1, 2)= 2.4946e-02 i get from the "displayed" values a 0 result....

So for me it is a rounding problem.

I appreciate your help, guys.

[Plamen Bonev]

Thank you Atsushi for your collaboration.

But I still am stuck with the rounding problem, so if you have any ideas i will really appreciate it.

I have never had this problem, i have never seen my values being rounded...

[Atsushi Ueno]

I thought you know the reason why denominator becomes zero because you declare that it is rounding problem. However, you explained that rounding is happening because the debug output shows 0 when the value should be 2.4946e-02. I think that's probably your assumption. There is no specification like "floating point value less than e-02 order is rounded to be 0".

Why don't you consider more about @Walter Roberson's comment? Indexing issue is not only for reading value but also for storing value.

[Plamen Bonev]

Yes actually I did see my display number more carefully and it appears that indexing is wrong...

For example for k=1, it set all the values of condGaz(i,j,g) for g = 2, 3 to 0. it only displays condGaz(i,j,g) for g=1

For k=2, it set the value of condGaz(i,j,g) for g=3 to 0. It only displays condGaz(i,j,g) for g=1, 2

For k=3, it set non of the values to 0 and shows them all...

Here are the results:

[Plamen Bonev]

Here is the full version of my code :

format shortE
%Constante des gaz parfaits en J/mol.K
Rg=8.314;

%---------------------- Températures en K -------------------------

%T=[200 : 100 : 1000];
%T=[1000 : 100 : 2000];
T=[600 : 200 : 2000];

% -------------------- Pressions en Pa -----------------------------

%------------ P = [0.01 Mpa ; 0.1 Mpa] ---------------------
% P=[0.1e5:0.1e5:1e5]
% ----------- P = [0.1 MPa (1 bar) ; 10 MPa (10 bar)] ---------
% P=[1e5:2e5:100e5];

% ----------- P = [40MPa, 45MPa] -----------
P=[400e5:50e5:450e5];

% ----------- P = [0.1MPa, 1MPa] -----------
% P=[1e5:1e5:10e5];

% P=[0.1e5 : 1e5 : 200e5];

% ---------- Echelle log pour la Pression (de 1e4 Pa (0.01 MPa) à 1e7 Pa (10MPa) avec 1000 points entre les
% deux) ---------
% P=logspace(4, 7, 1000)

%He indice 1
%Kr indice 2
%Xe indice 3

% ------------------- Nombre de gaz considéré -----------------------
NbGaz=[1, 2, 3];

%--------------Paramètres critiques pour chaque gaz (pour le modèle avancé) ----------

%Masses molaires en kg/mol
M=[4.003e-3, 83.798e-3, 131.293e-3];

%Pressions critiques Pc en Pa
Pc=[0.2275e6, 5.51e6, 5.84e6];

%Températures critiques Tc en K
Tc=[5.2, 209.4, 289.7];

%Densités critiques rhoc en kg/m^3
rhoc=[69.64, 908.4, 1110];

%Coefficients directeurs pour viscosité dynamique 
Anu=[3.063e-7, 6.963e-7, 7.568e-7];

%Coefficients de température pour viscosité dynamique en K
Tnu=[-21.33, 71.07, 112.31];

%Exposants pour viscosité dynamique 
n=[0.724, 0.667, 0.655];

%--------------- Constantes a et b pour chaque gas (pour le modèle
%simplifié) ---------------

%He indice 1
%Kr indice 2
%Xe indice 3

aURGAP=[176.32e-5, 9.8e-5, 4.6e-5];
bURGAP=[0.77227, 0.8007, 0.8425];
aROCHE=[284e-5, 8.11e-5, 5.14e-5];
bROCHE=[0.7, 0.83, 0.83];
aBISON=[263.9e-5, 8.247e-5, 4.351e-5];
bBISON=[0.7085, 0.8363, 0.8616];

%------------ Initialisation des vecteurs/matrices -----------------
visc0=zeros(length(T),1);
cond0=zeros(length(T),1);
Vetoile=zeros(length(NbGaz),1);
condPseudoCr=zeros(length(NbGaz),1);
rho=zeros(length(T), length(P), length(NbGaz));
B2=zeros(length(T), length(NbGaz));
B3=zeros(length(T), length(NbGaz));
theta=zeros(length(T), length(NbGaz));
coef1=zeros(length(T), length(P));
coef2=zeros(length(T), length(P), length(NbGaz));
coef3=zeros(length(T), length(P), length(NbGaz));
racines=cell(length(T), length(P), length(NbGaz));
vraiRacines=cell(length(T), length(P), length(NbGaz));
B=cell(length(T), length(P), length(NbGaz));
imagi=cell(length(T), length(P), length(NbGaz));
rhor=zeros(length(T), length(P), length(NbGaz));
psi=zeros(length(T), length(P), length(NbGaz));
condGaz=zeros(length(T), length(P), length(NbGaz));
condGazURGAP=zeros(length(T), length(NbGaz));
condGazROCHE=zeros(length(T), length(NbGaz));
condGazBISON=zeros(length(T), length(NbGaz));
rapportURGAP=zeros(length(T), length(P), length(NbGaz));
pourcentageAugm=zeros(length(T), length(P), length(NbGaz));
diffURGAP=zeros(length(T), length(P), length(NbGaz));
diffBISON=zeros(length(T), length(P), length(NbGaz));
diffROCHE=zeros(length(T), length(P), length(NbGaz));
Psigk=zeros(length(T), length(P), length(NbGaz), length(NbGaz));
Phigk=zeros(length(T), length(P), length(NbGaz), length(NbGaz));


for k=1:numel(NbGaz)
    %Volumes molaires caractéristiques
    Vetoile(k)=Rg*Tc(k)/Pc(k);
    condPseudoCr(k)=(0.201e-4*((Tc(k))^(0.277))*((M(k))^(-0.465)))/((0.291*Vetoile(k))^(0.415));
for i=1:numel(T)
    %----------Conductivité de chaque gaz (modèle simplifié)---------
        %Données a et b URGAP
    condGazURGAP(i,k)=aURGAP(k)*(T(i))^(bURGAP(k));
        %Données a e b ROCHE
    condGazROCHE(i,k)=aROCHE(k)*(T(i))^(bROCHE(k));
        %Données a et b BISON
    condGazBISON(i,k)=aBISON(k)*(T(i))^(bBISON(k));
    
    %----------Modèle avancé---------
    theta(i,k)=(T(i))/(Tc(k));
    visc0(i,k)=Anu(k)*(T(i)-Tnu(k))^(n(k));
    cond0(i,k)=visc0(i,k)*((15*Rg)/(4*M(k)));
  if k==2 | k==3
    B2(i,k)=(-102.6+(102.732-0.001*theta(i,k)-0.44*((theta(i,k))^(-1.22)))*tanh(4.5*sqrt(theta(i,k))))*Vetoile(k);
    B3(i,k)=(0.0757+(-0.0862+(-3.6e-5)*theta(i,k)+0.0237*((theta(i,k))^(-0.059)))*tanh(0.84*theta(i,k)))*((Vetoile(k))^2);
  else 
    B2(i,k)=(8.4-0.0018*T(i)+(115/(sqrt(T(i))))-835/(T(i)))*(10^(-6));
    B3(i,k)=0;
  end
        
for j=1:numel(P)
    %coef1, 2 et 3 corréspondent aux coefficients devant rho
coef1(i,j,k)=(Rg*T(i))/(P(j));
coef2(i,j,k)=(B2(i,k))*((Rg*T(i))/(P(j)));
coef3(i,j,k)=(B3(i,k))*((Rg*T(i))/(P(j)));
    % On cherche les racines de l'équation cubique d'état
racines{i,j,k}=((roots([coef3(i,j,k), coef2(i,j,k), coef1(i,j,k), -1])));
    % On récupère uniquement les racines réelles 
vraiRacines{i,j,k}=racines{i,j,k}(imag(racines{i,j,k})==0);
    % La densité corréspond à la valeur maximale de chaque racine trouvée
    % pour chaque jeu de (P, T, gaz)
rho(i,j,k)=(max(vraiRacines{i,j,k}))*M(k);
    % Densité réduite
rhor(i,j,k)=(rho(i,j,k))/rhoc(k);
    % Psi
% psi(i,j,k)=0.645*rhor(i,j,k)+0.331*((rhor(i,j,k))^2)+0.0368*((rhor(i,j,k))^3)-0.0128*((rhor(i,j,k))^4);
psi(i,j,k)=0.645+0.331*((rhor(i,j,k))^1)+0.0368*((rhor(i,j,k))^2)-0.0128*((rhor(i,j,k))^3);

%     % Calcul final de la conductivité d'un gaz
condGaz(i,j,k)=cond0(i,k)+(condPseudoCr(k)*psi(i,j,k)*(rhor(i,j,k)));

    %RapportsURGAP
rapportURGAP(i,j,k)=(condGaz(i,j,k))/(condGazURGAP(i,k));

    %Calcul d'erreur entre le modèle avancé/modèle simple
        % Erreur URGAP / avancé (en %)
diffURGAP(i,j,k)=(condGazURGAP(i,k)-condGaz(i,j,k))*10^2;
        % Erreur BISON/ avancé (en %)
diffBISON(i,j,k)=(condGazBISON(i,k)-condGaz(i,j,k))*10^2;
        % Erreur ROCHE/ avancé (en %)
diffROCHE(i,j,k)=(condGazROCHE(i,k)-condGaz(i,j,k))*10^2;

% POUR GUILHERME : LE PROBLEME COMMENCE ICI.
for g=1:numel(NbGaz)
    disp(['i,j,k,g = ' num2str([i j k g]) ' ; condGaz (i,j,g) = ' num2str(condGaz(i,j,g))])
    Phigk(i,j,k,g)=(condGaz(i,j,k))/(condGaz(i,j,g));
end
end
end
end

Answer 1

[Walter Roberson]

condGaz(i,j,k)=cond0(i,k)+(condPseudoCr(k)*psi(i,j,k)*(rhor(i,j,k)));

That sets condGaz(I,j,:) according to the largest k value that has been seen so far. For example at k = 1, condGaz(i,j,1) will be set, and condGaz(i,j,2) will not have been set yet.

for g=1:numel(NbGaz)
    disp(['i,j,k,g = ' num2str([i j k g]) ' ; condGaz (i,j,g) = ' num2str(condGaz(i,j,g))])

but numel(NbGaz) is greater than k = 1, so in this inner loop within for k, as soon as g is greater than the current k, you index condGaz(i,j,:) elements that are still 0, and that gets you a division by 0.

Perhaps that loop needs to be postponed until after the other arrays have been fully built.

[Plamen Bonev]

Thank you Walter for this clear answer. I really appreciate your collaboration.

So how do you recon I proceed? How can I postpone this loop? Sorry for those silly questions...

[Walter Roberson]

There are multiple ways to proceed. One of them is

for g=1:min(k, numel(NbGaz))
    disp(['i,j,k,g = ' num2str([i j k g]) ' ; condGaz (i,j,g) = ' num2str(condGaz(i,j,g))])
    Phigk(i,j,k,g)=(condGaz(i,j,k))/(condGaz(i,j,g));
    Phigk(i,j,g,k) = condGaz(i,j,g) ./ condGaz(i,j,k);  %symmetry
end

Another way is to remove the for g loop entirely, and then after the entire for k loop, do

Phigk = condGaz ./ permute(condGaz, [1 2 4 3]);

with no nested loops needed (but does need R2016b or later for this particular version of the code.)

[Plamen Bonev]

Thank you, Walter, you solved it all for me.

Cheers.