Using find command to find bifurcation points

9 visualizaciones (últimos 30 días)
Ross Mannion
Ross Mannion el 18 de Nov. de 2019
Comentada: Ross Mannion el 18 de Nov. de 2019
I have code which plots a graph showing the levels of Notch and Delta in a pair of cells. The graph is a bifurcation diagram with 2 bifurcation points. After plotting, I have attempted to find the 2 values of 'a' where the bifurcations happen. I have tried using the find command but only receive empty vectors when I run the code. Function file is at the end of the code.
% Specify initial conditions
dinit=[0 1];
ninit=[1 0];
initialconditions=[dinit; ninit]';
% Set values of a
a=logspace(-8,2,200);
for i=1:numel(a)
% Apply ode45
[t,y]=ode45(@(t,y)twocellfunct(t,y,a(i)),[0 200],initialconditions);
% Calculate maximum value of Notch in cells 1 and 2
mx=max([y(end,3), y(end,4)]);
% Calculate minimum value of Notch in cells 1 and 2
mn=min([y(end,3), y(end,4)]);
% Storing the max/min values of notch for each a
M(:,i)=[mx mn]';
end
% Plot
semilogx(a,M(1,:),'r');
hold on
semilogx(a,M(2,:),'b');
xlabel('a');
ylabel('notch level');
y=ylabel('notch level', 'rot', 90);
set(y, 'Units', 'Normalized', 'Position', [-0.07, 0.5, 0]);
a1=find(abs(M(1,:)-M(2,:))<eps,1,'first');
a2=find(abs(M(1,:)-M(2,:))<eps,1,'last');
function l = twocellfunct(t,y,a)
% Specifying parameters
b=100;
v=1;
k=2;
h=2;
% RHS functions
f=@(x)(x.^k./(a+x.^k));
g=@(x)(1./(1+b.*x.^h));
l = [v.*(g(y(3))-y(1)); v.*(g(y(4))-y(2)); f(y(2))-y(3); f(y(1))-y(4)];
end
  3 comentarios
the cyclist
the cyclist el 18 de Nov. de 2019
Without being able to run your code, I speculated that eps is too tight a tolerance. Have you tried making that larger?
Ross Mannion
Ross Mannion el 18 de Nov. de 2019
yes sorry, added it now, I have tried changing eps but seem to either get an empty vector or just 1 which also isnt right

Iniciar sesión para comentar.

Respuesta aceptada

the cyclist
the cyclist el 18 de Nov. de 2019
Editada: the cyclist el 18 de Nov. de 2019
Here is a loglog plot of the difference between M(1,:) and M(2,:).
I think what you are actually trying to do is find the first and last points where the difference is larger than some tolerance, not smaller. That's the bifurcation.
It's an inexact science, but the following tolerance will get close. You might be able to refine the estimate with additional rules.
tol = 5.e-3;
a1=find(abs(M(1,:)-M(2,:))>tol,1,'first');
a2=find(abs(M(1,:)-M(2,:))>tol,1,'last');
  3 comentarios
the cyclist
the cyclist el 18 de Nov. de 2019
a1 and a2 are the indices into a, not the values. You'll get what you want from
a(a1)
a(a2)
Ross Mannion
Ross Mannion el 18 de Nov. de 2019
Aah I see, thank you for your help!

Iniciar sesión para comentar.

Más respuestas (0)

Categorías

Más información sobre Startup and Shutdown en Help Center y File Exchange.

Etiquetas

Community Treasure Hunt

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

Start Hunting!

Translated by