# Stability of the equilibrium point (Matlab)

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!

Related Engineering and Comp Sci Homework Help News on Phys.org
Ok thought about this and yes, what I posted does give the solution. Good thing $$\beta$$ is bounded or else things would get complicated

Code:
 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:
Code:
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.

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

Code:
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')

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.

http://img831.imageshack.us/img831/1715/plotb.png [Broken]

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?

Last edited by a moderator: