Last active
May 22, 2018 12:54
-
-
Save samubernard/587a54c7ce1673590395050f2394d94e to your computer and use it in GitHub Desktop.
FKPP - version individu-centré 1D remplacement aléatoire
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
% FKPP - version individu-centre 1D remplacement aleatoire | |
% Des individus sont distribues sur un intervalle | |
% Certains mutants sont porteurs d'un gene favorable qui donne | |
% un avantage de croissance r | |
% Tous les individus de deplacent selon une marche aleatoire avec | |
% des sauts en sqrt(2*D) | |
% Les wild-types (non-porteurs) on un taux de croissance 0 | |
% A chaque reproduction d'un mutant, l'individu le plus | |
% proche est remplace par la progeniture | |
% Le nombre total d'individu est constant = N0 | |
% la croissance est limitee par le fait la probabilite de remplacer un | |
% wild-type est de 1 - umut, ou umut est | |
% la densite de population de mutant | |
% L'avantage de ce code est qu'il ne depend pas de la densite umut | |
% umut n'est calculee que pour l'affichage | |
% | |
% prob(mutant en x se reproduise) = dt*r | |
% prob(remplacer un wild-type) ~ (1-umut(x)) | |
%% parametre des equations | |
r = 0.5; % taux de croissance | |
D = 0.1; % coefficient de diffusion | |
%% parametres de simulation, espace | |
% x in [x0,x1] | |
h = 0.1; % intervalle de discretisation spatiale | |
S = 20.0; % longueur de l'intervalle | |
x0 = 0; % bord gauche de l'intervalle | |
x1 = S; % bord droit de l'intervalle | |
xi = x0:h:x1; % discretisation | |
J = length(xi); % nombre de points de discretisation | |
%% Individu-centre | |
N0 = 10000; % population totale | |
% chaque individu possede une position x et trait g | |
% g = 1 : gene favorable (mutants) | |
% g = 0 : pas de gene favorable | |
x = x0 + S*rand([1,N0]); % individus distribues sur [x0,x1] | |
g = false(1,N0); | |
g( x < S/10 ) = true; % individus sur la gauche de l'intervalle sont porteurs | |
w = zeros(1,N0); % probabilite pour l'individu de se reproduire | |
%% parametres de simulation, temps | |
t0 = 0; | |
tfinal = 30.0; | |
t = t0; | |
dt = 0.1 | |
%% Difference finie pour comparaison | |
% methode implicite Crank-Nicolson | |
ufd = zeros(J,1); % stocke seulement l'etat au temps t | |
newufd = zeros(J,1); | |
% Le Laplacien discretise est une matrice L de taille JxJ | |
% En 1D elle est symetrique et tridiagonale | |
L0 = sparse(1:J,1:J,-2); % matrice creuse, compacte en memoire | |
L0 = spdiags(ones(J,2),[-1 1],L0); % forme la matrice tridiagonale | |
L0(1,:) = 0; % u(x0) est donne par les conditions au bord | |
L0(J,:) = 0; % u(x1) est donne par les conditions au bord | |
% condition initiales | |
ufd = zeros(J,1); | |
ufd( xi < S/10 ) = 1.0; | |
% Schema implicite Crank-Nicolson | |
A = (speye(J) - dt/h^2*D/2*L0); | |
%% Initialisation | |
figure(1); clf; | |
umut = ksdensity(x(g),xi,'support',[x0,x1]); % kernel density approximation des mutants | |
plot(xi,sum(g)*umut*S/N0,xi,ufd) | |
%% BOUCLE PRINCIPALE | |
tic | |
while t < tfinal | |
drawnow; | |
% on calcule pour chaque mutant une probabilite de se reproduire | |
% (les individus non porteurs ne se reproduisent pas) | |
w = dt*r*g; % probabilite de se reproduire | |
irep = find(rand(1,N0) < w); % realisation | |
% un individu se reproduisant remplace l'individu le plus proche | |
% (autre que lui-meme) | |
% la probabilite de choisir un wild-type est environ 1-umut | |
for i = 1:length(irep) | |
[~,rempl] = min(abs(x - x(irep(i))) + 10*(x == x(irep(i)))); | |
g(rempl) = true; | |
end | |
% on deplace les individu en marche aleatoire | |
x = x + sqrt(dt)*sqrt(2*D)*randn(size(x)); | |
x = abs(x); % condition en x0 reflechissant | |
umut = ksdensity(x(g),xi,'support',[x0,x1]); % kernel density approximation des mutants | |
% solution differences finies implicite pour comparaison | |
newufd =A\(ufd + dt*r*ufd.*(1-ufd) + dt/h^2*D/2*L0*ufd); | |
newufd(1) = newufd(2); | |
newufd(J) = newufd(J-1); | |
ufd = newufd; | |
% affichage | |
plot(xi,sum(g)*umut*S/N0,xi,ufd); | |
axis([x0 x1 0 1.1]); | |
t = t + dt; | |
fprintf("t = %.5f, Nmut = %d\n",t, sum(g)); | |
end | |
toc |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment