Created
May 22, 2018 12:56
-
-
Save samubernard/5773316c09bc019af95ec6fb493ba47a to your computer and use it in GitHub Desktop.
FKPP - version individu-centre 1D croissance nonlinaire
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 croissance nonlineaire | |
% 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 probabilite qu'un mutant remplace un wild-type est modulee | |
% par (1 - umut), ou u est la densite de mutant: | |
% | |
% prob(mutant en position x remplace un wild-type) = dt*r*(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) | |
uxg = ksdensity(x(g),x); | |
w = dt*r*(1-sum(g)*uxg*S/N0).*g; % probabilite de se reproduire | |
irep = find(rand(1,length(w)) < w); % realisation | |
% un individu se reproduisant remplace l'individu wild-type le plus | |
% proche | |
for i = 1:length(irep) | |
[~,rempl] = min(abs(x - x(irep(i))) + 100*g); | |
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