Error for finding column index for iterative solution which uses bilinear interpolation - it works outside of being used for iterations, but not inside the iterative solver.

I am trying to complete an iterative solution to the idea gas equation. For this I need to complete two bilinear interpolations per iteration for the dependent variables in the governing equations. My script for the iterations are as follows:

mdot = 1; %kg/s - mass flow into vessel
V = 1*10^5; %m3 - volume of vessel
mmol = 0.0021588; %kg/mol - molar mass of H2
R = 8.314/mmol; %J/kg/K
%Run Bilinearinterpolator.m
he = 3.9445e+03*1000; %J/kg - specific enthalpy 
                    %of input H2 (condition inside of inlet)
cp0 = 14.3840*1000; %J/kg/K - Specific heat capacity
                    %of input H2 (condition inside of inlet)

%Numerical Method Parameters and data storage matrices

dt = 10000;
t1 = 0;
tu = 10^6;
n = round((tu-t1)/dt);
t = t1:dt:tu;
m = zeros(1,n+1);
T = zeros(1,n+1);
P = zeros(1,n+1);
Res = [t;m;T;P];
P1 = zeros(1,n+1);
P2  = zeros(1,n+1);
P1idx = zeros(1,n+1);
P2idx = zeros(1,n+1);
T1 = zeros(1,n+1);
T2 = zeros(1,n+1);
hf11 = zeros(1,n+1);
hf12 = zeros(1,n+1);
hf21 = zeros(1,n+1);
hf22 = zeros(1,n+1);
h = zeros(1,n+1);
cpf11 = zeros(1,n+1);
cpf12 = zeros(1,n+1);
cpf21 = zeros(1,n+1);
cpf22 = zeros(1,n+1);
cp = zeros(1,n+1);

%Initial Conditions

P(1) = 3*10^6; %Pa
T(1) = 298.15; %K
m(1) = P(1)*V/(R*T(1)); %kg
h(1) = he; %J/kg
cp(1) = cp0*1000 ; %J/kg/K

load('hn.mat')
load("cpn.mat")

for i=1:n
    m(i+1) = m(i) + dt*mdot;
    T(i+1) = T(i) + dt*mdot*(T(i)*V/mmol+he-he)/(m(i)*cp(i)-m(i)*V/mmol);
    P(i+1) = P(i) + dt*(mmol*m(i)*(T(i+1)-T(i))/dt+T(i)*mmol*m(i));
    
    P1(i+1) = 0.5*floor(P(i+1)/0.5);
    P2(i+1) = 0.5*ceil(P(i+1)/0.5);

    Pidx = [3;3.5;4;4.5;5 ;5.5; 6; 6.5; 7; 7.5; 8; 8.5; 9; 9.5; 10; 10.5; 11; 11.5; 12; 12.5; 13; 13.5; 14; 14.5; 15; 15.5; 16; 16.5; 17; 17.5; 18; 18.5; 19; 19.5; 20];

    P1idx(i+1) = find(Pidx==P1(i+1))+1;
    P2idx(i+1) = find(Pidx==P2(i+1))+1;

    T1(i+1) = 0.15+floor(T(i+1)-0.15);
    T2(i+1) = 0.15+ceil(T(i+1)-0.15);

    hf11(i+1) = hn{hn.K==T1(i+1),P1idx(i+1)};
    hf12(i+1) = hn{hn.K==T2(i+1),P1idx(i+1)};
    hf21(i+1) = hn{hn.K==T1(i+1),P2idx(i+1)};
    hf22(i+1) = hn{hn.K==T2(i+1),P2idx(i+1)};

    h(i+1) = 1/((P2(i+1)-P1(i+1))*(T2(i+1)-T1(i+1)))*[P2(i+1)-P(i+1) P(i+1)-P1(i+1)]*[hf11(i+1) hf12(i+1); hf21(i+1) hf22(i+1)]*[T2(i+1)-T(i+1);T(i+1)-T1(i+1)];

    cpf11(i+1) = cpn{cpn.K==T1(i+1),P1idx(i+1)};
    cpf12(i+1) = cpn{cpn.K==T2(i+1),P1idx(i+1)};
    cpf21(i+1) = cpn{cpn.K==T1(i+1),P2idx(i+1)};
    cpf22(i+1) = cpn{cpn.K==T2(i+1),P2idx(i+1)};

    cp(i+1) = 1/((P2(i+1)-P1(i+1))*(T2(i+1)-T1(i+1)))*[P2(i+1)-P(i+1) P(i+1)-P1(i+1)]*[cpf11(i+1) cpf12(i+1); cpf21(i+1) cpf22(i+1)]*[T2(i+1)-T(i+1);T(i+1)-T1(i+1)];
end

plot(t,m)
plot(t,T)
plot(t,P)

cpn.mat and hn.mat have been attached.

The script is able to run down to line 64 where there is an error:

>> thermalmodel
Unable to perform assignment because the left and right sides have a different number of elements.

Error in thermalmodel (line 64)
    P1idx(i+1) = find(Pidx==P1(i+1))+1;

However when running a non-iterative equivalent script for the bilinear interpolation for h, it works fine. Please see the non iterative script attached below:

Pi = 3.0+0.0000000000001;
Ti = 298.15+0.0000000000001;

Pi1 = 0.5*floor(Pi/0.5);
Pi2 = 0.5*ceil(Pi/0.5);

Pidx = [3;3.5;4;4.5;5 ;5.5; 6; 6.5; 7; 7.5; 8; 8.5; 9; 9.5; 10; 10.5; 11; 11.5; 12; 12.5; 13; 13.5; 14; 14.5; 15; 15.5; 16; 16.5; 17; 17.5; 18; 18.5; 19; 19.5; 20];

P1idx = find(Pidx==Pi1)+1;
P2idx = find(Pidx==Pi2)+1;

Ti1 = 0.15+floor(Ti-0.15);
Ti2 = 0.15+ceil(Ti-0.15);

load('hn.mat')

hf11 = hn{hn.K==Ti1,P1idx};
hf12 = hn{hn.K==Ti2,P1idx};
hf21 = hn{hn.K==Ti1,P2idx};
hf22 = hn{hn.K==Ti2,P2idx};

h = 1/((Pi2-Pi1)*(Ti2-Ti1))*[Pi2-Pi Pi-Pi1]*[hf11 hf12; hf21 hf22]*[Ti2-Ti;Ti-Ti1];
h
disp("kJ/kg")

load("cpn.mat")

cpf11 = cpn{cpn.K==Ti1,P1idx};
cpf12 = cpn{cpn.K==Ti2,P1idx};
cpf21 = cpn{cpn.K==Ti1,P2idx};
cpf22 = cpn{cpn.K==Ti2,P2idx};

cp = 1/((Pi2-Pi1)*(Ti2-Ti1))*[Pi2-Pi Pi-Pi1]*[cpf11 cpf12; cpf21 cpf22]*[Ti2-Ti;Ti-Ti1];
cp
disp("kJ/kg/K")

Please can someone identify why it works for the non-iterrative script and not the iterative script.

Kind regards and thanks for any help offered in advance

Answer 1

[Walter Roberson]

https://www.mathworks.com/matlabcentral/answers/1889267-wrong-output-using-find-while-filtering-a-matrix#answer_1142502

Do not compare floating point using ==