Stability of the equilibrium point (Matlab)

[PhMichael]

After analyzing a 3-DOF system, I've obtained the following 6th order characteristic polynomial:

P=\lambda^6+(6\beta-.8896)\lambda^5+(8\beta^2-3.5584\beta+6)\lambda^4+(\beta^3-.8896\beta^2+16\beta-3.5584)\lambda^3+(3\beta^2-1.7792\beta+8)\lambda^2+(3\beta-.8896)\lambda+1

Stability is determined by the real part of the roots of the this polynomial. Now, I need to vary numerically (matlab) the conditions for Re(\lambda), in other words, I want to plot the real part of the \lambda's as a function of \beta, where 0<\beta<1 ... How can I do that in matlab? - I tried to use the "solve" command, however, it doesn't even produce an explicit answer.

Thanks!

 

[viscousflow]

Ok thought about this and yes, what I posted does give the solution. Good thing \beta is bounded or else things would get complicated

 n=1000;
beta = linspace(0,1,n);

%coefficients of polynomial
coeff = [1 (6*beta-.8896) (8*beta.^2-3.5584*beta+6) ...
(beta.^3-.9986*beta.^2 +16*beta-3.5584) ...
(3*beta.^2-1.7792*beta+8) (3*beta*.8896) 1];

%solve for roots.
r=cell(n,1);
r_part = cell(n,1)

for i=1:n

r{i} = roots(coeff(i));
r_part{i} = real(r{i});%collect real part of roots

end

Now you have roots as a function of \beta

Also if you have it in symbolic format you do know how to do the same thing as above correct?
Replace the fourth line with:

coeff = wrev(coeffs(P))

where P is the symbolic equation in post #1, and of course substitute the value of beta for the linspace values.

 

[PhMichael]

Thanks a lot, viscousflow! ... After an intense struggle, I've managed to write a working code =) ... It goes like this:

G=0.8896;

beta=0:0.01:1;
A=zeros(6,size(beta));

a1 = 1 ;   
a2 = 6.*beta - G ;
a3 = 8.*beta.^2 - 4.*beta.*G + 6 ;
a4 = beta.^3 - G.*beta.^2 + 16.*beta - 4.*G ;
a5 = 3.*beta.^2 - 2.*beta.*G + 8 ;
a6 = 3.*beta - G ;
a7 = 1 ;

for i=1:length(beta)
    P=[a1 a2(i) a3(i) a4(i) a5(i) a6(i) a7];
    R=roots(P);
    Re=real(R);
    A(:,i)=Re'; 
end

plot(beta,A,'linewidth',3,'markersize',10)
grid;
ylabel ('Re(\lambda)')
xlabel ('\beta')

 

[viscousflow]

Thats good, however, I don't see how you'd get meaningful information from plotting Re(roots) versus Beta like that. You have 6 roots so I thought you'd want to see the trend of each root with beta. That's the advantage of the cell structure I had above, but that's just me.

Have a good day.

 

[PhMichael]

http://img831.imageshack.us/img831/1715/plotb.png

Uploaded with ImageShack.us

I need to plot Re(\lambda) as a function of \beta in order to find the stability range. According to the figure, the system is stable for \beta > 0.4457, while for \beta<0.4457, there exists at least one eigenvalue (root) with a non-negative real part and the system is unstable ... so far so good ...
Now, I want to determine \beta_H and the frequency \Omega_H, which correspond to the Hopf bifurcation ... How can I do that?