fsolve help

hi i am having a problem trying to use fsolve in order to have a solution for my problem. the problem is that although the programm runs for a single scenario when i try to use a 'for' in order to run 12 scenarios i get an error. My code is :

function F=opp_loop(x)
format long

Tc_in_pre=[90.56,90.48,90.33,90.25,90.52,90.44,90.17,90.09,90.44,90.29,90.09,90.05];
Th_in_sup=[503.5,459.3,394.2,363.1,483.7,423.3,353.6,334.5,435.2,382.4,317.6,291];
mh=[0.49,0.3947,0.3058,0.2054,0.5117,0.4408,0.3682,0.2457,0.4386,0.4095,0.3521,0.2690];  
mc=[0.0482,0.0328,0.0186,0.0106,0.0467,0.0304,0.0189,0.0119,0.0328,0.0234,0.0149,0.0096];  
P=[33,29,21,17,31,27,13,9,27,19,9,7];
u=sinartisi;
U_pre=u(:,1);
U_eva=u(:,2);
U_sup=u(:,3);

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(i));
Tc_in_eva=Tc_out_pre;
Tc_out_eva=Tc_in_eva;
Tc_in_sup=Tc_out_eva;

Cph_sup=airProp2(Th_in_sup(i)+273.15, 'cp');
Cpc_pre=(1000*XSteam('Cp_pT',P(i),Tc_in_pre(i)));
Cpc_sup=(1000*XSteam('CpV_p',P(i)));

Cph_pre =airProp2(x(1)+273.15, 'cp');

F(1) = x(3) - (mh(i)*Cph_pre*(x(1) - x(2)));
F(2) = x(3) - (mc(i)*Cpc_pre*(Tc_out_pre - Tc_in_pre(i)));

DT1_pre=x(1)- Tc_out_pre;
DT2_pre=x(2)- Tc_in_pre(i);
DTlm_pre=(DT1_pre-DT2_pre)/(log(DT1_pre/DT2_pre));

F(3) = x(3) - (U_pre(i) * x(4) * DTlm_pre);

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

F(4) = x(6) - (mh(i)*Cph_eva*(x(5) - x(1)));
F(5) = x(6) - (mc(i)*(hc_out_eva - hc_in_eva));

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(6) = x(6) - (U_eva(i) * x(7) * DTlm_eva);

F(7) = x(9) - (mh(i)*Cph_sup*(Th_in_sup(i) - x(5)));
F(8) = x(9) - (mc(i)*Cpc_sup*(x(8) - Tc_in_sup));

DT1_sup=x(5)-x(8);    %%%%%%%Th_in_eva=th_out_sup
DT2_sup=x(5)-Tc_in_sup;
DTlm_sup=(DT1_sup-DT2_sup)/(log(DT1_sup/DT2_sup));

F(9) = x(9) - (U_sup(i) * x(10) * DTlm_sup);
F(10) = A - x(4) - x(7) - x(10);

then in command window i type :

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

for i=1:12

options=optimset('Display','iter') ; 

[x] = fsolve(@opp_loop,x0) 

end

and i get the following error :

??? Subscript indices must either be real positive integers or logicals.

Error in ==> opp_loop at 22
Tc_out_pre=XSteam('Tsat_p',P(i));

Error in ==> fsolve at 254
            fuser = feval(funfcn{3},x,varargin{:});

Caused by:
    Failure in initial user-supplied objective function evaluation. FSOLVE
    cannot continue.

any help will be valuable

Answer 1

[Matt Tearle]

OK, the formatting is a bit hideous, but I think I see the problem: you're trying to reference the loop index i inside the function opp_loop. No can do. Function workspaces are local.

What you want to do is make opp_loop a function of two variables (x & i). Then to run:

x0=[310;240;31410;3.4;450;84000;7.2;380;17000;8];
x = zeros(10,12);
options=optimset('Display','iter') ;
for i=1:12
  f = @(x) opp_loop(x,i);
  x(:,i) = fsolve(f,x0);
end

(I'm assuming you want to save all the output x vectors, so I collected them in to a matrix, one column per x.) The point is that opp_loop is a function of x and i but f is an anonymous function handle of one variable ( x ), with the parameter i built into it.

Answer 2

[a black]

thanx for the help , i ve got one more thing. when i do what u tell me to do for my last scenario 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_provlimatiko_senario at 76
Cph_eva=airProp2(x(5)+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]=...

any ideas? airprop is a function that i call. i run the code without airprop and i get the same error .thanx in advance

[Matt Tearle]

It looks like T is going imaginary, but this is inside airProp2. So can you start a separate question about that, showing its code and the value of x(5)+273.15 that you call it with.