Iterating inside an iteration loop, while loop

I have two iteration loops and im untimtly trying to solve for e. The outer loop contain an inner iteration loop. The outer loop is dependent on the results form the inner loop. The inner loop should give a rhoj_2 value when conditions in the second while loop are meet:

while abs(h2_i-h2_j)>=10^(-4) & (h2_i>=h2_j) 

Then the outer loop should update e based using the new rhoj_2 until e converges via the first while loop

while abs((e(i+1)-e(i))/e(i)) >= 10^(-4) 

The inner loop will converge but the outer loop doesn't seem to update e with the new rhoJ_2 value for each iteration . Any help getting this outer iteration to converge/ update would be very much appreaicted.

Code:

clear;clc 
% KNOWN VARIABLES
rho1=5.3811*10^(-5); %kg/m3
P1= 3.2797; % N/m2---Do I CONVERT THIS??
T1=212.3; % K
v1=9385.0; %m/s
 
% Assume CPG before shock wave
cp=1.00; %kJ/kg
h1=cp*T1; 

%Set Up Outer Loop
i=1;
% e=[];
e(1)=0.08; % inital guess will change with each iteration
e(2)=0.1; % Need inital guess for e --error funcion set up

while abs((e(i+1)-e(i))/e(i)) >= 10^(-4) 
% Solve P and h2
P2_i(i)=P1+rho1*v1^2*(1-e(i));
h2_i(i)=h1+v1^2/2*(1-e(i)^2);

% INNER LOOP-iterate rho_2j start off with rho_2i as inital guess
% Read in curve fit data

%Brut Force curve fit data-CHECK WHY WOULDN'T READ IN DATA
h_data=[1.4000 0 0 0 0 0 0 0 0 0 0 1;
    1.4260 9.1800e-04 -0.0922 -0.0022 0.0198 -0.0366 -0.0775 0.0439 -15 -1 -1.0400 2;
    1.6469 -0.0621 -0.3350 0.0636 -0.0383 -0.0145 0.0734 -0.0024 -15 -1 -1.3600 3;
    1.4856 -0.4536 0.1521 0.3034 -0.4593 0.4484 0.2205 -0.2923 -10 -1 -1.6000 4;
    1.4000 0 0 0 0 0 0 0 0 0 0 5;
    1.4218 -3.6600e-04 -0.0836 6.7500e-04 0.0053 -0.1159 -0.0074 0.1462 -20 -1 -0.8600 6;
    1.7444 -0.0354 -0.4150 0.0619 0.0185 0.0436 0.0444 -0.0498 -20 -1.0400 -1.3360 7;
    1.4967 -0.0216 -0.1970 0.0309 -0.1577 -0.0092 0.1232 -0.0066 -10 -1.0500 -1.8950 8;
    1.1042 -0.0337 0.0318 0.0243 -0.1768 -0.0175 0.0804 0.0025 -15 -1.0800 -2.6500 9;
    1.4000 0 0 0 0 0 0 0 0 0 0 10
    1.4700 0.0079 -0.2442 -0.0256 0.8722 0.0495 -0.7642 1.4700e-04 -20 -1 -0.7420 11;
    3.1865 0.1379 -1.8953 -0.1035 -2.1457 -0.2727 2.0659 0.2230 -15 -1 -1.0410 12;
    1.6396 -0.0010 -0.3035 0.0165 -0.8522 -0.1012 0.5031 0.0436 -10 -1 -1.5440 13;
    1.5589 0.0559 -0.2118 -0.0235 -0.5490 -0.1018 0.2767 0.0460 -15 -1 -2.2500 14];


   

% Inital guess:
k=1;
%h2_j(1)=h1; % Inital guess for final enthalpy-should be less than inital?
h2_j(1)=h2_i-200; % Test will this converge-Inital guess for h2
rho2j(1)=rho1\e(1); % Inital guess for density
   
    
    while abs(h2_i-h2_j)>=10^(-4) & (h2_i>=h2_j) % Too Big--Doesnt converge-increase tol later

        % Curve Fit Constants
        Y=log10(rho2j(k)/1.292); % Changes with each iteration beacuse rho2j
        X=log10(P2_i/1.013^10^5);
        Z=X-Y;
        
         % This loop selects what constant to use base on X,Y,Z

         if Y>-0.5
             
             if Z<=0.3                 
                 c=h_data(1,1:11);
                 
             elseif (0.3 < Z)&&(Z<=1.15)                 
                 c=h_data(2,1:11);
                 
             elseif (1.15 < Z)&&(Z<= 1.60)                 
                 c=h_data(3,1:11);
                 
             elseif Z > 1.60                 
                 c=h_data(4,1:11);
                 
             end
             
             
             
         elseif (-4.5<Y)&&(Y<=-0.5)
             
             if Z <= 0.30                 
                 c=h_data(5,1:11);
                 
             elseif (0.30 < Z)&&(Z <= 0.98)                
                 c=h_data(6,1:11);
                 
             elseif (0.98 < Z)&&(Z <= 1.38)                 
                 c=h_data(7,1:11)
                 
             elseif (1.38 < Z)&&(Z <= 2.04)                 
                 c=h_data(8,1:11)
                 
             elseif Z > 2.04                 
                 c=h_data(9,1:11);
                 
             end
             
             
             
         elseif (-7 < Y)&&(Y< -4.5)
             
             if Z <= 0.398                 
                 c=h_data(10,1:11);
                 
             elseif (0.398 < Z)&&(Z <= 0.87)                 
                 c=h_data(11,1:11);
                 
             elseif (0.87 < Z)&&(Z <= 1.27)                 
                 c=h_data(12,1:11);
                 
             elseif (1.27 < Z)&&(Z <= 1.863)                 
                 c=h_data(13,1:11);
                 
             elseif  Z > 1.863                 
                 c=h_data(14,1:11);
                 
             end
             
                                         
         end
   
        % Enthaply
        g(k)=c(1)+c(2).*Y+c(3).*Z+c(4).*Y.*Z+(c(5)+c(6).*Y+c(7).*Z+c(8).*Y.*Z)./(1+exp(c(9).*(X+c(10).*Y+c(11))));
        % New enthapply calculated
        h2_j(k)=P2_i/rho2j(k)*(g(k)/(g(k)-1));
        % Does this converge if not then increase or decrease desnsity( Iterate)
       rho2j(k+1)=rho2j(k)-1e-8;
         k=k+1;
        

    end
e(i+2)=rho1/(rho2j(end));
i=i+1;
 
end

Answer 1

[Jan]

This is not a bug, but useful for efficient code: 10^(-4) is an expensive power operation while 1e-4 is a cheap constant.

if Z<=0.3                 
    c=h_data(1,1:11);
elseif (0.3 < Z)&&(Z<=1.15)
    %      ^^^^^^^^^^^^

This can be omitted. You have excluded z<= 0.3 before so the elseif can be reached only if Z > 0.3. This occurs mutliple times and you can simplify your code.

The code is easier to read, if you insert spaces around operators and use the elementwise multiplication .* only, if it is needed, not in general.

Is this a typo:

X=log10(P2_i/1.013^10^5);
%                 ^
% Do you mean:
X = log10(P2_i / 1.013e5);    

A simplified cleaned version (compare the readability):

% Inital guess:
k=1;
%h2_j(1)=h1; % Inital guess for final enthalpy-should be less than inital?
h2_j(1)  = h2_i - 200;  % Test will this converge-Inital guess for h2
rho2j(1) = rho1 \ e(1); % Inital guess for density
    
while abs(h2_i - h2_j) >= 10e-4 && (h2_i>=h2_j) % Too Big--Doesnt converge-increase tol later
    % Curve Fit Constants
    Y = log10(rho2j(k) / 1.292); % Changes with each iteration beacuse rho2j
    X = log10(P2_i / 1.013e5);
    Z = X - Y;
    
    % This loop selects what constant to use base on X,Y,Z
    if Y>-0.5   
         if Z <= 0.3
             d = 1;
         elseif Z <= 1.15                 
             d = 2;
         elseif Z <= 1.60                 
             d = 3;
         else                 
             d = 4
         end
         
     elseif -4.5 < Y
         if Z <= 0.3                 
             d = 5;
         elseif Z <= 0.98
             d = 6;
         elseif Z <= 1.38                 
             d = 7;
         elseif Z <= 2.04                 
             d = 8;
         else                
             d = 9;
         end
         
     elseif -7 < Y
         if Z <= 0.398                 
             d = 10;
         elseif Z <= 0.87                 
             d = 11;
         elseif Z <= 1.27                 
             d = 12;
         elseif Z <= 1.863                 
             d = 13;
         else                 
             d = 14;
         end
     end
     c = h_data(d, 1:11);
     
     % Enthaply
     g(k) = c(1) + c(2) * Y + c(3) * Z + c(4) * Y .* Z + ...
            (c(5) + c(6) * Y + c(7) * Z + c(8) * Y .* Z) ./ ...
            (1 + exp(c(9) * (X + c(10) * Y + c(11))));
     
     % New enthapply calculated
     h2_j(k) = P2_i / rho2j(k) * g(k) / (g(k)-1);
     % Does this converge if not then increase or decrease desnsity( Iterate)
     rho2j(k+1) = rho2j(k) - 1e-8;
     k = k+1;  
   end

   e(i+2) = rho1 / (rho2j(end));
   i      = i + 1;
end

[Cierra Keeling]

e still wont converge even if the tolerance is adjusted. e doesnt appear to be updating each iteration. Is the code getting cought in the inner loop? Any suggestions on how I can get e to converge

[Cierra Keeling]

Ive update my code but still am not getting a solution for e. I think the issue is the outer loop solving for e. When I put

% Recalculate e from new density from inner loop
    e(i+2)=rho1/rho2j(k);
    % Iterate
    i=i+1;

inside the inner loop where :

     while abs(h2_i - h2_j) >= 10e-4 & (h2_i>=h2_j)

e get updated for every iteration but the inner loop wont converge an thus the outer loop won't either. When I place the line in the outer loop where

if  abs((e(i+1)-e(i))/e(1))>=tolerance

e wont update for every iteration and only contains two values (the inital condtions)and thus the outer loop won't converge.

Any help on getting e to converge would be much appreaicted: Heres my update code with the code with e inside the outer loop.

clear;clc
% KNOWN VARIABLES
rho1=5.3811*10^(-5); %kg/m3
P1= 3.2797; % N/m2---Do I CONVERT THIS??
T1=212.3; % K
v1=9385.0; %m/s

% Assume CPG before shock wave
cp=1000; %J/kg
h1=cp*T1;

i=1
e(1)=0.08
e(2)=0.1
tolerance=1e-2 % Increase later

if  abs((e(i+1)-e(i))/e(1))>=tolerance
    
    P2_i(i)=P1+rho1* v1^2* (1-e(i));
    h2_i(i)=h1+ (v1^2)/2* (1-e(i)^2);
    
    % INNER LOOP-iterate rho_2j start off with rho_2i as inital guess
    % Read in curve fit data
    %Brut Force curve fit data-CHECK WHY WOULDN'T READ IN DATA
    h_data=[1.4000 0 0 0 0 0 0 0 0 0 0 1;
        1.4260 9.1800e-04 -0.0922 -0.0022 0.0198 -0.0366 -0.0775 0.0439 -15 -1 -1.0400 2;
        1.6469 -0.0621 -0.3350 0.0636 -0.0383 -0.0145 0.0734 -0.0024 -15 -1 -1.3600 3;
        1.4856 -0.4536 0.1521 0.3034 -0.4593 0.4484 0.2205 -0.2923 -10 -1 -1.6000 4;
        1.4000 0 0 0 0 0 0 0 0 0 0 5;
        1.4218 -3.6600e-04 -0.0836 6.7500e-04 0.0053 -0.1159 -0.0074 0.1462 -20 -1 -0.8600 6;
        1.7444 -0.0354 -0.4150 0.0619 0.0185 0.0436 0.0444 -0.0498 -20 -1.0400 -1.3360 7;
        1.4967 -0.0216 -0.1970 0.0309 -0.1577 -0.0092 0.1232 -0.0066 -10 -1.0500 -1.8950 8;
        1.1042 -0.0337 0.0318 0.0243 -0.1768 -0.0175 0.0804 0.0025 -15 -1.0800 -2.6500 9;
        1.4000 0 0 0 0 0 0 0 0 0 0 10
        1.4700 0.0079 -0.2442 -0.0256 0.8722 0.0495 -0.7642 1.4700e-04 -20 -1 -0.7420 11;
        3.1865 0.1379 -1.8953 -0.1035 -2.1457 -0.2727 2.0659 0.2230 -15 -1 -1.0410 12;
        1.6396 -0.0010 -0.3035 0.0165 -0.8522 -0.1012 0.5031 0.0436 -10 -1 -1.5440 13;
        1.5589 0.0559 -0.2118 -0.0235 -0.5490 -0.1018 0.2767 0.0460 -15 -1 -2.2500 14];
    
    %Insert Inner Loop
    
    % Inital guess:
    k=1;
    h2_j(1)=h1; % Inital guess for final enthalpy-should be less than inital?
    % h2_j(1)  = h2_i - 200;  % Test will this converge-Inital guess for h2
    rho2j(1) = rho1 \ e(1); % Inital guess for density
    
    while abs(h2_i - h2_j) >= 10e-4 & (h2_i>=h2_j) % Too Big--Doesnt converge-increase tol later
        % Curve Fit Constants
        Y = log10(rho2j(k) / 1.292); % Changes with each iteration beacuse rho2j
        X = log10(P2_i / 1.013e5);
        Z = X - Y;
        
        % This loop selects what constant to use base on X,Y,Z
        if Y>-0.5
            if Z <= 0.3
                d = 1;
            elseif Z <= 1.15
                d = 2;
            elseif Z <= 1.60
                d = 3;
            else
                d = 4
            end
            
        elseif -4.5 < Y
            if Z <= 0.3
                d = 5;
            elseif Z <= 0.98
                d = 6;
            elseif Z <= 1.38
                d = 7;
            elseif Z <= 2.04
                d = 8;
            else
                d = 9;
            end
            
        elseif -7 < Y
            if Z <= 0.398
                d = 10;
            elseif Z <= 0.87
                d = 11;
            elseif Z <= 1.27
                d = 12;
            elseif Z <= 1.863
                d = 13;
            else
                d = 14;
            end
        end
        c = h_data(d, 1:11);
        
        % Enthaply
        g(k) = c(1) + c(2) * Y + c(3) * Z + c(4) * Y .* Z + ...
            (c(5) + c(6) * Y + c(7) * Z + c(8) * Y .* Z) ./ ...
            (1 + exp(c(9) * (X + c(10) * Y + c(11))));
        
        % New enthapply calculated
        h2_j(k) = P2_i / rho2j(k) * g(k) / (g(k)-1);
        % Does this converge if not then increase or decrease desnsity( Iterate)
        rho2j(k+1) = rho2j(k) - 1; %1e-8;
        k = k+1;
        
    end
    
else
    % Recalculate e from new density from inner loop
    e(i+2)=rho1/rho2j;
    % Iterate
    i=i+1;
    disp(e(end))
end

[Johan]

I suggest you check your update condition:

rho2j(k+1) = rho2j(k) - 1;

this update quickly leads to a negative value for rho2j; your comment says this represents a density so this may not be the correct way to update your loop ?

For example updating rho2j by 99% of the previous value converges. Up to you to check what makes sense for your calculation.

 clear;clc
% KNOWN VARIABLES
rho1=5.3811*10^(-5); %kg/m3
P1= 3.2797; % N/m2---Do I CONVERT THIS??
T1=212.3; % K
v1=9385.0; %m/s

% Assume CPG before shock wave
cp=1000; %J/kg
h1=cp*T1;

i=1;
e(1)=0.08;
e(2)=0.1;
tolerance=1e-2; % Increase later

%this is not updated by the loop, can be set once outside
    % Read in curve fit data
    %Brut Force curve fit data-CHECK WHY WOULDN'T READ IN DATA
    h_data=[1.4000 0 0 0 0 0 0 0 0 0 0 1;
        1.4260 9.1800e-04 -0.0922 -0.0022 0.0198 -0.0366 -0.0775 0.0439 -15 -1 -1.0400 2;
        1.6469 -0.0621 -0.3350 0.0636 -0.0383 -0.0145 0.0734 -0.0024 -15 -1 -1.3600 3;
        1.4856 -0.4536 0.1521 0.3034 -0.4593 0.4484 0.2205 -0.2923 -10 -1 -1.6000 4;
        1.4000 0 0 0 0 0 0 0 0 0 0 5;
        1.4218 -3.6600e-04 -0.0836 6.7500e-04 0.0053 -0.1159 -0.0074 0.1462 -20 -1 -0.8600 6;
        1.7444 -0.0354 -0.4150 0.0619 0.0185 0.0436 0.0444 -0.0498 -20 -1.0400 -1.3360 7;
        1.4967 -0.0216 -0.1970 0.0309 -0.1577 -0.0092 0.1232 -0.0066 -10 -1.0500 -1.8950 8;
        1.1042 -0.0337 0.0318 0.0243 -0.1768 -0.0175 0.0804 0.0025 -15 -1.0800 -2.6500 9;
        1.4000 0 0 0 0 0 0 0 0 0 0 10
        1.4700 0.0079 -0.2442 -0.0256 0.8722 0.0495 -0.7642 1.4700e-04 -20 -1 -0.7420 11;
        3.1865 0.1379 -1.8953 -0.1035 -2.1457 -0.2727 2.0659 0.2230 -15 -1 -1.0410 12;
        1.6396 -0.0010 -0.3035 0.0165 -0.8522 -0.1012 0.5031 0.0436 -10 -1 -1.5440 13;
        1.5589 0.0559 -0.2118 -0.0235 -0.5490 -0.1018 0.2767 0.0460 -15 -1 -2.2500 14];
    
while  abs((e(i+1)-e(i))/e(1))>=tolerance
    
    P2_i(i)=P1+rho1* v1^2* (1-e(i));
    h2_i(i)=h1+ (v1^2)/2* (1-e(i)^2);
    
    % INNER LOOP-iterate rho_2j start off with rho_2i as inital guess

    % Inital guess:
    k=1;
    h2_j(1)=h1; % Inital guess for final enthalpy-should be less than inital?
    % h2_j(1)  = h2_i - 200;  % Test will this converge-Inital guess for h2
    rho2j(1) = rho1 \ e(1); % Inital guess for density
    
    while abs(h2_i(i) - h2_j(k)) >= 10e-4 & (h2_i(i)>=h2_j(k)) % Too Big--Doesnt converge-increase tol later
        % Curve Fit Constants
        Y = log10(rho2j(k) / 1.292); % Changes with each iteration beacuse rho2j
        X = log10(P2_i(i) / 1.013e5);
        Z = X - Y;
        
        % This loop selects what constant to use base on X,Y,Z
        if Y>-0.5
            if Z <= 0.3
                d = 1;
            elseif Z <= 1.15
                d = 2;
            elseif Z <= 1.60
                d = 3;
            else
                d = 4
            end
            
        elseif -4.5 < Y
            if Z <= 0.3
                d = 5;
            elseif Z <= 0.98
                d = 6;
            elseif Z <= 1.38
                d = 7;
            elseif Z <= 2.04
                d = 8;
            else
                d = 9;
            end
            
        elseif -7 < Y
            if Z <= 0.398
                d = 10;
            elseif Z <= 0.87
                d = 11;
            elseif Z <= 1.27
                d = 12;
            elseif Z <= 1.863
                d = 13;
            else
                d = 14;
            end
        end
        c = h_data(d, 1:11);
        
        % Enthaply
        g(k+1) = c(1) + c(2) * Y + c(3) * Z + c(4) * Y .* Z + ...
            (c(5) + c(6) * Y + c(7) * Z + c(8) * Y .* Z) ./ ...
            (1 + exp(c(9) * (X + c(10) * Y + c(11))));
        
        % New enthapply calculated
        h2_j(k+1) = P2_i(i) / rho2j(k) * g(k) / (g(k)-1);
        % Does this converge if not then increase or decrease desnsity( Iterate)
        rho2j(k+1) = rho2j(k)*0.99; %1e-8;
        k = k+1;
        if k>1e4 % added escape condition to prevent infinite loop
            disp('inner loop break')
            break
        end
    end
    

    % Recalculate e from new density from inner loop
    e(i+2)=rho1/rho2j(end);
    % Iterate
    i=i+1;
%     disp(e(end))
    if i > 1e3 % added escape condition to prevent infinite loop
        disp('outer loop break')
        break
    end
end
plot(e)

k = 1423

i = 3