problem here

% Note :- The code working perfect with single value of of P but %when I want to use multi values like 1:0.5:8 , its says (??? Error %using ==> mpower %Inputs must be a scalar and a square matrix.) help plz

clear all

 L = 1;

 P = 5.0 % <<<-- 

 T(L) = 270 ;%K

 R = 8.314e-3;% MPa m3 / mole.K

 Tc = [126.19 304.1 190.4 305.3 369.8 408.05 425.15 460.4 469.8...  507.6 540.2 568.7]';

 Tr = T(L)./Tc;

 Pc = [3.3978 7.38 4.6 4.8839 4.87 3.648 3.796 3.385326 3.36  3.02... 2.81 2.49]';

 omega = [0.0377 0.2236 0.0115 0.0995 0.1523 0.177 0.2002 0.2270... 

 0.2515 0.3013 0.3495 0.3996]';

 s = 0.480 + 1.574.*omega - 0.176.*omega.^2;

 alpha = (1 + s.*(1-sqrt(Tr))).^2;

 a = 0.42747.*(((R*Tc).^2)./Pc);

 b = 0.08664.*(R.*Tc./Pc);

 %-----------------------------------

 %For Gas phase , using y

 %--------------------------------------

 q=0;

 y = [0.05651 0.00284 0.833482 0.07526 0.02009 0.00305 0.0052  0.0012...

0.000144 0.00068 0.000138 0.00011]';

 x = [0.025 0.001 0.8819 0.078 0.01 0.000525 0.0019 0.0006  0.00062...
  0.00034 0.000067 0.000054]';

 while q==0

     % for Vapor Phase using y

     sumofA = 0;

     rou = zeros(12,1);

    for i = 1:12

        for j = 1:12

            sumofA = sumofA + 
y(i).*y(j)*sqrt(a(i).*a(j)*alpha(i).*alpha(j));

            rou(i) = rou(i) + 
y(j).*sqrt(a(i).*a(j).*alpha(i).*alpha(j));

        end

    end

    sumofB = sum(y.*b);

    %----------------------------------------

    % For Liquid Phase,  using x

    %---------------------------------------

    sumofAA = 0;

    mew = zeros(12,1);

    for i=1:12

        for j=1:12

            sumofAA = sumofAA + 
x(i).*x(j)*sqrt(a(i).*a(j)*alpha(i).*alpha(j));

            mew(i) = mew(i) + 
x(j).*sqrt(a(i).*a(j).*alpha(i).*alpha(j));

        end

    end

    sumofBB = sum(x.*b);

    % Calculate the Values of A and B for gas Phase

    AV = (sumofA .* P)./(R.^2.*T(L).^2);

    BV =  (sumofB.*P)./(R.*T(L));

    % Calculate the Values of A and B for Liquid Phase

    AL = (sumofAA .* P)./(R.^2.*T(L).^2);

    BL =  (sumofBB.*P)./(R.*T(L));

    % Find the compressibility factor for gas Phase

    ZV=roots([1 -1  AV-BV-BV^2 -AV*BV]);

    Zv=max(ZV);

    Zvc=isreal(Zv);

    if Zvc == false

        Zv=abs(Zv);

        disp('error 1')

    end

    % Find the compressibilty factor for Liquid

    ZL=roots([1 -1  AL-BL-BL^2 -AL*BL]);

    Zl=min(ZL);

    Zlc=isreal(Zl);

    if Zlc == false

        Zl=abs(Zl);

        disp('error 2')

    end

    %Fagucity of Liquid

    phil =exp((b.*(Zl-1)./sumofBB) - log(Zl-BL)-(AL/BL).*
((2*mew/sumofAA)-(b/sumofBB)).*log(1 + BL/Zl));

    %disp(phil)

    %Fagucity of Vapor

    phiv =exp((b.*(Zv-1)./sumofB) - log(Zv-BV)-(AV/BV).*
((2*rou/sumofA)-(b/sumofB)).*log(1 + BV/Zv));

    %disp(phiv);

    %  equilibirum calculation Ki

    K = (phil./phiv);

    % disp(K)

    raw(1)= sum(y./K)-1;

    if (abs(raw(1)) < 1e-5 & abs((y./K)-x) < 1e-5)

        % disp(K) 

        break

    else

        if L==1

            L=L+1;

            T(L) = 1.01 * T(L-1);

            raw(2)=raw(1);

        else

            L=L+1;

            T(L)=T(L-1)-raw(1)/(raw(1)-raw(2)).*(T(L-1)-T(L-2));

            raw(2)=raw(1);

        end

        %adjustment

        N = y./K;

        x = (1-0.5).*N + 0.5.*x;

    end

    n=0;

    n=n+1;

    if n>1000

        break

    end

end disp(T(L));

[Jan]

Please copy the full error message and mention the line, which causes the error.

Answer 1

[Jan]

At first your program fails in :

ZV = roots([1 -1  AV-BV-BV^2 -AV*BV]);

This can be fixed - I've inserted commas to improve the readability:

ZV = roots([1, -1,  AV-BV-BV.^2, -AV.*BV]);

You need the elementwise power ".^", not the matrix power "^".

But later on the script fails in:

phil = exp((b.*(Zl-1)./sumofBB), ...
     -log(Zl-BL)-(AL/BL) .* ...
     ((2*mew./sumofAA)-(b./sumofBB)) .* log(1 + BL/Zl));

I do not understand, what you want to calculate. This line contains vectors of length 12 and 15. Perhaps you want to use bsxfun to get a [12 x 15] matrix as result, but I cannot know this detail.

You can use the debugger to find such problems by your own:

dbstop if error

Then Matlab stops, when the error occurs and you can inspect the values of the variables in the command window.

[Nasir Qazi]

wht do u mean by 'bsxfun' and secondly how to use elementwise power?