Off[General::"spell"];
Off[General::"spell1"];

Hopf Bifurcation Theorem

This notebook is designed to carry out the calculations required to distinguish between supercritical and subcritical Hopf bifurcations.  All you need to do in enter the model information as requested and evaluate the expressions. The example given is that of the text for the model of sexual selection:

Gz (a p - c z (1 + (z^2))) - 
      B b p (1 + (p^2))
B (a p - c z (1 + (z^2))) - 
      Gp b p (1 + (p^2))

In the text we use the composite parameter, γ = -b Gp - c Gz + a B, as the bifurcation parameter. Therefore, we must rewrite the above equations in terms of this parameter. To match the notation of Recipe 11.2, however, let's call this parameter μ instead of γ. You can rewrite the model in terms of this bifurcation parameter in a number of different ways, and here we will simply solve the equation μ = -b Gp - c Gz + a B for B, and then substitute this into the equations of the model. We get B = ((b Gp + c Gz + μ)/a), and therefore, in terms of the bifurcation parameter, μ, the model can be written as

Gz (a p - 
            c z (1 + 
                  (z^2))) - ((b Gp + c Gz + μ)/a) b p (1 + 
            (p^2))
((b Gp + c Gz + μ)/a) (a p - c z (1 + (z^2))) - 
      Gp b p (1 + (p^2))

In the series of 3 steps that follow, you are required to input model specific information where the text is BLUE. You are required to evaluate these, and all other expressions, in the order that they appear.

Step 1: First, you need to make sure that the equilibirum of interest in the model is (0,0).  If the equilibrium of interest is not the origin, then the model should be re-written in terms of two new variables whose equilibrium values ARE (0,0). This is done by taking the original model, say dx/dt = f(x,y) and dy/dy = g(x,y), and the equilibirum of interest in this model,  x=x* and y=y*, and then defining the two new variables, v=x-x* and w=y-y*. Then you derive differential equations for these new variables as dv/dt = dx/dt = f(x,y) = f(v+x*,w+y*) and dw/dt = dy/dt = g(x,y) = g(v+x*,w+y*). This gives a version of the original model, but rewritten so that the critical value of the bifurcation parameter is μstar, and the equilibrium of interst is (0,0). For the sexual selection model of the text, we can skip this step since the equilibirum of interest in the original variables is already (0,0). 

Step 2: Input Model Details: Define two function, f and g, as the right-hand side of the differential equations in the model. For the present case we have:

f[x_, y_, μ_] := 
        Gz (a y - 
                c x (1 + 
                      (x^2))) - ((b Gp + c Gz + μ)/a) b y (1 + 
                (y^2))
g[x_, 
          y_, μ_] := ((b Gp + c Gz + μ)/a) (a y - 
                c x (1 + (x^2))) - Gp b y (1 + (y^2))

Next, input the critical value of the bifurcation parameter where stability changes (in this case it is μ*=0):

μstar = 0;

The following then serves as a check that you have done everything above correctly. f and g should equal zero at the origin (if they don't, you need to go back and make sure that the equilbirum of interest in your model is the origin):

Simplify[{f[x, y, μ], g[x, y, μ]} /. {x -> 0, 
          y -> 0}]

The Jacobian matrix for the model is then:

jac[x_, y_, μ_] := 
      Outer[D, {f[x, y, μ], g[x, y, μ]}, {x, y}]

At the equil point of interest (i.e., the origin) the Jacobian simplifies to:

jaceq[μ_] := 
      jac[x, y, μ] /. {x -> 0, y -> 0}

and its eigenvlaues are given by

Eigenvalues[jaceq[μ]]

Now if you entered everthing correctly, you should see that these eigenvalues are purely imgainary at the critical value of μ that you entered above. As a check, we can calculate them at this critical value as

Simplify[
      Eigenvalues[jaceq[μ] /. μ -> μstar]]

So far, this does the basic calculations for the Hopf Bifurcation Theorem that area readily calculated by hand. The above eigenvalues at μ=μ* should be purely imaginary. 

To check if the bifurcation is supercritical or subcritical, and thus if the limit cycle is stable or unstable (or if there is a collection of neutral cycles around te equilibirum) we need to does some more work. This is where this notebook comes in handy since the calculations are quite involved and tedious.

Step 3: To proceed further with the calculations, we typically first have to make another change of coordinates so that the Jacobian matrix is in a simple form. Specifically, we want the Jacobian evaluated at the origina with μ=μ* to have the form(({{0, 1}, {-K, 0}})).

This is readily automated by simply evaluating the following commands:

This first command calculates the appropriate change in coordinates:

p = {{aa, bb}, {1, 1}} /. Solve[
            jaceq[0] . {{aa, bb}, {1, 1}} == {{aa, bb}, {1, 1}} . {{0, 
                    1}, {-dden, 0}}, {aa, bb, 
              dden}][[1]]
den = 
      dden /. Solve[
            jaceq[0] . {{aa, bb}, {1, 1}} == {{aa, bb}, {1, 1}} . {{0, 
                    1}, {-dden, 0}}, {aa, bb, 
              dden}][[1]]

The model, written in these new coordinates is then given by

fnew[xnew_, 
        ynew_, μ] := Simplify[(Inverse[p] . {f[x, y, μ], 
                  g[x, y, μ]}) /. {x -> (p . {xnew, 
                        ynew})[[1]], 
              y -> (p . {xnew, 
                        ynew})[[2]]}][[1]]
gnew[xnew_, 
        ynew_, μ] := Simplify[(Inverse[p] . {f[x, y, μ], 
                  g[x, y, μ]}) /. {x -> (p . {xnew, 
                        ynew})[[1]], 
              y -> (p . {xnew, 
                        ynew})[[2]]}][[2]]

The Jacobian in these new coordinates is given by

jacnew[xnew_, ynew_, μ_] := 
      Outer[D, {fnew[xnew, ynew, μ], gnew[xnew, ynew, μ]}, {xnew, 
          ynew}]

As a check, let's see if this Jacobian, at the origin and μ=μstar, is in the desired form:

FullSimplify[
      jacnew[xnew, ynew, μ] /. {xnew -> 0, 
          ynew -> 0, μ -> μstar}]

Simlarly, we can check that the model in these new coordinates, has eigenvlaues that become purely imaginary when μ=μstar

FullSimplify[
      Eigenvalues[
        jacnew[xnew, ynew, μ] /. {xnew -> 0, 
            ynew -> 0}]]

Then, the final required calculation is obtained by evaluating the following (Marsden and McCracken 1976 The Hopf Bifurcation and Its Applications.  Springer-Verlag, NY):

v[xnew_, 
        ynew_, μ_] := ((3 Pi)/(4  den)) (D[
                fnew[xnew, ynew, μ], xnew, xnew, xnew] + 
              D[fnew[xnew, ynew, μ], xnew, ynew, ynew] + 
              D[gnew[xnew, ynew, μ], xnew, xnew, ynew] + 
              D[gnew[xnew, ynew, μ], ynew, ynew, 
                ynew]) + ((3 Pi)/(4 (den^2))) (D[
                  fnew[xnew, ynew, μ], xnew, 
                  ynew] (D[fnew[xnew, ynew, μ], xnew, xnew] + 
                    D[fnew[xnew, ynew, μ], ynew, ynew]) + 
              D[gnew[xnew, ynew, μ], xnew, 
                  ynew] (D[gnew[xnew, ynew, μ], xnew, xnew] + 
                    D[gnew[xnew, ynew, μ], ynew, ynew]) + 
              D[fnew[xnew, ynew, μ], xnew, xnew] D[
                  gnew[xnew, ynew, μ], xnew, xnew] - 
              D[fnew[xnew, ynew, μ], ynew, ynew] D[
                  gnew[xnew, ynew, μ], ynew, ynew])

And if the next expression is negative, then we have a supercritical bifurcation. If it is positive we have a subcritical bifurcation.

FullSimplify[
      v[xnew, ynew, μ] /. {xnew -> 0, 
          ynew -> 0, μ -> 0}]

It might require some work to determine the sign of the resulting expression.