Code for iteration doesn't work

Hi everybody, I need support for Matlab iteration code. In particular, I'd like to calculate "lambda_v_ratio" parameter in this code:

     m=1.5; %politropic exponent
     pc_ad=120000; %intake pressure
     z_ad=5000; %%quota di adattamento
     eta_c=0.8; %compressor efficiency 
     eta_mc=0.9; %compressor mechanical efficiency
     eta_mt=0.9; %turbine mechanical efficiency
     Lc = cpa*T_ad/eta_c*((pc_ad/p_ad)^((k-1)/k)-1); %compressor work

     for j=1:length(z)
       if  z(j)<=z_ad 
           Tc(j) = T(j)+(Lc/cpa); %intake temperature
           pc(j)=pc_ad; %intake pressure
       else
           Tc(j) = T(j)+(Lc/cpa);
           beta(j) = (1+((Lc*eta_c./(cpa.*T(j))))).^(k/(k-1)); %pressure ratio
           pc(j)=beta(j)*p(j);
       end
       rho_c(j) = pc(j)/(R*Tc(j)); %intake density
       mu(j) = (pc_ad/p0)*((T0./Tc(j)).^0.5); %correct density
       lambda_v(j) = lambda_v0* ((T(j)/T0).^0.5); %altitude volumetric efficiency w/o turbocompressor

       lambda_v_ratio(j)=0.8;
       lambda_v_tc(j) = lambda_v(j)*lambda_v_ratio(j); %volumetric efficiency due to turbocharging
       air_flow_rate(j) = lambda_v_tc(j)*rho_c(j)*iV*n/(2*60);
       air_mass(j) = lambda_v_tc(j)*rho_c(j)*V; %per cylinder
       fuel_mass(j) = air_mass(j)/alfa;
       v1(j)= (V+V0)/air_mass(j);
       v2(j)=v1(j)/r;
       p2(j)=pc(j)*(r^k);
       T2(j) = p2(j)*v2(j)/R;
       T3(j) = (T2(j)*alfa*cpa+eta_b*Hi*1000000)/(cpg*(1+alfa));
       v3(j) = V0/(air_mass(j)+fuel_mass(j));
       p3(j) = R*T3(j)/v3(j);
       v4(j) = (V+V0)/(air_mass(j)+fuel_mass(j));
       p4(j) = p3(j)*((v3(j)/v4(j))^k2);
       T4(j) = v4(j)*p4(j)/R;
       ps(j)= p4(j)/4;
       ps2(j)=((-lambda_v_ratio(j)+1)*m*pc(j)*(r-1))+pc(j);

I'd like to calculate "lambda_v_ratio(j)" so that ps(j) and ps2(j) are almost equal (ps(j) and ps2(j) are both functions of lambda_v_ratio(j) and they have to converge). I tried with this code:

while (abs(ps(j)-ps2(j)))<100
      lambda_v_ratio(j)=lambda_v_ratio(j)+0.001

but it didn't work. Start value of lambda_v_ratio(j) is 0.8. Any suggestion? Thanks for your support.

[Walter Roberson]

That code phrase would continue until the difference became more than 100, which would be the opposite of wanting them "almost equal".

[Antonio Tricarico]

Thanks, but actually I've tried with >100 and the calculation never stops...it reaches very high values for lambda_v_ratio, whereas it should converge near values of 1 (I've already calculated in excel)

[Walter Roberson]

It would probably make more sense to put most of the code into a function and then fzero on the difference between the two values.

Answer 1

[Antonio Tricarico]

I've tried with this *.m function:

function f=usingfzero1_fun(lambda_v_ratio)

global m pc R p2 lambda_v V V0 rho_c r alfa cpa eta_b Hi cpg k2

f=(((-lambda_v_ratio+1)*m*pc*(r-1))+pc)-(((((R*(((((p2*(((((V+V0)/(((lambda_v*lambda_v_ratio*rho_c*V))))/r))/R))*alfa*cpa+eta_b*Hi*1000000)/(cpg*(1+alfa))))/((V0/(((lambda_v*lambda_v_ratio)*rho_c*V)+((((lambda_v*lambda_v_ratio)*rho_c*V)/alfa))))))))*((((((V0/(((lambda_v*lambda_v_ratio)*rho_c*V)+((((lambda_v*lambda_v_ratio)*rho_c*V)/alfa)))))))/((((V+V0)/((((lambda_v*lambda_v_ratio)*rho_c*V))+((((lambda_v*lambda_v_ratio)*rho_c*V)/alfa)))))))^k2))/4))

The main *.m file recalls this function through code below:

c0=0.5;

[c,fval,exitflag,output]=fzero('usingfzero1_fun',c0)

but it calculates only one value of c (i.e. lambda_v_ratio).

I need a vector of n-solutions, because I'm working inside a for cycle with n indices.

Could you take me on the right way?

Thank you.

[Walter Roberson]

We need to see the values for m pc R p2 lambda_v V V0 rho_c r alfa cpa eta_b Hi cpg k2 to test with. But we recommend against global. http://www.mathworks.com/help/matlab/math/parameterizing-functions.html

For your vector of n solutions, which variable will take on different values?

[Antonio Tricarico]

m, R, V, V0, r, alfa, cpa, eta_b, Hi, cpg and k2 are fixed values. pc, p2, lambda_v and rho_c have different values in each of the n cases, because they vary with the index.

[Antonio Tricarico]

The example b = 2; c = 3.5; cubicpoly = @(x) x^3 + b*x + c; x = fzero(cubicpoly,0) can be solved if b and c are vectors? That's my case.

[Walter Roberson]

fzero() cannot be used with vector-valued functions. fsolve() can be used with vector-valued functions, in which case it would attempt to solve all of the equations simultaneously. In situations where you have an equation to be solved over a range of values, you would iterate fzero using one value each time.