Stability of the equilibrium point (Matlab)

  • Thread starter PhMichael
  • Start date
  • #1
134
0
After analyzing a 3-DOF system, I've obtained the following 6th order characteristic polynomial:

[tex]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[/tex]

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

Thanks!
 

Answers and Replies

  • #2
272
0
Ok thought about this and yes, what I posted does give the solution. Good thing [tex]\beta [/tex] 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 [tex]\beta [/tex]

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:
  • #3
134
0
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')
 
  • #4
272
0
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.
 
  • #5
134
0
http://img831.imageshack.us/img831/1715/plotb.png [Broken]

Uploaded with ImageShack.us

I need to plot [tex]Re(\lambda)[/tex] as a function of [tex]\beta[/tex] in order to find the stability range. According to the figure, the system is stable for [tex]\beta > 0.4457 [/tex], while for [tex]\beta<0.4457[/tex], 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 [tex]\beta_H[/tex] and the frequency [tex]\Omega_H[/tex], which correspond to the Hopf bifurcation ... How can I do that?
 
Last edited by a moderator:

Related Threads on Stability of the equilibrium point (Matlab)

  • Last Post
Replies
9
Views
975
Replies
5
Views
457
Replies
2
Views
334
  • Last Post
Replies
4
Views
2K
Replies
3
Views
1K
  • Last Post
Replies
2
Views
3K
Replies
2
Views
5K
Top