bifurcation for Van der Pol osscillator

5 views (last 30 days)
Sreeshankar P
Sreeshankar P on 26 Nov 2021
Commented: Sreeshankar P on 6 Dec 2021
  1 Comment
Sreeshankar P
Sreeshankar P on 26 Nov 2021
Edited: Walter Roberson on 27 Nov 2021
%%Initialize the environment
close all;
clear all;
clc;
%%Model parameters
r = 0.806;
a = 15;
b = 16;
c = 17;
e = 0.333;
d = 0.3;
h = 0.01;
K = 200;
%%Time parameters
dt = 0.01;
N = 10000;
%%Set-up figure and axes
figure;
ax(1) = subplot(2,1,1);
hold on
xlabel ('m');
ylabel ('x');
ax(2) = subplot(2,1,2);
hold on
xlabel ('m');
ylabel ('y');
%%Main loop
for m=6:1:22
x = zeros(N,1);
y = zeros(N,1);
t = zeros(N,1);
x(1) = 0.7;
y(1) = 0.11;
t(1) = 0;
for i=1:N
t(i+1) = t(i) + dt;
x(i+1) = x(i) + dt*( r*x(i)*(1-x(i)/K)-m*x(i)*y(i)/(a*x (i)+b*y(i)+c) );
y(i+1) = y(i) + dt*( e*m*x(i)*y(i)/(a*x(i)+b*y(i)+c)-d*y(i)-h*y(i)^2 );
end
plot(ax(1),m,x,'color','blue','marker','.');
plot(ax(2),m,y,'color','blue','marker','.');
end
I am not getting the bifurcation diagram
please help

Sign in to comment.

Answers (2)

Walter Roberson
Walter Roberson on 28 Nov 2021
I am not sure what you are expecting ?
%%Initialize the environment
%%Model parameters
r = 0.806;
a = 15;
b = 16;
c = 17;
e = 0.333;
d = 0.3;
h = 0.01;
K = 200;
%%Time parameters
dt = 0.01;
N = 1000; %10000;
%%Set-up figure and axes
figure;
ax(1) = subplot(3,1,1);
hold(ax(1),'on');
xlabel(ax(1),'m');
ylabel(ax(1),'x');
ax(2) = subplot(3,1,2);
hold(ax(2),'on')
xlabel(ax(2),'m');
ylabel(ax(2),'y');
ax(3) = subplot(3,1,3);
hold(ax(3),'on')
xlabel(ax(3),'x');
ylabel(ax(3),'y');
cmap = parula(22);
%%Main loop
for m=6:1:22
x = zeros(N,1);
y = zeros(N,1);
t = zeros(N,1);
x(1) = 0.7;
y(1) = 0.11;
t(1) = 0;
for i=1:N
t(i+1) = t(i) + dt;
x(i+1) = x(i) + dt*( r*x(i)*(1-x(i)/K)-m*x(i)*y(i)/(a*x (i)+b*y(i)+c) );
y(i+1) = y(i) + dt*( e*m*x(i)*y(i)/(a*x(i)+b*y(i)+c)-d*y(i)-h*y(i)^2 );
end
plot(ax(1),m,x,'color','blue','marker','.');
plot(ax(2),m,y,'color','blue','marker','.');
plot(ax(3),x,y,'color',cmap(m,:),'Marker','.', 'DisplayName', string(m));
end
legend(ax(3), 'show')
  7 Comments
Walter Roberson
Walter Roberson on 2 Dec 2021
When global is needed at all (not clear that it is necessary here) then it must be followed by a space separated list of pure identifiers with no indexing. Equations are not permitted in a "global" statement.
Sreeshankar P
Sreeshankar P on 2 Dec 2021
I dont get it sir..
please help in sloving the problem..
I have no idea what to do now..

Sign in to comment.


Sreeshankar P
Sreeshankar P on 6 Dec 2021
function res = pport();
clear;
gamma=0.17;
alpha = 1;
f = 15;
omega = (2: 0.001: 4.1);
x = [0 1];
ts= 25;
tspan=[0,ts];
% opt=odeset('RelTol',le-6,'AbsTol',le-9);
figure;
hold on;
for omega =(3.9 : 0.001 :4.1 )
Y0=[v0,x0];
[~,Y]=ode45(@(t,y)VDP(f,omega,alpha,t,gamma),tspan,[v0]);
v=Y(:,1);
x=Y(:,2);
plot(x,v,'k')
drawnow
end
hold off;
% FUNCTION
function res= VDP(y,f,omega,alpha,t,gamma)
v=y(1);
x=y(2);
dv= f*cos(omega.*t)-gamma*x*(1-(x^2))-alpha.*x;
t;
res[dv,v];
Can some one please help, i am not getting the bifurcation plot for the codes...
  6 Comments
Walter Roberson
Walter Roberson on 6 Dec 2021
The earlier code used 1e-6 but your latest code uses le-6 . 1e-6 begins with the digit for the number One, but le-6 begins with lower-case L
Sreeshankar P
Sreeshankar P on 6 Dec 2021
could you please correct the error , i have been trying to to do this for a long time..

Sign in to comment.

Categories

Find more on Modeling in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!