4
$\begingroup$

I am working with a complex set of two differential equations with many parameters. A simplified minimal working example would be the following system:

$$x'[t]=x^2[t]+y[t]+b$$ $$y'[t]=x[t]+y^2[t]-a$$

,where $a\in [0,1]$ and $b\in [0,1]$ are two parameters.

Now, depending on parameters the system may have 0,1 or 2 fixed points. What I want to do is Regionplot/Contourplot the parameterspace $[0,1]$x$[0,1]$ and shade regions/plot the boundary of regions for which there exist a different number of equilibria.

What I started with is plotting the nullclines as well as computing the fixed points of the system. I could just count the number of fixed points but I dont know how much that will help.

a = 0.5;
b = 0.5;
eq1=x^2 + y + b;
eq2=x + y^2 - a; 
ContourPlot[{eq1 == 0, eq2 == 0}, {x, -3, 3}, {y, -3, 
  3}, ContourStyle -> {Red, Blue}, ImageSize -> Scaled[0.3], 
 PlotLegends -> 
  LineLegend[{"\!\(\*SubscriptBox[\(\[PartialD]\), \(t\)]\)x=0", 
"\!\(\*SubscriptBox[\(\[PartialD]\), \(t\)]\)y=0"}], 
 FrameLabel -> {x, y}]

a = NSolve[{eq1 == 0 && eq2 == 0 }, {x, y}, Reals]
Length[a]

Ideally, I would like to end up with a plot similar to the following, where either regions are shaded or the boundary-line is plotted:

enter image description here

Edit: I want to emphasize that I am looking for numerical methods for a more complex problem.

Here is my current code, which works but still has 2 major issues:

1.I was only able to ListPlot3D the list and I would like to have a 2D projection of the data.

2.Empty sets, i.e. parameter configurations for which there exist no fixed point are still counted. I would like to get rid of those entries, while still preserving the value 0 in the plot.

eq1 = x^2 + y + b;
eq2 = x + y^2 - a;

c = Table[{a, b, x} /. 
NSolve[{eq1 == 0 && eq2 == 0 }, {x, y}, Reals], {a, 0, 1, 
0.005}, {b, 0, 1, 0.005}];

d = Flatten[Map[Length, c, {2}]];
e = Table[{a, b}, {a, 0, 1, 0.005}, {b, 0, 1, 0.005}];
f = Flatten[e, 1];
g = Flatten[Riffle[f, d]];
par = Partition[g, 3];

ListPointPlot3D[par, ColorFunction -> "Rainbow", Boxed -> False,  
 Axes -> {True, True, False}, ViewPoint -> {0, 0, 250}, 
 PlotLegends -> SwatchLegend[{Red, Blue}, {0, 2}]]

enter image description here

$\endgroup$
1
  • 2
    $\begingroup$ Your a has two meanings, which is usually a bad idea. $\endgroup$
    – corey979
    Jan 24, 2017 at 16:10

2 Answers 2

8
$\begingroup$

Update

Inspired by Chris K and showing larger parameter range (full range of real critical points: 0,1,2,3,4:

fun[a_, b_] := 
  Length[{x, y} /. 
     NSolve[{f[a, b, x, y] == {0, 0}}, {x, y}, Reals] /. {x, y} :> {}];
g[a_, b_] := (27 + 256 a^3 + 288 a b - 256 a^2 b^2 - 256 b^3);
h[a_, b_] := If[g[a, b] < 0, fun[a, b], 2];
lcp = ContourPlot[h[a, b], {a, -8, 8}, {b, -8, 8}, 
   ContourShading -> {Red, Green, Blue}, 
   PlotLegends -> SwatchLegend[{Red, Green, Blue}, {0, 2, 4}]];
DynamicModule[{p = {1/2, 1/2}}, 
 Row[{LocatorPane[Dynamic[p], 
    Show[lcp, Graphics[{PointSize[0.03], Yellow, Dynamic[Point[p]]}], 
     ImageSize -> 300]], 
   Dynamic@ContourPlot[
     Evaluate@Thread[(f[##, x, y] & @@ p) == 0], {x, -4, 4}, {y, -4, 
      4}, ContourStyle -> {Red, Blue}, ImageSize -> 300]}]]

enter image description here

Original Answer

For these quadratic relations it is relatively straight forward to use the discriminant, e.g.:

f[a_, b_, x_, y_] := {x^2 + y + b, x + y^2 - a}
DynamicModule[{p = {1/2, 1/2}},
 Row[{LocatorPane[Dynamic[p],
    ContourPlot[
     27 + 256 a^3 + 288 a b - 256 a^2 b^2 - 256 b^3, {a, 0, 1}, {b, 0,
       1}, Contours -> {{0}}, ContourShading -> {Red, Blue}, 
     ContourStyle -> Directive[Yellow, Thickness[0.02]], 
     Epilog -> {PointSize[0.03], Yellow, Dynamic[Point[p]]}, 
     ImageSize -> 300]],
   Dynamic@
    ContourPlot[
     Evaluate@Thread[(f[##, x, y] & @@ p) == 0], {x, -2, 2}, {y, -2, 
      2}, ContourStyle -> {Red, Blue}, ImageSize -> 300]
   }]
 ]

An example of no real zeroes:

enter image description here

"Approximately" 1: enter image description here

Two real zeros:

enter image description here

$\endgroup$
2
  • 1
    $\begingroup$ Thank you very much. This solves the problem for the minimal working example in the description. However, the problem I am dealing with contains polynomials of 5th order and hence I was not able to derive an analytical solution. It was kind of my bad in posting this overly simplified working example above. Could anyone try to get this working numerically? $\endgroup$
    – user45841
    Jan 25, 2017 at 8:50
  • $\begingroup$ @ParetoWilli yes, I'd also like to see a purely numerical approach! $\endgroup$
    – Chris K
    Jan 25, 2017 at 13:23
7
$\begingroup$

For this particular example you can find the bifurcations analytically by also setting the slopes of the isoclines to be equal.

eq1 := x^2 + y + b;
eq2 := x + y^2 - a;
bif = Solve[{eq1 == 0, eq2 == 0, D[eq1, x]/D[eq1, y] == D[eq2, x]/D[eq2, y]}, {x, y, a, b}]

(* {{y -> 1/(4 x), a -> (1 + 16 x^3)/(16 x^2), b -> (-1 - 4 x^3)/(4 x)}} *)

ParametricPlot[{a, b}/. bif, {x, -10, 10}, PlotRange -> {{-4, 8}, {-8, 4}}, Frame -> True, Axes -> False]

Mathematica graphics

I couldn't resist looking at a larger range of a and b, because there is another bifurcation between 2 and 4 equilibria. I don't have a way to automatically label the number of equilibria, but in the upper left there are 0, in the middle there are 2, and in the lower right there are 4.

Taking a look at the phase planes at three bifurcation points where a=2:

bif2 = NSolve[{eq1 == 0, eq2 == 0, D[eq1, x]/D[eq1, y] == D[eq2, x]/D[eq2, y], a == 2}, {x, y, a, b}]

(* {{a -> 2., x -> 1.98412, y -> 0.126, b -> -4.06275}, {a -> 2., x -> -0.169722, y -> -1.473, b -> 1.44419}, {a -> 2., x -> 0.185598, y -> 1.347, b -> -1.38144}} *)

GraphicsRow[Table[{a, b} = {a, b} /. bif2[[i]]; ContourPlot[{eq1 == 0, eq2 == 0}, {x, -3, 3}, {y, -3, 3},
ContourStyle -> {Red, Blue}], {i, 3}]]

Mathematica graphics

$\endgroup$
2
  • $\begingroup$ best answer. +1 an ashamed of mine though if I extend my parameter range my answer agrees with yours though my 'discriminant' does not allow different shading for 3 and 4 root case... $\endgroup$
    – ubpdqn
    Jan 25, 2017 at 4:16
  • 1
    $\begingroup$ @ubpdqn There is still plenty of room for improvement to really answer the question: 1) counting # of equilibria, 2) shading regions, 3) numerical methods for more complex models. Hopefully someone takes a stab at it. $\endgroup$
    – Chris K
    Jan 25, 2017 at 4:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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