bifurcation for Van der Pol osscillator

[Sreeshankar P]

%%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

Answer 1

[Walter Roberson]

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')

[Sreeshankar P]

Sorry,

This is not the result i am looking for.

I Would like to get the codes for bifurcation diagram and FFT on the equation given using the given parameters..

Any kind of help is appreciated.

[Walter Roberson]

http://www.math.le.ac.uk/people/rld8/ma1251/lab3.html

[Sreeshankar P]

global alpha=1 f=15 gamma=0.17 omega=[3.9,4.1];
x0=[1,0];
t=zeros(x0)
for omega = 3.9:0.001:4.1
    x(1) = 0.5; 
    for n = 1:1000
        x(1) = f*cos(omega*t)-gamma*x(2)*(1-(x(1))^2)-alpha*x(1); 
    end 
    for n = 1:1000
        x(n+1) =  fcos(omega*t)-gamma*y(1)*(1-(x(n))^2)-alpha*x(n); 
    end 
    plot(r*ones(Nplot,1), x, '.', 'markersize', 2); 
    hold on; 
end 
title('Bifurcation diagram of the logistic map'); 
xlabel('r'); ylabel('x_n'); 
set(gca, 'xlim', [2.5 4.0]); 
hold off;

I alredy tried that but i am not getting the plot for bifurcation and fft..

can you please help

[Walter Roberson]

global alpha=1 f=15 gamma=0.17 omega=[3.9,4.1];

That line is not valid MATLAB syntax.

It is also not at all clear why you would want those variables to be global variables.

Note that your

for omega = 3.9:0.001:4.1

line would overwrite omega each iteration of the loop.

[Sreeshankar P]

Sir,

I am soo new to this matlab and have been trying to solve this equation for more than a week..

Could you please help me ...

[Walter Roberson]

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]

I dont get it sir..

please help in sloving the problem..

I have no idea what to do now..

Answer 2

[Sreeshankar P]

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

[Walter Roberson]

Your function VDP expects y as its first input, but your anonymous function is not passing y into VDP.

[Walter Roberson]

 res[dv,v];

That is not valid syntax. You are missing an assignment operation. Also you must return a column vector, not a row vector.

[Sreeshankar P]

so that should be res [ dv , y]??

[Sreeshankar P]

pport

Error using <=

Not enough input arguments.

Error in pport (line 11)

opt=odeset('RelTol',le-6,AbsTol,le-9);

this is the error showing up

[Walter Roberson]

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]

could you please correct the error , i have been trying to to do this for a long time..