Bifurcation diagram using numerical solutions (e.g., ODE45)????

How plot the phases of a dynamical system?

I am trying to generate a bifurcation diagram for a predator prey interaction but I am struggling to find a way to plot it. This is the problem: Suppose the solution for the differential equations that describes the dynamic of the predator and the prey after a fixed number of iteration steps (to avoid transient) is unique, the dynamics are stable. If you plot the last few values of the iteration procedure will get a single point, which is good so far.

But lets say that the dynamics oscillate, a 2 phase limit cycle, for each component (predator or prey) they fluctuate around an equilibrium point with a constant max and min. If you plot the last steps of the iteration you will get many points close to each other in between the max and min of the behavior, not showing that the dynamic is a phase 2 limit cycle. The solution would be plot the max and min of a list of the last steps of the iteration. This does work, but only if the dynamic is stable of for a phase 2 cycle. If the system becomes a phase 4 or even if it has a chaotic behavior it will not be shown in my bifurcation diagram if I just plot max and min. The max and min plot will give you the amplitude of the dynamic but not if it is pahse 2, 4 , 8 or chaotic like in a classical logistic map.

So how to make a plot to show just the phases of the cycles or all the points like in a chaotic behavior?

To illustrate an example ploting max and min I will use the code provided by Dr. Young Kuang os his website. But how do I make a plot showing the phases of the cycles?

the first mfile (predprey.m) should contain the fuction describing the ODEs

 function [Ydot] = predprey(t,Y) 
r=2; K=1.4; s=1; c=0.5; a=2; d=0.1; 
Ydot(1) = r*Y(1)*(1-Y(1)/K)-s*Y(1)*Y(2)/(a*Y(1)+1); 
Ydot(2) = Y(2)*(c*s*Y(1)/(a*Y(1)+1)-d); 
Ydot=Ydot';

The second mfile will do the solution and the bifurcation diagram

 rect = [200 80 700 650]; %fix the window size and position
set(0, 'defaultfigureposition',rect);
global time IC r K s c a d
r=2; s=1; c=0.5; a=2; d=0.1;
option=odeset('AbsTol',1e-11,'RelTol',1e-11);
inc=[0.01:1:220]';
time=[ inc ] ;
limit=[ 200:1:220];
for K=0.1:.05 : 1.8,
IC=[0.5 .2];
[t,U]=ode15s('predprey',time,IC);
u1=U(limit,1);
u2=U(limit,2);
xmin=min(u1);
xmax=max(u1);
ymin=min(u2);
ymax=max(u2);
subplot(211),
plot(K, xmin,'*','MarkerSize',2);
plot(K, xmax,'+','MarkerSize',2);
%axis(ax);
hold on
xlim([.1 1.8]); ylim([0 1.6]);
end
title('Bifurcation diagram for the Holling type II predator-prey model','FontSize',12)
ylabel('x','Rotation',0);

[David Gonzalez]

Rafael,

In order to create a bifurcation diagram for a continuous system you have to create a Poincaré section for each variation of the parameter and then plot the result (removing transients). It is not about plotting the output variable vs parameter variation (in the 1D Logistic Map the output is the same Poincaré map).

Answer 1

[Steven Lord]

There are a couple problems with your code.

1. I see that you call ODESET, but you don't then pass that options structure into ODE15S. The ODESET function doesn't change the "global options for all ODE solver calls" or anything like that; if you want to use the changed options, you need to pass them into your call.
2. Your loop variable is K, and you have declared K as global in your script, but you have NOT declared K as global in predprey. In fact, since you haven't declared any of the variables global in predprey, that function will use the hard-coded values. I suspect you instead intended to have the ODE solver solve your ODEs for a different value of K each loop iteration. If that is correct, take a look at the link describing how to pass additional parameters into the ODE solvers on the documentation page for ODE15S.
3. Related to point 2, I recommend passing a function handle to predprey (@predprey) into the solver rather than the name of the function. It will greatly simplify the process of passing in additional inputs (which may eliminate the need for GLOBAL.)
4. Not directly related to your question, but you may eventually want to specify the OutputFcn option as @odephas2 to plot the 2-D phase plane for your model. See for example the orbitode example.

Answer 2

[PUSPENDU ROY]

How to plot birfurcation diagram of van der pol oscillators in scilab/matlab?