All the Maple work for this project can be found on our (relatively messy) Maple Worksheet - right-click to save

Our first step in analyzing this model was to reduce the two-equation system into its state equation. This was done by solving the first equation for u.

    a1v3 + lambda
u = -------------
     Bv2 + lambda

We then substituted this equation into the second equation, resulting in a single cubic state equation. Using this state equation, we could determine the steady states of the equation by solving and plotting the bifurcation diagram using the implicitplot Maple function. As before, our variable parameter was gamma2.

As we varied gamma2 from .2 to .1, we noticed that the graph turned in on itself, indicating the presence of 2 or 3 solutions to the equation for certain values of lambda:

Graph for gamma2 = .20 Graph for gamma2 = .10

To find these solutions, we solved the equation with the initial conditions:

v(0) = 0

Dv(0) = 0 

By plotting a number of the stable solutions, we discovered that as gamma2 increased, the two stable steady states got closer and closer together. We continued to plot solutions until we discovered the point in which the solutions converge: gamma2 = .122. You can see that at around lambda = 4.04 there is only one solution; the graph shows a perfectly vertical drop.

In order to find any hopf bifurcation points, we needed to define the jacobian of the system:

Hopf bifurcation points are found when the trace of the jacobian is equal to zero and the determinant of the jacobian is greater than zero. These conditions guarantee (according to the quadratic formula) that the eigenvalues will be pure imaginary complex conjugates. Using maple, we chose a value for v in the jacobian to get an approximate value for lambda. Then using these two values, we could approximate an area where the hopf point might be. Then we solved the jacobian for values of v and lambda in that area, looking for values that would satisfy the hopf condtions.

We discovered that there were 2 hopf bifurcation points for values of gamma2 below .157355. at that point, the hopf points converge at 10 decimal places. At .157356, there are again 2 hopf points, but when we looked at .1573578, we saw that the two solutions converged to 10 decimal places. We began to try larger and larger values of gamma2. At around gamma2 = .1574 the hopf points seemed to dissappear completely.

However, Maple does have inconsistencies in the commands that we used to find the hopf points:

Notice that the ranges are very nearly the same, but three different hopf points appear. As v varies, Maple may make errors in its floating point math. It seems that as the difference in the position of those hopf points gets smaller, there are errors in calculating the exact position.

For the purposes of this analysis, we shall consider the existence of hopf points for gamma values greater than .157355 to be Maple errors. This value is consistent with the "official" value.

Below are some XRK plots of solutions for specific values of gamma2 and lambda:

gamma2 = .12, lambda = .1343740933

This plot is centered around a hopf point when gamma2 = .12. It is interesting to note that this periodic solution decays extremely slowly as it approaches the limit cycle. After running XRK simulations for about 1.5 million seconds according to the time on XRK, we finally found a stable limit cycle:

The ranges on this graph are:

.6345 <= x <= .635 ; .0176 <= y <= .072 

gamma2 = .1, lambda = 3.10771

This is a phase plot of two stable steady states for a gamma value of .1. If you look carefully, you will notice that there is a saddle point in between the two solutions. The saddle corresponds to the unstable solution for that particular value of lambda.

gamma2 = .157355, lambda = .1490012434

This plot is of the values of gamma2 and lambda in which the hopf points converge.

gamma2 = .18, lambda = 4

This shows a single spiral sink solution.

This shows the solution for the same gamma2, but lambda = 1. These two graphs differ slightly, but show the same behavior.

gamma2 = .122, lambda = 4.041203429

Exactly at the interface where the two stable steady states converge, we get solutions that look like this. Notice that this solution is a little less "spiral" than the dual steady state graph above.

Ethan Deneault (eand@wpi.edu)