2D reaction-diffusion equation, explicit method

reaction_diffusion.m

%% 2D reaction-diffusion equation, explicit method

% equation parameters
r = 1.0; % some model parameter
D = 0.1; % diffusion coefficient

% Simulation parameters: space
% Domain: rectangle of size Sx by Sy
h = 0.2; % grid step size
Sx = 10.0;
Sy = 6.0;
x0 = 0;
x1 = Sx;
y0 = 0;
y1 = Sy;
x = x0:h:Sx; % x-space discretisation
y = y0:h:Sy; % y-space discretisation
[X,Y] = meshgrid(x,y); % build grid for display and initial conditions
J1 = length(x);
J2 = length(y);
J  = J1*J2; % total number of grid points

% dynamical (dependent) variables
u = zeros(J,1); % initialize the solution vector u to zero
unew = zeros(J,1); % initialize auxiliary solution vector

% Discretisation of the Laplacian
% L is a JxJ matrix (J is large!!)

% first define the borders 
% The grid is indexed with point (1,1) in upper left 
% corner and (J2,J1) in lower right corner:
% 1  J2+1  .  .  .  (J1-1)*J2+1
% 2  J2+2           (J1-1)*J2+2
% .  .              . 
% .  .              . 
% J2 2*J2           J1*J2

upperleftcorner = 1;
lowerleftcorner = J2;
upperrightcorner = (J1-1)*J2+1;
lowerrightcorner = J; % = J1*J2
leftborder = 2:J2-1;
rightborder = J2*(J1-1)+2:J-1;
upperborder = J2+1:J2:J2*(J1-2)+1;
lowerborder = 2*J2:J2:J2*(J1-1);
border = [upperleftcorner, lowerleftcorner, ...
    upperrightcorner, lowerrightcorner, ...
    leftborder, rightborder, upperborder, lowerborder];
interior = setdiff(1:J, border);

% L matrix
L = sparse(interior,interior, -4, J,J); % sparse matrix with -4 on the diagonal entries corresponding to interior points, make it JxJ
L = L + sparse(interior,interior+1,1,J,J);
L = L + sparse(interior,interior-1,1,J,J);
L = L + sparse(interior,interior+J2,1,J,J);
L = L + sparse(interior,interior-J2,1,J,J);
% all non-interior rows are zeros by default

% Initial conditions
% u = exp( - X(:).^2/0.5 - Y(:).^2/2 ); % smooth init conds

% time parameters
t0 = 0;
tfinal = 5;
t = t0;
dt = min(0.5, 0.9*( h^2/4/D ));


figure(1); clf;
image(x,y,255*reshape(u,J2,J1)); % reshape u to matrix, and plot image on a 0, 255 RGB scale
axis equal
pause

% MAIN LOOP
tic % time the simulation
while t < tfinal
    % update the reaction diffusion eq
    newu = u + dt*(-r)*((X(:)-1).^2 + (Y(:)-2).^2 < 1 ) + dt/h^2*D*L*u;
    

    % boundary conditions: Neuman no flux du/dt = 0
    newu(leftborder) = newu(leftborder+J2);
    newu(rightborder) = newu(rightborder-J2);
    newu(upperborder) = newu(upperborder+1);
    newu(lowerborder) = newu(lowerborder-1);
    
    u = newu;
    
    image(x,y,255*abs(reshape(u,J2,J1)));
    axis equal
    drawnow;
    
    t = t + dt;
    fprintf("t = %.5f\n", t);
end
toc

% the end