fsolve problem

hi i am trying to run my programm and i get the following error:

    ??? Error using ==> interp1 at 165
  The interpolation points XI should be real.
Error in ==> airProp2 at 62
    varargout{1}=interp1(airTab(:,1),airTab(:,col),T);

Error in ==> opp_NTU at 31
Cph_pre =airProp2(Th_out_eva+273.15, 'cp');

Error in ==> trustnleqn at 199
    F = feval(funfcn{3},x,varargin{:});

Error in ==> fsolve at 378
    [x,FVAL,JACOB,EXITFLAG,OUTPUT,msgData]=...

and the code that i run :

function F=opp_NTU(x)

Th_out_eva=x(1);  Th_out_pre=x(2);  Q_pre=x(3);  A_pre=x(4);  Th_out_sup=x(5);
Q_eva=x(6);  A_eva=x(7);  Tc_out_sup=x(8);  Q_sup=x(9);  A_sup=x(10);
NTU_pre=x(11);
NTU_sup=x(12);

format long
A=18.24;%m^2
mh=0.49;%kg/s
mc=0.0482;%kg/s
P=33;%bar
Tc_in_pre=90.56;
Tc_out_pre=XSteam('Tsat_p',P);
Tc_in_eva=Tc_out_pre;
Tc_out_eva=Tc_in_eva;
Tc_in_sup=Tc_out_eva;
Th_in_sup=503.5;
U_pre=80.2589;  
U_eva=85.4191; 
U_sup=16.7160;
Cph_sup=airProp2(503.5+273.15, 'cp');
Cpc_pre=(1000*XSteam('Cp_pT',33,90.56));
Cpc_sup=(1000*XSteam('CpV_p',33));

%%%%%%%%%%%%%%%%%%%%%%%preaheater%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Cph_pre =airProp2(Th_out_eva+273.15, 'cp');

F(1) = Q_pre - (mh*Cph_pre*(Th_out_eva - Th_out_pre));
F(2) = Q_pre - (mc*Cpc_pre*(Tc_out_pre - Tc_in_pre));

Ch_pre=mh*Cph_pre;
Cc_pre=mc*Cpc_pre;

if Ch_pre<=Cc_pre
    Cmin_pre=Ch_pre;
    Cmax_pre=Cc_pre;
else
    Cmin_pre=Cc_pre;
    Cmax_pre=Ch_pre;
end

C_pre=Cmin_pre/Cmax_pre;
Qmax_pre=Cmin_pre*(Th_out_eva-Tc_in_pre);
%A_pre=
F(3)=NTU_pre-((U_pre*A_pre)/(Cmin_pre*1000));
e1_pre=exp(-C_pre*(NTU_pre^0.78))-1;
e2_pre=((NTU_pre^0.22)/C_pre)*e1_pre;
e_pre=1-exp(e2_pre);

F(4)= Q_pre -(e_pre*Qmax_pre);

%%%%%%%%%%%%%%%%%%%%%%%%%%evaporator%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

hc_out_eva=(1000*XSteam('hV_p',P));
hc_in_eva=(1000*XSteam('hL_p',P));
Cph_eva=airProp2(Th_out_sup+273.15, 'cp');

F(5) = Q_eva - (mh*Cph_eva*(Th_out_sup - Th_out_eva));
F(6) = Q_eva - (mc*(hc_out_eva - hc_in_eva));

Cpc_eva=(1000*XSteam('Cp_pT',33,Tc_out_pre));

DT1_eva=x(5) - Tc_out_eva;%%%%%%%Th_in_eva=th_out_sup
DT2_eva=x(1) - Tc_in_eva;
DTlm_eva=(DT1_eva-DT2_eva)/(log(DT1_eva/DT2_eva));

F(7) = Q_eva - (U_eva * A_eva * DTlm_eva);

%%%%%%%%%%%%%%%%%%%%%%%%%superheater%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
F(8) = Q_sup - (mh*Cph_sup*(Th_in_sup - Th_out_sup));
F(9) = Q_sup - (mc*Cpc_sup*(Tc_out_sup - Tc_in_sup));

Ch_sup=mh*Cph_sup;
Cc_sup=mc*Cpc_sup;

if Ch_sup<=Cc_sup
    Cmin_sup=Ch_sup;
    Cmax_sup=Cc_sup;
else
    Cmin_sup=Cc_sup;
    Cmax_sup=Ch_sup;
end

C_sup=Cmin_sup/Cmax_sup;
Qmax_sup=Cmin_sup*(Th_in_sup-Tc_in_sup);
%A_pre=
F(10)=NTU_sup-((U_sup*A_sup)/(Cmin_sup*1000));
e1_sup=exp(-C_sup*(NTU_sup^0.78))-1;
e2_sup=((NTU_sup^0.22)/C_sup)*e1_sup;
 e_sup=1-exp(e2_sup);

F(11)=Q_sup-(e_sup*Qmax_sup);
F(12) = A - A_pre - A_eva - A_sup;

% x0=[310;240;31410;3.4;450;84000;7.2;380;17000;8;0.6;0.6];

finnaly the problems seems to be something with the funtion that i call inside the program.i donwload the function from matlab and the code is :

function varargout=airProp2(T, prop)
  %----------------------------------------------
  % Interpolates thermodynamic air properties
  % Temp. range:  100-2500 K
  % According to Eckert & Drake, Analysis of Heat
  % and Mass Transfer, p. 780 
  %
  % Values are in SI-units:
  %
  %   col-#  prop.   units
  %   ------------------------
  %  1    T    K
  %  2    rho    kg/m^3
  %  3    cp    J/(kg K)
  %  4    my    kg/ms
  %  5    ny    m^2/s
  %  6    k    W/(m K)
  %  7    alpha  m^2/s
  %  8    Pr    -
  %
  % Example 1:  out=airProp2(296, 'ny')
  % Example 2:  
  %  [cp, ny]=airProp2([333 444],{'cp' 'ny'})
  %----------------------------------------------
  % (c)2004 by Stefan Billig
  %----------------------------------------------
  % Last Change:  04-Jun-2004
  %----------------------------------------------

    % check # of input arguments
    if ~isequal(nargin,2)
      error('airProp2 requires 2 input arguments!')
      return
      % check temperature request
    elseif find(T<100)  | ~isnumeric(T)
      error('Valid temperature range: 100 <= T[K] <= 2500')
      return
    end
    % get table
    load propTabAir2
    % if multi property request
    if iscell(prop)
      % scan along cells
      for idx=1:length(prop)
        % identify property column
        col=find(strcmp(propInfo,prop(idx)));
        if isempty(col)
          disp(['Property "' char(prop(idx)) '" not recognized!'])
        else
          % create output
          varargout{idx}=interp1(airTab(:,1),airTab(:,col),T);
        end
      end
    % single property request
    else
      % identify property column
      col=find(strcmp(propInfo,prop));
      if isempty(col)
        disp(['Property "' prop '" not recognized!'])
      else
        % create output
        varargout{1}=interp1(airTab(:,1),airTab(:,col),T);
      end
    end

any ideas?

Answer 1

[Matt Tearle]

Given that airProp2 doesn't appear to modify T, my guess is the problem is with the value you're passing in, which is Th_out_eva+273.15.

Th_out_eva appears to be the first element of the input x to opp_NTU. I don't see how that could become complex unless the output (ie the function value F) went complex.

To check this, add

if ~isreal(x)
 disp('x went complex!')
 x
end

to the beginning of opp_NTU. And/or a similar check for F at the end.

Also, try adding the following options to fsolve: set Display to 'iter' and FunValCheck to 'on'.

Answer 2

[a black]

general things that i observed are: that if i run only airprop2 for the Cph_pre i get a real number. Also if in the code i dsplay Cph_pre its also a real number. Moreover i dont think its a problem only for the particular Cph_pre i think that if i fix this one in the next part of the code when i call again airprop2 i will have the same problem. I did what u suggested:add at the begininng

if ~isreal(x)
 disp('x went complex!')
 x
end

and at the end

if ~isreal(F)
 disp('F went complex!')
 F
end

and i get first tha F went coplex and then that x went coplex i dont seem to understand why since the value of Cph_pre or the value of Th_out_eva goes complex since the starting value is a real number

Moreover when i get the values of F for columns one to twelve i get a complex number for column eleven which is NTU_pre , and this happens for all the repetations that fsolve does

F went complex!

F =

    1.0e+003 *
    Column 1
   -1.119046171479280                     
    Column 2
    0.364779214609764                     
    Column 3
    0.000458754150674                     
    Column 4
    0.421579040080665                     
    Column 5
    2.514320427960091                     
    Column 6
   -0.352760403182503                     
    Column 7
    1.599735803776945                     
    Column 8
   -3.107807138104990                     
    Column 9
   -2.287813142223574                     
    Column 10
    0.000133099687982                     
    Column 11
    2.095889001905194 - 0.003273726807226i
    Column 12
   -0.000273738036145

Answer 3

[a black]

Dear Matt, I think i found the error thanks to you. Since i had a problem with NTU_pre i go and checked there and i found the error. Thanx once again for your valuable help. For the record the error had something to do with the Cp that Xsteam gives and the Cp that airprop gives. The firsts are given in KJ/Kg*K and the others in J/Kg*K

Answer 4

[Walter Roberson]

NTU_pre and NTU_sup are both raised to fractional powers. If either one of them went negative (F(11) and F(12) respectively) then on the next cycle complex numbers would be created, and those complex values would probably propagate quickly.

Try checking for x(11) or x(12) being negative at the top of the cycle.

Note: whether a negative number raised to a fractional power has a real or complex result mathematically depends upon the representation of the fraction. For example, (x^2)^(1/4) = x^(2/4) calls for squaring the negative first and then calculating the principle square root of that result, which would be positive. But if you reduce 2/4 to 1/2, then x^(1/2) would be the square root of a negative number, which would be imaginary. In practice, to raise values to a fractional power, Matlab uses the log formulation, generating the (complex) log of a negative number and only ending up with a real result if exp() of the complex number times the power happens to end up as a real number -- something that is not especially likely for 0.78 and 0.22 .