9 - Unstable Branches

The purpose of this program is to plot branches of fixed points and/or periodic orbits when they are not sinks. In this situation we cannot use the program for the bifurcation diagram (Program 7) and, unless the branches can be determined analytically, we must use a different approach. Here is one possibility.

f[a_, x_] := -Exp[x] + a

It is not difficult to verify that for every a there is one and only one fixed point of the system.  In fact the equation f(a,x)=x is equivalent to 
-(e^x)=x–a. The line y=x–a intersects the graph of  -(e^x) in one and only one point. The solution is negative for a<1 and positive for a>1. Consequently, the fixed point is stable for a<1 and unstable for a>1. We can use this program to plot the branch of fixed points for values of a>1 and <1.

a = 0;

q1 = 
      Table[{a, x /. FindRoot[f[a, x] - x == 0, {x, 0}], a = a + 0.01}, {
          200}];

q2 = q1[[Range[1, 200], Range[1, 2]]];

q3 = ListPlot[q2, PlotStyle -> PointSize[0.01]];

For a>1 there is a periodic orbit of period 2 which is a sink. We can plot the two branches of periodic points of period 2.

a = 1;

q4 = 
      Table[{a, x /. FindRoot[f[a, f[a, x]] - x == 0, {x, 1.5}], 
          a = a + 0.005}, {100}];

q5 = q4[[Range[1, 100], Range[1, 2]]];

q6 = ListPlot[q5, PlotStyle -> PointSize[0.01]];

a = 1;

q7 = 
      Table[{a, x /. FindRoot[f[a, f[a, x]] - x == 0, {x, -1.25}], 
          a = a + 0.005}, {100}];

q8 = q7[[Range[1, 100], Range[1, 2]]];

q9 = ListPlot[q8, PlotStyle -> PointSize[0.01]];

Show[q3, q6, q9];

We can see that two branches of periodic points of period 2, which exist for a>1 . The periodic orbit is a sink for a very close to 1. Here is the bifurcation diagram obtained with Program 7.

f[x_] := -Exp[x] + r

OrbitDiagram[start_, finish_, n_,x0_, init_, final_] := Show[
	Graphics[{PointSize[0.001], 
		Table[
			Map[Point[{r, #}] &, NestList[f, Nest[f, x0, init], final]],
			{r, finish, start, (start - finish)/n}
		]}
	], 
	AxesOrigin->{start,0},Axes->True
];

p1 = OrbitDiagram[0,  1,  100, -.2,  400, 4];

p2 = OrbitDiagram[1,  1.4,  100, .5,  400,  4];

Show[p1, p2,  Axes -> True];