Thread Subject:

Positive real root

I have vector with complex and real roots - P. I have to write code to find only positive real root(in this calculation its only one). This is my code:
R1=1.230; X1=1.500; Xm=53.700; R2=0.787; X2=2.490;
Vn=220; Tn=36.340;T=10; Pn=5.500; I1n=11.500; I0n=19.5; f=30;
fn=50; P=4; sn=0.036; ns=1500; P1n=6451.5;E1n=209; m=T/Tn; alfa=f/fn;
ws=157; nn=1446; wn=151.348; J=0.017; a=0.0327; b=3.112; k=1.400;
Pag=Tn*ws; Pmag=P1n-Pag-(3*R1*I1n^2+0.005*Pn);
B=3*R1*I0n^2; C=((R1+R2)*Pag^2)/(3*E1n^2); D=Pmag;

a0=1;
a1=2*b/(1-b);
a2=(b/(b-1))^2+(2*B*(b-a)-C*m^2*(b-1)^2)/(D*alfa^k*(b-1)^2);
a3=(2*C*m^2*b)/(D*alfa^k*(b-1));
a4=-(C*m^2*b)/(D*alfa^k*(b-1)^2);

Pol=[a0 0 a1 0 a2 0 a3 0 a4];
P= roots(Pol)

Here the positive real root is 0.2968. But I don't know how to do this by code ( not to see all roots and to say - This is real positive root). Please help me.

"Kiril " <kkirqkov@gmail.com> wrote in message <gs7uda$fkc$1@fred.mathworks.com>...
> I have vector with complex and real roots - P. I have to write code to find only positive real root(in this calculation its only one). This is my code:
> R1=1.230; X1=1.500; Xm=53.700; R2=0.787; X2=2.490;
> Vn=220; Tn=36.340;T=10; Pn=5.500; I1n=11.500; I0n=19.5; f=30;
> fn=50; P=4; sn=0.036; ns=1500; P1n=6451.5;E1n=209; m=T/Tn; alfa=f/fn;
> ws=157; nn=1446; wn=151.348; J=0.017; a=0.0327; b=3.112; k=1.400;
> Pag=Tn*ws; Pmag=P1n-Pag-(3*R1*I1n^2+0.005*Pn);
> B=3*R1*I0n^2; C=((R1+R2)*Pag^2)/(3*E1n^2); D=Pmag;
>
> a0=1;
> a1=2*b/(1-b);
> a2=(b/(b-1))^2+(2*B*(b-a)-C*m^2*(b-1)^2)/(D*alfa^k*(b-1)^2);
> a3=(2*C*m^2*b)/(D*alfa^k*(b-1));
> a4=-(C*m^2*b)/(D*alfa^k*(b-1)^2);
>
> Pol=[a0 0 a1 0 a2 0 a3 0 a4];
> P= roots(Pol)
>
> Here the positive real root is 0.2968. But I don't know how to do this by code ( not to see all roots and to say - This is real positive root). Please help me.

  Life might be easier for you if you take advantage of the fact that in this case you are solving a quartic (fourth order) polynomial equation in the square of the unknown. In such a case the eight roots are paired roots with opposite signs. Just select any real root and then take their positive square roots.

 r = roots([a0 a1 a2 a3 a4];
 P = sqrt(r(isreal(r)));

  Actually you might want to make that a little more robust by requiring only that the imaginary part of any root be smaller than some appropriate tolerance value. In some cases 'roots' might return a root that really should be entirely real as one having a very small imaginary part due to round-off error.

Roger Stafford

"Kiril " <kkirqkov@gmail.com> wrote in message <gs7uda$fkc$1@fred.mathworks.com>...
> I have vector with complex and real roots - P. I have to write code to find only positive real root(in this calculation its only one). This is my code:
> R1=1.230; X1=1.500; Xm=53.700; R2=0.787; X2=2.490;
> Vn=220; Tn=36.340;T=10; Pn=5.500; I1n=11.500; I0n=19.5; f=30;
> fn=50; P=4; sn=0.036; ns=1500; P1n=6451.5;E1n=209; m=T/Tn; alfa=f/fn;
> ws=157; nn=1446; wn=151.348; J=0.017; a=0.0327; b=3.112; k=1.400;
> Pag=Tn*ws; Pmag=P1n-Pag-(3*R1*I1n^2+0.005*Pn);
> B=3*R1*I0n^2; C=((R1+R2)*Pag^2)/(3*E1n^2); D=Pmag;
>
> a0=1;
> a1=2*b/(1-b);
> a2=(b/(b-1))^2+(2*B*(b-a)-C*m^2*(b-1)^2)/(D*alfa^k*(b-1)^2);
> a3=(2*C*m^2*b)/(D*alfa^k*(b-1));
> a4=-(C*m^2*b)/(D*alfa^k*(b-1)^2);
>
> Pol=[a0 0 a1 0 a2 0 a3 0 a4];
> P= roots(Pol)
>
> Here the positive real root is 0.2968. But I don't know how to do this by code ( not to see all roots and to say - This is real positive root). Please help me.

% If I understand, P is real & complex and you just want the real & > 0 part.
% So, add this to the end of your code (this works for this case),
% BUT THERE COULD BE NUMERICAL ACCURACY PROBLEMS IN GENERAL.

R = P(imag(P)==0);
PosR = R(R>0)

> % If I understand, P is real & complex and you just want the real & > 0 part.
> % So, add this to the end of your code (this works for this case),
> % BUT THERE COULD BE NUMERICAL ACCURACY PROBLEMS IN GENERAL.

Yes, but is in this calculation i always have only one real positive root i don't think that where is a problem with ACCURACY or if there is, please could u describe me why?

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <gs7vqa$q3i$1@fred.mathworks.com>...
> "Kiril " <kkirqkov@gmail.com> wrote in message <gs7uda$fkc$1@fred.mathworks.com>...
> > I have vector with complex and real roots - P. I have to write code to find only positive real root(in this calculation its only one).

> Just select any real root and then take their positive square roots.

  Kiril, my thinking wasn't straight in my previous response. You could very well have a negative real root to the quartic polynomial and the square root of that would be imaginary. (Blush!)

  Just do this:

 r = roots([a0 a1 a2 a3 a4]);
 P = sqrt(r(isreal(r)&real(r)>=0));

Roger Stafford

"Kiril " <kkirqkov@gmail.com> wrote in message <gs80et$daa$1@fred.mathworks.com>...
>
> > % If I understand, P is real & complex and you just want the real & > 0 part.
> > % So, add this to the end of your code (this works for this case),
> > % BUT THERE COULD BE NUMERICAL ACCURACY PROBLEMS IN GENERAL.
>
> Yes, but is in this calculation i always have only one real positive root i don't think that where is a problem with ACCURACY or if there is, please could u describe me why?

In general. roots might return a value of:

1000000 + 0.000001 i

Is that "real" enough for you?

If so, P(imag(P)==0 will fail.

That is all I was trying to covey.
In general, doing an equality test for floating point number is a bad idea.
They should be tested to within a tolerance (that you define).

Yes now i get it...(i hope so) :-)

"Kiril " <kkirqkov@gmail.com> wrote in message <gs89bu$3j4$1@fred.mathworks.com>...
> Yes now i get it...(i hope so) :-)
I did this
 p=[10 0 0 0 0 -1]
>> r=roots(p)
r =

  -0.5105 + 0.3709i
  -0.5105 - 0.3709i
   0.1950 + 0.6001i
   0.1950 - 0.6001i
   0.6310
for i=1:size(r,1);a=imag(r(i));if a==0;sol=r(i);end;end
>> sol
sol =
  0.6310