Message Boards Message Boards

0
|
5446 Views
|
7 Replies
|
2 Total Likes
View groups...
Share
Share this post:

Find the value of b for which one of the eigenvalues crosses the imaginary?

Dear members,

I have a Matrix that that depends on a parameter "b"

A={{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
  0}, {0, 0, 0, 0, 0, 0, -2 \[Pi], 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0,
   0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, -4 \[Pi], 0, 0, 0, 0, 0, 1, 0, 0, 
  0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, -6 \[Pi], 0, 0, 0, 0, 0,
   1, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, -8 \[Pi], 0, 
  0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 
  0, -10 \[Pi], 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0}, {0, 2 \[Pi], 0, 0, 
  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0}, {0, 0, 
  4 \[Pi], 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 
  0}, {0, 0, 0, 6 \[Pi], 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
  1, 0, 0}, {0, 0, 0, 0, 8 \[Pi], 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
  0, 0, 0, 1, 0}, {0, 0, 0, 0, 0, 10 \[Pi], 0, 0, 0, 0, 0, 0, 0, 0, 0,
   0, 0, 0, 0, 0, 0, 1}, {-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.1, 0, 
  0, 0, 0, 0, -(b/2), 0, 0, 0, 0}, {0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
  0, -0.1, 0, 0, 0, 0, -2 \[Pi], -(b/2), 0, 0, 0}, {0, 0, -1, 0, 0, 0,
   0, 0, 0, 0, 0, 0, 0, -0.1, 0, 0, 0, b/2, -4 \[Pi], -(b/2), 0, 
  0}, {0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.1, 0, 0, 0, b/
  2, -6 \[Pi], -(b/2), 0}, {0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
   0, -0.1, 0, 0, 0, b/2, -8 \[Pi], -(b/2)}, {0, 0, 0, 0, 0, -1, 0, 0,
   0, 0, 0, 0, 0, 0, 0, 0, -0.1, 0, 0, 0, b/2, -10 \[Pi]}, {0, 0, 0, 
  0, 0, 0, -1, 0, 0, 0, 0, -b, 2 \[Pi], b/2, 0, 0, 0, -0.1, 0, 0, 0, 
  0}, {0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, -(b/2), 4 \[Pi], b/2, 0, 
  0, 0, -0.1, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 
  0, -(b/2), 6 \[Pi], b/2, 0, 0, 0, -0.1, 0, 0}, {0, 0, 0, 0, 0, 0, 0,
   0, 0, -1, 0, 0, 0, 0, -(b/2), 8 \[Pi], b/2, 0, 0, 0, -0.1, 0}, {0, 
  0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, -(b/2), 10 \[Pi], 0, 0, 
  0, 0, -0.1}}

To find the value of b for which one of the eigenvalues crosses the imaginary I am doing this

FindRoot[Max[Re[Eigenvalues[A]]] == 0, {b, 1}]

varies errors rise when executing FindRoot. For intance, "Unable to find all roots of the characteristic polynomial" It seem to me that FindRoot does not take "b" as a variable. Am I using FindRoot correctly?

Thanks for any help you can provide.

Jesús

7 Replies

Daniel,

I am trying to determine the transition curves, i.e. the lines that separate unstable from stable regions in space parameters, so yes it is to do with things such as bifurcations. I am going to check your interesting answer!

Thanks!

Jesus

I guess you are looking for Hopf bifurcations of some system? Anyway, one can use some tricks with Resultant to get the right polynomial. For this I will recommend using an exact input since the underlying polynomial algebra might otherwise get into trouble.

mat = Rationalize[{{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0,
      0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, -2 \[Pi], 0, 0, 0, 0, 0, 1, 
     0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, -4 \[Pi], 0, 0,
      0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 
     0, -6 \[Pi], 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0,
      0, 0, 0, 0, 0, -8 \[Pi], 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 
     0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -10 \[Pi], 0, 0, 0, 0, 0, 1, 
     0, 0, 0, 0, 0}, {0, 2 \[Pi], 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
     0, 0, 0, 1, 0, 0, 0, 0}, {0, 0, 4 \[Pi], 0, 0, 0, 0, 0, 0, 0, 0, 
     0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0}, {0, 0, 0, 6 \[Pi], 0, 0, 0, 0, 
     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0}, {0, 0, 0, 0, 8 \[Pi], 
     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0}, {0, 0, 0, 0, 
     0, 10 \[Pi], 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
     1}, {-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.1, 0, 0, 0, 0, 
     0, -(b/2), 0, 0, 0, 0}, {0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
     0, -0.1, 0, 0, 0, 0, -2 \[Pi], -(b/2), 0, 0, 0}, {0, 0, -1, 0, 0,
      0, 0, 0, 0, 0, 0, 0, 0, -0.1, 0, 0, 0, b/2, -4 \[Pi], -(b/2), 0,
      0}, {0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.1, 0, 0, 0, 
     b/2, -6 \[Pi], -(b/2), 0}, {0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 
     0, 0, 0, -0.1, 0, 0, 0, b/2, -8 \[Pi], -(b/2)}, {0, 0, 0, 0, 
     0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.1, 0, 0, 0, 
     b/2, -10 \[Pi]}, {0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, -b, 2 \[Pi], 
     b/2, 0, 0, 0, -0.1, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 
     0, 0, -(b/2), 4 \[Pi], b/2, 0, 0, 0, -0.1, 0, 0, 0}, {0, 0, 0, 0,
      0, 0, 0, 0, -1, 0, 0, 0, 0, -(b/2), 6 \[Pi], b/2, 0, 0, 0, -0.1,
      0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, -(b/2), 
     8 \[Pi], b/2, 0, 0, 0, -0.1, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 
     0, -1, 0, 0, 0, 0, -(b/2), 10 \[Pi], 0, 0, 0, 0, -0.1}}];

We compute the Jacobian and then do symmetrizing and antisymmetrizing.

jac = CharacteristicPolynomial[mat, t];
jacOdd = Expand[(jac - (jac /. t -> -t))/t];
jacEven = jac + (jac /. t -> -t);

Take the resultant of those. This gives a polynomial such that real zeros correspond to eigenvalues with real parts of zero.

AbsoluteTiming[resjac = Resultant[jacOdd, jacEven, t];]
LeafCount[resjac]
Variables[resjac]

(* Out[117]= {61.904456, Null}

Out[118]= 39656

Out[119]= {b}

In[120]:= Exponent[resjac, b]

Out[120]= 312 *)

See if there is a nontrivial factorization.

fax = FactorList[resjac];
Length[fax]
fax[[All, 2]]

(* Out[194]= 2

Out[195]= {-1, 2} *)

The first factor is always numeric. The exponent of 2 indicates the resultant is a square. We extract zeros from the squarefree part (faster, and more accurate since it avoids multiplicity).

zeros = b /. NSolve[fax[[2, 1]] == 0, b];
realzeros = Select[zeros, FreeQ[#, Complex] &]

(* Out[145]= {-224.800771278, -55.3181307708, -32.2609903819, \
-11.6818472022, -11.6632745351, -11.3588528657, -10.5718118805, \
-6.85508765395, 6.85508765395, 10.5718118805, 11.3588528657, \
11.6632745351, 11.6818472022, 32.2609903819, 55.3181307708, \
224.800771278} *)

How good a job did this do? We check on the minimal in magnitude real eigenvalues for each of these values of b.

Map[Min[Abs[Re[Eigenvalues[mat /. b -> #]]]] &, realzeros]

(* Out[197]= {9.00261133707*10^-6, 0.0000726273146539, 0.000232560670565,
  7.7200827675*10^-10, 8.02005961642*10^-10, 4.31832347658*10^-11, 
 1.42108547152*10^-14, 4.4408920985*10^-14, 1.95399252334*10^-14, 
 1.06581410364*10^-14, 4.91602314412*10^-11, 8.42769298792*10^-10, 
 8.31033228343*10^-10, 0.000232560670582, 0.0000726273146503, 
 9.00261133652*10^-6} *)

None are dreadful and some are quite good.

We can cut down on the degree by observing that the Jacobian is a polynomial in b^2. We reduce to a polynomial in b, again extract a resultant, solve, and take square roots from the eventual result.

jac2 = jac /. b -> Sqrt[b];
jacOdd = Expand[(jac2 - (jac2 /. t -> -t))/t];
jacEven = jac2 + (jac2 /. t -> -t);
res = Resultant[jacOdd, jacEven, t];

Again look for a nontrivial factorization.

faxB = FactorList[res];
Length[faxB]
faxB[[All, 2]]

(* Out[215]= 2

Out[216]= {-1, 2} *)

zerosB = b /. NSolve[faxB[[2, 1]] == 0, b];
realzerosB = Sqrt[Select[zerosB, FreeQ[#, Complex] &]]

(* Out[222]= {6.85508765395, 10.5718118805, 11.3588528657, \
11.6632745351, 11.6818472022, 32.2609903819, 55.31813zerosB = b /. NSolve[faxB[[2, 1]] == 0, b]07708, \
224.800771278} *)

Again check the eigenvalues for each.

Map[Min[Abs[Re[Eigenvalues[mat /. b -> #]]]] &, realzeros]

(* Out[223]= {5.3290705182*10^-15, 1.49880108324*10^-14, 
 4.24180690572*10^-11, 8.24576129599*10^-10, 
 7.98321200812*10^-10, 0.000232560670574, 0.0000726273146503, 
 9.00261133729*10^-6} *)

No improvement to speak of, but there is a better chance now of doing higher precision root finding in finite time, should that be desired. I did that actually, and it showed that the last three values do not give eigenvalues where the real part vanishes. I'm not sure what goes wrong here, other than that this resultant vanishing might be a necessary but not sufficient condition. Offhand I would expect the "bad" values to comprise a factor of the resultant, over the rationals (since that's how resultants usually show spurious roots to systems). Apparently that isn't the case here.

POSTED BY: Daniel Lichtblau

Hi again

Where can I find more about the procedure you used in your answer?. Some of the values of b you have found are very accurate.

I could not figure out if the values in between are all stable or unstable.

Best regards

Jesus

One reference is a tech report from 2000 by Karin Gatermann: "Counting stable solutions of sparse polynomial systems in chemistry"

https://www.google.com/url?sa=t&rct=j&q=&esrc=s&source=web&cd=1&cad=rja&uact=8&ved=0ahUKEwjiwaOm5eDbAhWp44MKHQ5vAPUQFggtMAA&url=https%3A%2F%2Fopus4.kobv.de%2Fopus4-zib%2Ffrontdoor%2Fdeliver%2Findex%2FdocId%2F600%2Ffile%2FZR-00-32.ps&usg=AOvVaw3nKM2IUGHWqcbuhyii8BEL

See section 4, "Finding Hopf bifurcations with resultants"

I also did this computation in a different way, basically substituting re+I*im for the char poly variable, setting 're' to zero, and taking the resultant of the explicit real and imaginary parts with respect to im to get a polynomial in the parameter b. It gives the same result(ant).

As for assessing stability in between the critical values of b, it suffices to take one numeric value between any adjacent pair (as well as one positive value beyond either extreme). Evaluate the eigenvalues and see if all real parts are less than zero. In principle you would also need to add the ordinary critical points, that is, values of b for which an eigenvalue is zero. As it turns out (see below), there are no positive values of b for which this happen.

In[253]:= NSolve[jac /. t -> 0, b]

Out[253]= {{b -> -0.202419223967 - 
    38.001292257 I}, {b -> -0.202419223967 + 
    38.001292257 I}, {b -> -0.155569770065 - 
    14.9249163595 I}, {b -> -0.155569770065 + 14.9249163595 I}, {b -> 
   0.155569770065 - 14.9249163595 I}, {b -> 
   0.155569770065 + 14.9249163595 I}, {b -> 
   0.202419223967 - 38.001292257 I}, {b -> 
   0.202419223967 + 38.001292257 I}}
POSTED BY: Daniel Lichtblau

Yes, it works!!

Thanks

Jesus

Make the dependency explicit, please

In[31]:= Clear[A]
A[b_?NumericQ] := 
  N[{{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
     0}, {0, 0, 0, 0, 0, 0, -2 \[Pi], 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0,
      0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, -4 \[Pi], 0, 0, 0, 0, 0, 1, 
     0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, -6 \[Pi], 0, 0,
      0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 
     0, -8 \[Pi], 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0,
      0, 0, 0, 0, 0, -10 \[Pi], 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0}, {0, 
     2 \[Pi], 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0,
      0}, {0, 0, 4 \[Pi], 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
      1, 0, 0, 0}, {0, 0, 0, 6 \[Pi], 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
      0, 0, 0, 0, 1, 0, 0}, {0, 0, 0, 0, 8 \[Pi], 0, 0, 0, 0, 0, 0, 0,
      0, 0, 0, 0, 0, 0, 0, 0, 1, 0}, {0, 0, 0, 0, 0, 10 \[Pi], 0, 0, 
     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1}, {-1, 0, 0, 0, 0, 0, 0,
      0, 0, 0, 0, -0.1, 0, 0, 0, 0, 0, -(b/2), 0, 0, 0, 0}, {0, -1, 0,
      0, 0, 0, 0, 0, 0, 0, 0, 0, -0.1, 0, 0, 0, 0, -2 \[Pi], -(b/2), 
     0, 0, 0}, {0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.1, 0, 0, 0,
      b/2, -4 \[Pi], -(b/2), 0, 0}, {0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0,
      0, 0, 0, -0.1, 0, 0, 0, b/2, -6 \[Pi], -(b/2), 0}, {0, 0, 0, 
     0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.1, 0, 0, 0, 
     b/2, -8 \[Pi], -(b/2)}, {0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 
     0, 0, 0, -0.1, 0, 0, 0, b/2, -10 \[Pi]}, {0, 0, 0, 0, 0, 0, -1, 
     0, 0, 0, 0, -b, 2 \[Pi], b/2, 0, 0, 0, -0.1, 0, 0, 0, 0}, {0, 0, 
     0, 0, 0, 0, 0, -1, 0, 0, 0, 0, -(b/2), 4 \[Pi], b/2, 0, 0, 
     0, -0.1, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 
     0, -(b/2), 6 \[Pi], b/2, 0, 0, 0, -0.1, 0, 0}, {0, 0, 0, 0, 0, 0,
      0, 0, 0, -1, 0, 0, 0, 0, -(b/2), 8 \[Pi], b/2, 0, 0, 0, -0.1, 
     0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, -(b/2), 
     10 \[Pi], 0, 0, 0, 0, -0.1}}];

In[50]:= Clear[rico]
rico[x_?NumericQ] := Max[Re[Evaluate[Eigenvalues[A[y] //. y -> x]]]]

In[52]:= FindRoot[rico[q], {q, 1., .9}]
Out[52]= {q -> 6.85509}

and

Plot[rico[x], {x, -1, 9}

Shows it seems to work. With Eigenvalues[] one can have trouble in numerical functions because it can also generate symbolic results.

POSTED BY: Dent de Lion

The error msg comes from applying Eigenvalues to A.
Eigenvalues::eival: Unable to find all roots of the characteristic polynomial. It returns unevaluated.

POSTED BY: Frank Kampas
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract