Say I have a matrix that depends on a parameter $u$. The following may be an example:

a = {{1, u^3 + 2, 1}, {0, 1, u}, {1/u^2, 0, 1}};

My real problem is 12x12 matrix with difficult parametric dependencies (it won't make much sense to paste it here). I would like to find the values of $u$ such that there exist imaginary eigenvalues and no eigenvalues with positive real parts. How may I achieve this with Mathematica?

  • $\begingroup$ In at least your $3\times3$ case, it's easy: evaluate the discriminant of your characteristic polynomial, and determine the relevant conditions for it to be negative. $12\times12$ would be symbolically difficult, I think. $\endgroup$ – J. M.'s discontentment Jan 7 '17 at 14:55
  • $\begingroup$ is the real problem just one parameter? $\endgroup$ – george2079 Jan 7 '17 at 14:56
  • $\begingroup$ @george2079 yes one parameter only. $\endgroup$ – Mirko Aveta Jan 7 '17 at 15:00
  • $\begingroup$ @J.M. Indeed! In the real problem it couldn't evaluate the symbolic expression of the eigenvalues in the first place. $\endgroup$ – Mirko Aveta Jan 7 '17 at 15:02
  • 2
    $\begingroup$ at 12x12 I'm sure you need to go about it numerically. construct a function that returns the max of the real part of the eigenvalues for a numeric u. Simply plotting that will get you started. $\endgroup$ – george2079 Jan 7 '17 at 15:37

Based on @george2079's suggestion, let's start by plotting the maximum real eigenvalue:

a = {{1, u^3 + 2, 1}, {0, 1, u}, {1/u^2, 0, 1}};
maxev[v_?NumericQ] := Max[Re[Eigenvalues[a /. u -> v]]]
Plot[maxev[u], {u, -10, 10}, PlotRange -> {0, 10}]

Mathematica graphics

Your 3x3 example matrix doesn't seem to have any points with Re[eigenvalue]=0, so let me modify it:

a = {{-1, u^3 + 2, 1}, {0, -1, u}, {1/u^2, 0, -1}};
Plot[maxev[u], {u, -10, 10}]

Mathematica graphics

This seems to have two roots, which we can find with FindRoot:

bif1 = FindRoot[maxev[u], {u, -1}]
bif2 = FindRoot[maxev[u], {u, -2}]

(* {u -> -0.621916} *)
(* {u -> -1.38647} *)

Now it remains to see if either of these have imaginary parts:

Chop[Eigenvalues[a /. bif1]]
Chop[Eigenvalues[a /. bif2]]

(* {-3., 0. + 0.643852 I, 0. - 0.643852 I} *)
(* {-1.5 + 0.479364 I, -1.5 - 0.479364 I, 0} *)

Thus, in this example, bif1={u -> -0.621916} seems to be what you're looking for (a Hopf bifurcation point?)

| improve this answer | |
  • $\begingroup$ Yes I am looking for a bifurcation point. The difficulty lays in the fact that my system has a large dimension. It would be nice to have a general code that globally searches such points. This is what I was really looking for. May I ask a silly question? Why did you ask your code to search for tha Max real eigenvalue? Is it not the eigenvalue with smallest Real part the one that gets closest in crossing the imaginary axis? $\endgroup$ – Mirko Aveta Jan 15 '17 at 19:04
  • $\begingroup$ @MirkoAveta FindRoot will work well if you have a good initial guess. Plotting can help. I'm sure there are other examples of code to find all roots on this site too. As for why Max[Re[eigenvalues]], when Max[Re[eigenvalues]]=0 you know that all the other eigenvalues have real parts that are <=0. $\endgroup$ – Chris K Jan 15 '17 at 19:17
  • $\begingroup$ thanks for your patient reply. I'll take a look around for global root solutions. Your answer was what I was looking for! $\endgroup$ – Mirko Aveta Jan 15 '17 at 19:21

This is what I was looking for. Precision is given through eps. For my needs eps=0.1 is sufficient.

EigQ[v_] := 
 AllTrue[Re[Eigenvalues[matrixA /. u -> v]], # < 0 &]
ustar[kbeta_] := 
 With[{}, u = 1; While[EigQ[u] == True, u = u + eps]; u]

This codes works until one knows that the bigger is u the closer is the wanted condition. A general solution for this problem I believe is much more complicated to search.

| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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