bifurcation for Van der Pol osscillator

8 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
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
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.

Community Treasure Hunt

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

Start Hunting!