Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

When we do a StreamPlot, I want to show the bifurcation when $a = 0$ transitions to $a >0$, but do not see a better way to do this than the following.

    Animate[StreamPlot[{1,-x^4 + 5 a x^2-4 a^2}, {y,-5,5}, {x,-4,4},
                        PlotLabel -> Row[{"a = ", a}]], {a,-1,1}]

Is there a better way to show this bifurcation?

Maybe an animated gif like this Bifurcation diagrams for multiple equation systems) (this is a system and my example is not a system)?

Update

To be more clear, the original system is

$$x'(y) = -(x(y))^4 + 5a(x(y))^2 - 4a^2$$

share|improve this question

2 Answers 2

up vote 3 down vote accepted

Since you have just two variables, $x$ and $a$, you can make a 2D plot:

f[x_, a_] := -x^4 + 5 a x^2 - 4 a^2;

Show[StreamPlot[{f[x, a], 0}, {x, -4, 4}, {a, -1, 1}, StreamPoints -> Fine], 
 ContourPlot[f[x, a] == 0, {x, -4, 4}, {a, -1, 1}, ContourStyle -> ColorData[1, 2]]]

enter image description here

My explanatory comment below was incorrect, so I'm including a corrected version here. The parameter $a$ is on the vertical axis. For $a=1$ we have four equilibrium points. These are visible as intersections of the red curve with the slice $a=1$ parallel to the horizontal axis. Note that if you have $x'=f_a(x)=−x^4+5ax^2−4a^2$ with $a=1$, then for $|x|>2$ you have $x'<0$, for $2<|x|<1$ you have $x'>0$, and for $|x|<1$ you have $x'<0$, all of which are visible in the diagram in terms of the direction of the arrows.

Nevertheless, the version in Mark's answer is much more beautiful.

share|improve this answer

I agree with Rahul that, since you've got a single, one dimensional equation indexed by a single parameter, a bifurcation diagram is a natural way to visualize this situation. If you want an animation or dynamic image, you might highlight one particular phase line as a function of $a$. You could also display the slope field of the system right along side of the phase diagram. The result looks something like so.

enter image description here

Note that I prefer to draw the solution through a simple slope field, rather than through the groovy StreamPlot because I think that's a simpler concept - particularly, for undergraduate students who I interact with a lot. I've also oriented the vertical $x$ axis in both plots in the same direction and used graphics primitives there to grab a little more control over the image.

Here's code for an interactive version of this that also includes a locator on the slope field allowing you to specify the initial condition.

Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];
Clear[slopeField, slopeFieldWithSol];
slopeField[a_?(# < 0 &)] := VectorPlot[
   {1, -x^4 + 5 a x^2 - 4 a^2}, {t, -4, 4}, {x, -4, 4},
   VectorScale -> {0.03, Automatic, None},
   VectorStyle -> {Gray, Arrowheads[0]}];
slopeField[a_?(# >= 0 &)] := Show[{
    VectorPlot[
     {1, -x^4 + 5 a x^2 - 4 a^2}, {t, -4, 4}, {x, -4, 4},
     VectorScale -> {0.03, Automatic, None},
     VectorStyle -> {Gray, Arrowheads[0]}],
    Plot[{-2 Sqrt[a], -Sqrt[a], Sqrt[a], 2 Sqrt[a]}, {t, -4, 4},
     PlotStyle -> Directive[Black, Dashed]]
    }];
slopeFieldWithSol[a_, p_] := 
  Module[{t, x, t0, x0, eq, ic, eqs, sol, xf, ifd, tRange, plot},
   {t0, x0} = p;
   eq = x'[t] == -x[t]^4 + 5 a*x[t]^2 - 4 a^2;
   ic = x[t0] == x0;
   eqs = {eq, ic};
   Quiet[sol = First[NDSolve[eqs, x, {t, -4, 4}]], NDSolve::ndsz];
   xf = x /. sol;
   ifd = InterpolatingFunctionDomain[xf];
   tRange = {t, ifd[[1, 1]], ifd[[1, 2]]};
   plot = Plot[xf[t], Evaluate[tRange], 
     PlotStyle -> {Thick, Black}, PlotRange -> {{-4, 4}, {-4, 4}}];
   Show[{slopeField[a], plot}, PlotRange -> {{-4, 4}, {-4, 4}}]
   ];     
Clear[phaseLine, phaseDiagram, phaseDiagramWithLine];
phaseLine[a_?(# < 0 &), ___] := {ColorData[1, 2], 
   Arrow[{{a, 4}, {a, -4}}]};
phaseLine[0, eps_] = phaseLine[0.0, eps_] = {ColorData[1, 2],
    Arrow[{{0, 4}, {0, eps}}], Arrow[{{0, -eps}, {0, -4}}]
    };
phaseLine[a_?(# > 0 &), eps_] := {
   Arrowheads[Medium],
   ColorData[1, 2], Arrow[{{a, 4}, {a, 2 Sqrt[a] + eps}}],
   ColorData[1, 1], 
   Arrow[{{a, Sqrt[a] + eps}, {a, 2 Sqrt[a] - eps}}],
   ColorData[1, 2], Arrow[{{a, Sqrt[a] - eps}, {a, -Sqrt[a] + eps}}],
   ColorData[1, 1], 
   Arrow[{{a, -2 Sqrt[a] + eps}, {a, -Sqrt[a] - eps}}],
   ColorData[1, 2], Arrow[{{a, -2 Sqrt[a] - eps}, {a, -4}}]
   };
phaseDiagram = ParametricPlot[{{x^2/4, x}, {x^2, x}}, {x, -4, 4},
   PlotRange -> {{-1, 4}, {-4, 4}},
   PlotStyle -> Directive[Thickness[0.007], Black],
   Epilog -> {Opacity[0.4], 
     Table[phaseLine[a, 0.1], {a, -0.7, 4, 0.2}]}];
phaseDiagramWithLine[a_] := Show[{phaseDiagram,
    Graphics[{Thick, phaseLine[a, 0.1], If[a >= 0, {PointSize[Large],
        Point[{a, #} & /@ {2 Sqrt[a], 
           Sqrt[a], -Sqrt[a], -2 Sqrt[a]}]}, {}]}]}];
Manipulate[Grid[{{
    Show[slopeFieldWithSol[a, p], ImageSize -> 300],
    Show[phaseDiagramWithLine[a], 
     ImageSize -> {Automatic, 300}]
    }}, Spacings -> 3],
 {{a, 0}, -1, 4}, {{p, {0, 1}}, {-4, -4}, {4, 4}, Locator}]

The animated gif was generated in a manner almost identical to the Manipulate. I just specified a particular p value $(-1.2,1.2)$, changed Manipulate to Table (which required me to change the a specification slightly), and Exported the result to a GIF.

pics = Table[Grid[{{
  Show[slopeFieldWithSol[a, {-1.2, 1.2}], ImageSize -> 300],
  Show[phaseDiagramWithLine[a], 
   ImageSize -> {Automatic, 300}]
  }}, Spacings -> 3],
 {a, -1, 4, 0.1}];
Export["temp.gif", pics]
share|improve this answer
1  
Absolutely gorgeous! Thank you, I need to study this in greater detail! +1 May I ask what I think is a simple question? How do you get the animated GIF? –  Amzoti Mar 1 at 13:41
    
@Amzoti Glad you like it! I've added just a little commentary that might explain a few things. Have fun! –  Mark McClure Mar 1 at 13:44
    
@Amzoti I also just noticed that you asked about the GIF, so I added that as well. –  Mark McClure Mar 1 at 13:49
    
Greatly appreciated, I have so much to learn regarding Mathematica! Regards –  Amzoti Mar 1 at 13:50

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.