Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
Positive real root

Subject: Positive real root

From: Kiril

Date: 16 Apr, 2009 18:46:02

Message: 1 of 9

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.

Subject: Positive real root

From: Roger Stafford

Date: 16 Apr, 2009 19:10:03

Message: 2 of 9

"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

Subject: Positive real root

From: someone

Date: 16 Apr, 2009 19:13:02

Message: 3 of 9

"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)

Subject: Positive real root

From: Kiril

Date: 16 Apr, 2009 19:21:02

Message: 4 of 9


> % 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?

Subject: Positive real root

From: Roger Stafford

Date: 16 Apr, 2009 21:02:02

Message: 5 of 9

"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

Subject: Positive real root

From: someone

Date: 16 Apr, 2009 21:26:02

Message: 6 of 9

"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).

Subject: Positive real root

From: Kiril

Date: 16 Apr, 2009 21:53:02

Message: 7 of 9

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

Subject: Positive real root

From: v

Date: 9 Jun, 2010 06:30:22

Message: 8 of 9

"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

Subject: Positive real root

From: DeeDeeDee

Date: 16 Nov, 2013 20:44:05

Message: 9 of 9

"v " <vijays.nitt@gmail.com> wrote in message <huncdu$3bk$1@fred.mathworks.com>...
> "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

Hi,
Can anyone explain how to find all of the real and positive solutions? For example, if

> r =
>
> -0.5105 + 0.3709i
> -0.5105 - 0.3709i
> 0.1950 + 0.6001i
> 0.1950 - 0.6001i
> 0.6310
> 0.5783

How can this be modified to output both of the real and positive solutions?
  

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us