Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
pairwise invasibility plot

Subject: pairwise invasibility plot

From: J G

Date: 8 May, 2012 03:33:06

Message: 1 of 2

Does anyone know how to do this in Matlab?

Subject: pairwise invasibility plot

From: jan

Date: 16 Jun, 2013 21:58:17

Message: 2 of 2

"J G" wrote in message <joa45i$c2i$1@newscl01ah.mathworks.com>...
> Does anyone know how to do this in Matlab?

Hi,

I solved the problem like this:

1. You need to think of a trade off function. it may be arbitrary as long as it serves your needs, for example, the uptake-rate should be decreasing with increasing yield.

2. You get a steady state function of the resource depending on your trait. Set a resource-zero point at which point the species is able to live. this is a constant and should cross the steady state function at least twice.

3. look for trait values at which the resource steady state is given. plot these values.

I wrote this code, it is straightforward from 2. . warning: it is not commented.

%% Initialize

% % ex1
% x = linspace(0,10);
% S_0 = 2.8;

% ex2
x = linspace(2,5);
S_0 = 7;

ns = zeros(2,1);
inex = zeros(1,2);
step = 0.1;
counter = 0;
runner = 3;
options = optimset('TolX', 1e-13);
R_0 = 5;
%% get roots of S_bar

% % ex1
% delta = 0.4;
% mu_x = 0.4;
% S = @(x)(x+mu_x)./log(x+1) - S_0;

% ex2
r0 = 5;
delta = 0.7;
% S = @(x) r0./(1-delta*exp(-(5-x).*(x-3)));
S = @(x) 1/20 *(r0./(1-delta*exp((x-5).^3.*(x-3))) - S_0);

fplot(S,[2,5])
ns(1) = fzero(S,runner);

while counter ~= 2 && runner <= ns(1)+10
    ns(2) = fzero(S,runner);
    
    if ns(2) >= 1.2*ns(1)
        counter = counter+1;
    end
    runner = runner+step;
end
ns

%% Find min of S

min_S = fminbnd(S,ns(1),ns(2));

%% Evaluate S_bar

counter = 1;

for r = ns(1):0.001:min_S
    for i = ns(2):-0.001:min_S
        if abs(S(r) - S(i)) <= 1/300000
            inex(counter,1) = r;
            inex(counter,2) = i;
            counter = counter+1;
        end
    end
end

%% Evaluate invasion exponent
%
% counter = 1;
%
% for r = ns(1):0.000001:min_S
% for i = ns(2):-0.000001:min_S
% if abs(log(i+1)/(log(r+1)) * (0.4+r) - 0.4 - i) <= 1/10000
% inex(counter,1) = r;
% inex(counter,2) = i;
% ns(2) = i;
% counter = counter+1;
% end
% end
% end

%% plot PIP, inex(:,1) is resident

figure;
hold on
plot(inex(:,1),inex(:,2),'-','Linewidth',2)
plot(inex(:,2),inex(:,1),'-','Linewidth',2)
plot(inex(:,1),inex(:,1),'-','Linewidth',2)
plot(inex(:,2),inex(:,2),'-','Linewidth',2)
hold off
set(gca, 'DataAspectRatio', [1 1 1]);
set(gca, 'XTickLabelMode', 'manual', 'XTickLabel', []);
set(gca, 'YTickLabelMode', 'manual', 'YTickLabel', []);

% ex1
% set(gca, 'XTick', [1,4.24]);
% set(gca, 'YTick', [4.24]);
% xlabel('\rightarrow \gamma_{res}');
% ylabel('\rightarrow \gamma_{inv}');
% text(-0.36,4.24,'\gamma_{max}');
% text(-0.2,-0.2,'\gamma_{min}');
% text(4.24,-0.2,'\gamma_{max}');
% text(0.5,1,'+');
% text(1,1.5,'-');
% text(1.5,1,'+');
% text(1,0.5,'-');
% text(1,-0.2,'\gamma*');

I hope this helps.

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us