Equation FKPP difference finie 2D explicite

FKPP_2D_fd_explicite.m

% Equation FKPP difference finie 2D explicite
% FKPP equation finite-difference explicite 2D

% parametre des equations/equation parameters
r = 0.5; % taux de croissance
D = 0.1; % coefficient de diffusion

% parametres de simulation, espace/space simulation parameters
% x in Omega [0,Sx]x[0,Sy]
h = 0.025;  % intervalle de discretisation spatiale/space step
Sx = 10.0;
Sy = 6.0;
x0 = 0;     % bord gauche de l'intervalle/left border
x1 = Sx;    % bord droit de l'intervalle/right border
y0 = 0;
y1 = Sy;
x = x0:h:x1; % discretisation/discretization
y = y0:h:y1; 
[X,Y] = meshgrid(x,y);
J1 = length(x); 
J2 = length(y);
J  = J1*J2;  % nombre de points de discretisation/number of points

% variable dynamiques/dynamical variables
u = zeros(J,1); % stocke seulement l'etat au temps t
newu = zeros(J,1);

% discretisation du Laplacien avec conditions aux bord de Neumann
% no flux, i.e. du/dx = 0 en x0 et en x1
% Neumann, no flux boundary conditions

% Le Laplacien discretise est une matrice L de taille JxJ
% Discretized Laplacian, 

cornertopleft = 1;
cornerbottomleft = J2;
cornertopright = (J1-1)*J2+1;
cornerbottomright = J;
sideleft = 2:J2-1;
sidetop = J2+1:J2:J2*(J1-2)+1;
sidebottom = 2*J2:J2:J2*(J1-1);
sideright = J2*(J1-1)+2:J-1;
side = [cornertopleft, cornertopright, cornerbottomleft, cornerbottomright, ...
    sideleft, sidetop, sidebottom, sideright];
interieur = setdiff(1:J, side);

% interieur/inner domain
L = sparse(interieur,interieur,-4,J,J); % matrice creuse, compacte en memoire
L = L + sparse(interieur,interieur+1,1,J,J);
L = L + sparse(interieur,interieur-1,1,J,J);
L = L + sparse(interieur,interieur+J2,1,J,J);
L = L + sparse(interieur,interieur-J2,1,J,J);

% condition initiales/initial condtions
% u( X(:).^2 + Y(:).^2 < 3 ) = 0.8; % condition initiale constante par
% morceaux
u = exp( - X(:).^2/0.5 - Y(:).^2/2 ); % condition initiale lisse

% parametres de simulation, temps/time simulation parameters
t0 = 0;
tfinal = 20; 
t = t0;
tp = t0;
% must have k < h^2/4/D
k = min(0.5,0.99*( h^2/4/D ));

% figures
figure(1); clf;
image(x,y,255*reshape(u,J2,J1));
axis equal;

% uncomment the next two lines to pause the simulation at the initial condition
% disp('appuyer sur une touche pour continuer');
% pause

% BOUCLE PRINCIPALE/MAIN LOOP
while t < tfinal
    drawnow;
    newu = u + k*r*u.*(1-u) + k/h^2*D*L*u;  % time-stepping
    newu(sideleft) = newu(sideleft+J2);       % no flux - left
    newu(sideright) = newu(sideright-J2);     % no flux - right
    newu(sidetop) = newu(sidetop+1);          % no flux - top
    newu(sidebottom) = newu(sidebottom-1);    % no flux - bottom
    % There is no specific rule for updating the corners, but they
    % are not strictly part of the discretized boundary

    u = newu;
    
    if ( fix(10*t) > fix(10*tp) )
        image(x,y,255*reshape(u,J2,J1));
        axis equal;
    end
    tp = t;
    t = t + k;
    fprintf("t = %.5f\n",t);
end