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.

I am interested in constructing a bifurcation diagram for some of my parameters (especially for β and γ) in the dynamical system given in the code below. I want to see how parameter changes affect the stability of the system, but I am having some trouble in estimating the critical points in Mathematica and then using them to get the critical point orbits. I found this post helpful,

Mathematica code for Bifurcation Diagram

but the post refers to a logistic-equation-like case. I do not how to apply it to my multiple equation case. Any help or pointers would be most appreciated

The following code incldes some wrong code intended to get the critical points, so my apologies in advance.

f = {L, ψ, d} /. 
  NDSolve[{d'[
        t] == (1/z) (α - β (d[t]/Y[t]) - γ L[t] - δ ψ[
             t] - τ Y[t]) + 0.75 (d[t]/Y[t]) d[t], 
      L'[t] == (j (Y[t]/K) - e)/N, ψ'[t] == p L[t] - 0.35, 
      Y[t] == (I + E + α - β (d[t]/Y[t]) - γ L[t] - δ ψ[
            t])/(1 - (1 - τ) (ψ[t] + (1 - s) (1 - ψ[t]) - m z)), 
      L[0] == 0.50, ψ[0] == 0.50, 
      d[0] == 0} /. {α -> 0.05, γ -> 0.75, δ -> 0.75, 
      X -> 1, β -> 0.05, j -> 0.45, τ -> 0.35, I -> 1, s -> 0.85, p -> 0.75, 
      E -> 1, m -> 1.2, e -> 0.035, z -> 5.5, K -> 3.5, N -> 4},
      {L, ψ, d, Y}, {t, 0, 300}, MaxSteps -> 1000000][[1]]

l = {L[t], ψ[t], d[t]};    
cps[L] = f /. Quiet[Solve[f[[L'[t]]] == 0], All];

Show[GraphicsArray[
  Table[Plot[f[[i]][t], {t, 0, 300}, PlotRange -> All, PlotStyle -> Blue, Filling -> 0, 
    AxesLabel -> TraditionalForm /@ {t, l[[i]]}, DisplayFunction -> Identity], {i, 3}]]]

data2 = Table[Evaluate[f[[#]][t] & /@ {1, 2, 3}], {t, 0, 300}];
Export["data2.csv", data2]

ListLinePlot[Transpose@data2, Filling -> 0, AxesLabel -> {t, {L, ψ, d}}, PlotRange -> All]

Graphics3D[Point[data2], BoxRatios -> 1, AxesLabel -> TraditionalForm /@ l]
ParametricPlot3D[Evaluate[Append[#[t] & /@ f, Red]], {t, 0, 300}, BoxRatios -> {1, 1, 1}, 
 PlotRange -> All, PlotPoints -> 1500, AxesLabel -> TraditionalForm /@ l]
share|improve this question
    
Vyperultra: Thanks Mr. Abbasi! Sure, this equation is a Rossler-inspired system with three state equations \dot{d}, \dot{L} and \dot{\psi} which determine an output variable (which is Y) which in turn feedbacks into the equations. The rest are parameters. As the code shows, it gives dampened three dimensional cyclical dynamics. –  vyperultra Dec 6 '12 at 18:46
add comment

1 Answer

First, I have a few general comments. The simplest bifurcation diagrams for differential equations involve a single parameter in a single equation and there are several illustrations of these on the Wolfram Demonstrations site. More generally, of course, a bifurcation occurs when we see a qualitative change in the behavior of a system as some parameter changes. For a system of equations, we could use the eigensystem of the linearization of the system in the neighborhood of an equilibrium to identify bifuractions. You, however, appear to have four differential-algebraic equations with 16 parameters. Also, a number of these (K, I, E, N) are actually reserved symbols; I'd avoid that.

I don't think I'm going to wade into this system that I know nothing about but I can present some ideas to study this kind of thing in the context of an example that I understand, namely the Selkov model presented on the Demonstrations site: http://demonstrations.wolfram.com/HopfBifurcationInTheSelkovModel/

That demonstration shows you how to study the bifurcation using NDSolve. It might make sense to study it using the groovy new ParametricNDSolve. First, the system is the following:

$$x' = -x + a y + x^2 y$$ $$y' = b - a y - x^2 y$$

We can solve this as follows.

pf = ParametricNDSolveValue[{
  x'[t] == -x[t] + a*y[t] + x[t]^2 y[t],
  y'[t] == b - a*y[t] - x[t]^2*y[t],
  x[0] == 0, y[0] == 2},
 {x, y}, {t, 0, 100}, {a, b}];

We can investigate the behavior with respect to the parameters as follows. If you hold $a=0.1$ and let be range, a two Hopf bifurcations (changes from fixed point to cycle and back) should be evident.

Manipulate[
 Block[{$PerformanceGoal = "Quality"},
  Show[{
    ContourPlot[-x + a*y + x^2 y == 0, {x, 0, 3}, {y, 0, 3},
     ContourStyle -> Dashed],
    ContourPlot[b - a*y - x^2*y == 0, {x, 0, 3}, {y, 0, 3},
     ContourStyle -> Dashed],
    StreamPlot[{-x + a*y + x^2 y, b - a*y - x^2*y},
     {x, 0, 3}, {y, 0, 3},
     StreamStyle -> Directive[Opacity[0.5]]],
    ParametricPlot[
     Evaluate[Through[pf[a, b][t]]],
     {t, 0, 100}, PlotRange -> {{0, 3}, {0, 3}},
     PlotStyle -> Directive[Thickness[0.007], Black]]}]],
 {{a, 0.1}, 0, 1}, {b, 0, 1}]

Portrayed as an animation, this looks like so:

enter image description here

share|improve this answer
add comment

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.