Equation FKPP difference finie 2D implicite Crank-Nicolson

FKPP_2D_fd_implicite.m

% Equation FKPP difference finie 2D implicite Crank-Nicolson
% FKPP equation finite-difference implicite Crank-Nicolson 2D

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

% parametres de simulation, espace
% 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 end of the interval
x1 = Sx;  % bord droit de l'intervalle/right end of the interval
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/store only current sol
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
L0 = sparse(interieur,interieur,-4,J,J); % matrice creuse, compacte en memoire
L0 = L0 + sparse(interieur,interieur+1,1,J,J);
L0 = L0 + sparse(interieur,interieur-1,1,J,J);
L0 = L0 + sparse(interieur,interieur+J2,1,J,J);
L0 = L0 + sparse(interieur,interieur-J2,1,J,J);

% condition initiales
% 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
t0 = 0;
tfinal = 35; 
t = t0;
tp = t0;
k = 0.5;

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

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

% Schema implicite Crank-Nicolson implicit scheme
A = (speye(J) - k/h^2*D/2*L0);

% BOUCLE PRINCIPALE/MAIN LOOP
while t < tfinal
    drawnow;
    newu =A\(u + k*r*u.*(1-u) + k/h^2*D/2*L0*u);
    newu(sideleft) = newu(sideleft+J2);
    newu(sideright) = newu(sideright-J2);
    newu(sidetop) = newu(sidetop+1);
    newu(sidebottom) = newu(sidebottom-1);
    u = newu;
    if ( fix(10*t) > fix(10*tp) )
        image(x,y,255*reshape(u,J2,J1));
    end
    tp = t;
    t = t + k;
    fprintf("t = %.5f\n",t);
end