Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

So I have this dynamical system given by: $$ \left\{\begin{aligned} x' &= a(y-\phi(x))\\ y' &= x-y+z\\ z' &= -by \end{aligned}\right. $$ where $\phi(x) = \mu x^3 - \nu x$ and $a,b,\mu,\nu$ are positive real parameters. Obviously, the origin is an equilibrium point. I now want to determine what kind of equilibrium this is, which depends on the signs of the real parts of the eigenvalues of the matrix of the linearized system. The linearized system is given by: $$ X'=AX\\ A=\begin{pmatrix} a\nu & a & 0\\ 1 & -1 & 1\\ 0 & -b & 0 \end{pmatrix} $$ The characteristic polynomial is then given by: $$ p(\lambda)_{a,b,s} = -\lambda^3 + (a\nu-1)\lambda^2 + (a\nu-b+a)\lambda + ab\nu $$ which is pure hell to solve (which may have something to do with the fact that this family of systems is chaotic...). Anyway, I tried throwing this into Mathematica and it just spits out about 4 pages of symbols. Then I tried Reduce in order to find out whether these eigenvalues had negative or positive real parts. Sadly, this gives me about 15 pages of symbols, which doesn't really help. Does anybody know a good way to deal with this?

Tl;Dr

How can I compute the signs of the real parts of above matrix A, for all possible values of $a,b,\nu>0$? Criteria for them to be positive will also do.

Edit

I tried the following input:

A[a_, b_, \[Sigma]_] := {{a*\[Sigma], a, 0}, {1, -1, 1}, {0, -b, 0}}
eig[a_, b_, \[Sigma]_] := Eigensystem[A[a, b, \[Sigma]]]
Table[Reduce[eig[a, b, \[Sigma]][[1]][[i]] > 0, {a, b, \[Sigma]}, Reals], {i, 1, 3}]

The output is displayed by Mathematica as:

A very large output was generated. Here is a sample of it:
{<<1>>, <<1>>, <<1>>}

Basically it is a list with a lot of different cases. This is of course to be expected, but the real problem is that there are tons of Roots in it...

share|improve this question
    
It would be helpful if you posted the code you tried and some part of the output. –  David Carraher Apr 1 '13 at 22:22
    
A visualization suggests that except along the boundaries, there are always two negative and one positive real parts: RegionPlot3D[And @@ MapThread[Equal, {{-1, -1, 1}, Sign[Re[#]] & /@ (l /. Solve[-l^3 + (a n - 1) l^2 + (a n - b + a) l + a b n == 0, {l}])}], {a, 0, 5}, {b, 0, 5}, {n, 0, 5}]. –  whuber Apr 1 '13 at 22:30
    
An analysis with the trace and determinant confirms that there is one negative and two positive real parts when $0<a\nu \leq 1$. But I can't show it for the other case... –  Eric Spreen Apr 1 '13 at 22:42

2 Answers 2

up vote 2 down vote accepted

The strategy pursued in a comment is to examine the signs of the coefficients of the characteristic polynomial for information about the signs of its roots. Mathematica can help.

There are two cases to consider: all roots are real or two of them are complex conjugates of each other and the third is real.

In the first case, let the roots be $x$, $y$, and $z$. From the positivity of the determinant $a b \nu$ it is immediate that an even number of the roots are negative. Can exactly two be negative? We can find out by matching coefficients

m = {{a n, a, 0}, {1, -1, 1}, {0, -b, 0}};
match = And @@ (Equal @@@ Transpose[{CoefficientList[Expand[(l - x) (l - y) (l - z)], l], 
      CoefficientList[-CharacteristicPolynomial[m, l], l]}])

$-x y z==-a b n\ \&\&\ x y+x z+y z==-a+b-a n\ \&\&\ -x-y-z==1-a n$

and looking for an instance that satisfies all the inequalities simultaneously:

FindInstance[match && a > 0 && b > 0 && n > 0 && x > 0 && y > 0 && z > 0, {a, b, n, x, y, z}]

$\{\}$

There are none. In the second case, let the roots be $x\pm i y$ and $z$. This time, matching coefficients gives

match = And @@ (Equal @@@ Transpose[
  {CoefficientList[Expand[(l - (x+I y)) (l - (x-I y)) (l - z)], l], 
   CoefficientList[-CharacteristicPolynomial[m, l], l]}])

$-x^2 z-y^2 z==-a b n\ \&\&\ x^2+y^2+2 x z==-a+b-a n\ \&\&\ -2 x-z==1-a n$

Noticing that $y$ enters only as $y^2$, we may retain the inequality $y\ge 0$ (because the more conditions there are, the easier it may be to prove that no instance exists). Once again we search for an instance with all positive real parts; that is, $x \gt 0$:

FindInstance[match && a > 0 && b > 0 && n > 0 && x > 0 && y >= 0 && z > 0, {a, b, n, x, y, z}]

$\{\}$

There are none. (This method also works by leaving out y >= 0 but it takes about a minute to execute, rather than a couple of seconds.)

In either case we conclude that exactly two of the real parts are negative and one of them is positive.

share|improve this answer
1  
Thank you so much! This really does the trick. I never knew of the existence of FindInstance. I guess I should really plough through the documentation when I get the chance once. ;) –  Eric Spreen Apr 2 '13 at 10:03

I'm kind of shocked no one mentioned this, its perhaps one of the best kept secrets amongst people who study dynamical systems (surprisingly it rarly gets mentioned in an intro DE class) but by the Routh–Hurwitz stability criterion this equilibrium is unstable.

$$ p(\lambda)_{a,b,s} = -\lambda^3 + (a\nu-1)\lambda^2 + (a\nu-b+a)\lambda + ab\nu, $$

In this case, $a_3=-ab\nu<0$, which means the characteristic equation has at least one root with a positive real part. This may help you save a lot of time in the future, although whuber's solution actually gets at some of the reasoning why the Routh–Hurwitz stability criterion is true.

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.