% solving the Nash equilibrium for multiple DVM strategies of predator and prey % code Jérôme Pinti - associated with AmNat paper "Predator prey games in % multiple habitats reveal mixed strategies in diel vertical migration." % parameters clear; n = 30; % the number of levels H = 300; % [m] total depth dH = H/n; % [m] depth bin n cN = 1E-5; % [# m^-1 day^-1] cost of migrating 1 m for prey cP = 1E-5; % [# m^-1 day^-1] cost of migrating 1 m for predator gmax = 0.1; % [# day^-1] maximum sigma = 0.65; % day light fractio growth rate for prey b1 = linspace(-2,2,41); bmaxA = 10.^b1; %[m3 day-1] maximum clearance rate for predator Pavg = 1; % [# m^-3] average density of predators Navg = 10000; % [# m^-3] average density of prey mu = 0.01; % [day^-1] mean density dependent mortality of predator epsilon = 0.01; % [] how much predator reproduce on 1 prey rho = 10^-3; % reduction of voracity from day to night k = 0.07; % [m^-1] light attenuation coefficient zo = 50; zm = 10; L0 = 1; % [W/m^2] half saturation light intesity Lmax = 100; % [W/m^2] surface max day light intensity Niter = 2000000; %decrease the number for a coarser (but faster) output Iavg = 1000; dtfact = 0.1; fact = 1; % Initialization and equations l = @(z) exp(k*z); % functional form of light gfun = @(z) (1 - tanh((-zo - z)/zm))/2; % functional form of growth bfun = @(z) Lmax*l(z)./(L0 + Lmax*l(z)); % functional form of clearance rate ofun = @(z) rho*Lmax*l(z)./(L0 + rho*Lmax*l(z)); %functional form of clearance rate at night Vi = (1:n); %Water layers Ai = Vi'*ones(1,n); zi = -Vi*dH + dH/2; NdayA = zeros(n,length(bmaxA)); NnigA = zeros(n,length(bmaxA)); PdayA = zeros(n,length(bmaxA)); PnigA = zeros(n,length(bmaxA)); N = ones(n,n); Nsum = sum(sum(N)); N = N/Nsum; P = N; FNmaxA = zeros(size(bmaxA)); FPmaxA = zeros(size(bmaxA)); FNmaxM = zeros(size(bmaxA)); FPmaxM = zeros(size(bmaxA)); FNmaxV = zeros(size(bmaxA)); FPmaxV = zeros(size(bmaxA)); FNmax0 = zeros(size(bmaxA)); FPmax0 = zeros(size(bmaxA)); FNmax1 = zeros(size(bmaxA)); FPmax1 = zeros(size(bmaxA)); FNmaxIavg = zeros(1,Iavg); FPmaxIavg = zeros(1,Iavg); NIavg = zeros(1,Iavg); PIavg = zeros(1,Iavg); FNIavg = zeros(1,Iavg); FPIavg = zeros(1,Iavg); FNmaxIavgA = zeros(length(bmaxA),Iavg); FPmaxIavgA = zeros(length(bmaxA),Iavg); resetNP = true; saveall = 0; % Solver for j = 1:length(bmaxA) bmax = bmaxA(j); if resetNP N = ones(n,n); Nsum = sum(sum(N)); N = N/Nsum; P = N; end v = bmax*bfun(zi); %[m3 day^-1] o = bmax*ofun(zi); %[m3 day^-1] g = gmax*gfun(zi); %[day^-1] Nday = sum(N,1); Nnig = sum(N,2); Pday = sum(P,1); Pnig = sum(P,2); CN = cN*dH*abs(Ai-Ai'); GN = (1-sigma)*g'*ones(1,n) + sigma*ones(n,1)*g - CN ; MN = n*Pavg*((1-sigma)*(o'.*Pnig)*ones(1,n) + sigma*ones(n,1)*(Pday.*v)) ; CP = cP*dH*abs(Ai-Ai'); GP = epsilon*n*Navg*((1-sigma)*(o'.*Nnig)*ones(1,n) + sigma*ones(n,1)*(Nday.*v)) - CP; MP = mu*n*Pavg*((1-sigma)*Pnig*ones(1,n) + sigma*ones(n,1)*Pday) ; FN = GN - MN; FP = GP - MP; i = Niter; notdone = true; dNmax0 = 1E-11; dPmax0 = 1E-11; Ncum = zeros(n,n); Pcum = Ncum; while notdone i = i - 1; %Nold = N; Pold = P; %FNold = FN; FPold = FP; FNmax = max(max(FN)); FNmin = min(min(FN)); FPmax = max(max(FP)); FPmin = min(min(FP)); fact = dtfact/max([FPmax,FNmax,-FNmin,-FPmin]); % fact is dynamic so that maximum increase is at most 20% per time step N = N + fact*FN.*N; P = P + fact*FP.*P; N(N 0); end % Ncum = Ncum / Iavg; Pcum = Pcum / Iavg; N(N