Chapter 8, Section 8.8.4, Figure 8.15

Define the function h to be the Hamiltonian for the mixed resonance system (see Eq.(8.186)).

h[x_, y_, δ_, c_]:=
   (1/2)*δ*(x^2+y^2)+(1/4)*(x^2+y^2)^2+
   2*x*Sqrt[x^2+y^2+c]

Define the functions f and g as given in Eq.(8.185).

f[δ_,X_]:=δ+X^2;
g[i_, X_, c_]:=
(-1)^(i+1)*2*(Sqrt[X^2+c]/X + X/Sqrt[X^2+c]);

Calculate the derivatives of f and g.

df =  D[f[δ, X], X]

dg  =  D[g[0, X, c], X]

Plot both derivatives (for c = –5) as a function of X = R.  Note that the singularity causes problems if start with lower values of X.

Plot[{df, (dg  /.  c -> -5)}, {X, 2.3, 10}]

Find the point of intersection.  At this point the tangents to f and g have equal gradient.

Xb  =  
      X  /.  FindRoot[
          df == (dg  /.  c -> -5),  {X, 2.5}]

Find the corresponding values of g and δb.

gval  =  g[0, Xb, -5]

δbval  =  gval  -  Xb^2

Derive th plot f for \(TraditionalForm\`δ\&_\) = –15, –10, δbval, 0, 5 and 10, and g for c = –5, i = 0, 1.  This shows the tangential touching for \(TraditionalForm\`δ\&_\) =δbval, the bifurcation value (see Fig.8.15a).  Note the problems that Plot has when dealing with some values of X.

\(p = 

      Plot[{f[\(-15\), X], f[\(-10\), X], f[δbval, \ X], \ f[0, X], \ 
          f[5, X], f[10, X], g[0, X, \(-5\)], g[1, X, \(-5\)]}, {X, 0, 
          6}, \ 
PlotRange -> {{0, 6}, {\(-15\), 15}}, 
        AxesOrigin -> {0, \(-15\)}, PlotPoints -> 100, AspectRatio -> 1, \ 
        Frame -> True, \ DisplayFunction -> Identity]\)

Create a graphics object to denote the singularity line at X = (Sqrt[5]).

pline  =  
        Graphics[Line[{{Sqrt[5], -15}, {Sqrt[5], 15}}]];

Display the plot and the singularity line (see Fig.8.15a).

Show[p, pline, DisplayFunction -> $DisplayFunction]

Define a function, plotit, that uses Eq.(8.187) to plot the relationship between \(TraditionalForm\`δ\&_\) andthe value of x associated with the equilibrium points for the specified value c.  The function also plots a horizontal line for the \(TraditionalForm\`δ\&_\) associated with the bifurcation value and the singularity lines at x = ±(Sqrt[5]).

\(plotit[c_, \ dbv_] := 

      Module[{}, 

        p1 = Plot[\((\(-x^3\) - 2*Sqrt[x^2 + c] - 2*x^2/Sqrt[x^2 + c])\)/
              x, {x, \(-6\), 6}, \ AspectRatio -> 1, \ Frame -> True, 
            PlotRange -> {{\(-6\), 6}, {\(-30\), 10}}, 
            DisplayFunction -> Identity]; 

        p2\  = \ Graphics[Line[{{\(-6\), dbv}, {6, dbv}}]]; 

        p3\  = Graphics[
            Line[{{\(-Sqrt[\(-c\)]\), \(-30\)}, {\(-Sqrt[\(-c\)]\), 10}}]]; 

        p4\  = \ Graphics[
            Line[{{Sqrt[\(-c\)], \(-30\)}, {Sqrt[\(-c\)], 10}}]]; 

        Show[p1, p2, \ p3, p4, 
          DisplayFunction -> $DisplayFunction]]\)

Plot the function at c = –5 (see Fig.8.15b).  Note the problems encountered by plotit for certain values.

plotit[-5, δbval]

Calculate the value of the Hamiltonian associated with the bifurcation point.

hval0 = h[Xb, 0, δbval,  -5]

Plot the curves of constant Hamiltonian for h = –55, –50, –40, –31.72847, hval0, –10, 10 and 30 (see Fig.8.15c).  Again note the problems that ContourPlot encounters because of the singularities.

cp = 
      ContourPlot[
        h[x, y, δbval,  -5], {x, -7, 5}, {y, -6, 6},  
        ContourShading -> False, AxesOrigin -> {0, 0}, 
        Contours -> {-55, -50, -40, -31.72847, 
            hval0, -10, 10, 30}, 
AspectRatio -> 1, 
        PlotPoints -> 150]

The forbidden region is a circle of radius (Sqrt[5]) centred on the origin.  Create a graphics object for this circle ...

circ = Graphics[Circle[{0, 0}, Sqrt[5]]];

... and plot it together with the contour plot.

Show[cp, circ, PlotRange -> {{-7, 5}, {-6, 6}}, 
      AspectRatio -> 1]

Plot h as a function of x for y = 0 and c = –5.

Plot[h[x, 0, δbval, -5], {x, -7, 5}]