Crank-Nicholson in 2D with MATLAB

• MATLAB
I have the code which solves the Sel'kov reaction-diffusion in MATLAB with a Crank-Nicholson scheme.

I would love to modify or write a 2D Crank-Nicolson scheme which solves the equations:

##u_t = D_u(u_{xx}+u_{yy})-u+a*v+u^2*v##
##v_y = D_v(v_{xx}+v_{yy}) +b-av-u^2v##

Where ##D_u, D_v## are the diffusive constants and ##a ,b## are just positive constants.

Matlab:
% Modelo de Sel'kov en 1D

% Este código en MATLAB simulara el modelo de Sel'kov para la glucólisis en 1D
clear all
close all

%%%%%%%%%%%%%%%%%%%%%
%%%Inicialización %%%
%%%%%%%%%%%%%%%%%%%%%

% Parámetros
bc  = 0.6; %0.98; %0.90033; %0.9502; % %0.6;
ac = 0.08; %0.08; %0; %
Dv = 0.001; %0.02; %0.002;
Du = 0.001;  %Constantes difusivas

% w = 10;  %sin patrón
w = 1;  % patrón

Nx = 1000; % Total de puntos discretizados en el dominio [0,L]
x = linspace(0,w,Nx);
dx = x(2) - x(1);
%pause();

dt = 0.0001; % tamaño del paso
t = 0:dt:100;
Nt = length(t); %numero de puntos en el tiempo

% Condiciones para la superficie
[X, T] = meshgrid(x, t);
%U = 0*X;
%V = 0*X;

% Vectores columna (más fáciles)
x = x(:);
t = t(:);

%Condición inicial: pequeña perturbación lejos del estado estable
u =  bc.*ones(length(x),1) + 0.5.*rand(Nx, 1); %0.01.*(cos(pi.*x)) +    %exp(-(x-0.5).^2);
v =  (bc/(ac + bc^2)).*ones(length(x),1) + 0.5.*rand(Nx,1); %0.01.*(cos(2.*pi.*x)) +  % exp(-(x-0.5).^2);
U(1,:) = u;
V(1,:) = v;

%%%%%%%%%%%%%%%%%%%%%%%%%
%%%Generando la matriz%%%
%%%%%%%%%%%%%%%%%%%%%%%%%
% Usando un esquema implícito
%  Para u
lambda1 = dt.*Du./(2.*dx.^2);
Au = eye(Nx).*(1+2.*lambda1);
Au(1+1:Nx+1:end) = -lambda1;
Au(Nx+1:Nx+1:end) = -lambda1;
Au([1,end]) = 1+lambda1;
Bu = eye(Nx).*(1-2.*lambda1);
Bu(1+1:Nx+1:end) = lambda1;
Bu(Nx+1:Nx+1:end) = lambda1;
Bu([1,end]) = 1-lambda1;

%  Para v
lambda2 = dt.*Dv./(2.*dx.^2);
Av = eye(Nx).*(1+2.*lambda2);
Av(1+1:Nx+1:end) = -lambda2;
Av(Nx+1:Nx+1:end) = -lambda2;
Av([1,end]) = 1+lambda2;
Bv = eye(Nx).*(1-2.*lambda2);
Bv(1+1:Nx+1:end) = lambda2;
Bv(Nx+1:Nx+1:end) = lambda2;
Bv([1,end]) = 1-lambda2;

%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%

axis([0 1 0 1])
uma(1)=u(Nx/2);
vma(1)=v(Nx/2);
for j = 1:Nt-1
% f y g son los términos no lineales en el modelo de Sel'kov

f= -u+ac.*v+u.^2.*v;
g= bc-ac.*v-u.^2.*v;
%f = u.^2./v-bc*u;
%g = u.^2 - v;
%size(Bu), size(u),
% En cada paso resolvemos el sistema de ecuaciones
u = Au\(Bu*u + dt.*f); %u = Bu\(u + dt.*f);
v = Av\(Bv*v + dt.*g); %v = Bv\(v + dt.*g);

uma(j+1)=u(Nx/2);
vma(j+1)=v(Nx/2);
%u = up;
%v = vp;

% Gráficas
plot(x,u,'g.-', 'linewidth',1);
hold on;
plot(x,v,'r.-', 'linewidth',1);
hold off;
axis([0 1 0 2])
title(['t = ', num2str(j*dt)],'fontsize',24)
drawnow;
%pause();
% Datos para la superficie
U(j+1,:) = u;
V(j+1,:) = v;
end
%%%% Gráfica de la superficie %%%%
figure(1);
s = surf(x, t, U);
set(s, 'EdgeColor', 'none', 'FaceColor', 'interp');
xlabel('<----x---->')
ylabel('<----t---->')
zlabel('<----u---->')
figure(2);
s = surf(x, t, V);
set(s, 'EdgeColor', 'none', 'FaceColor', 'interp');
xlabel('<----x---->')
ylabel('<----t---->')
zlabel('<----v---->')

%%%% contour plot %%%
figure(3);
p = pcolor(x, t, U);
set(p, 'EdgeColor', 'none', 'FaceColor', 'interp');

figure(4);
q = pcolor(x, t, V);
set(q, 'EdgeColor', 'none', 'FaceColor', 'interp');

figure(5);
plot(t,uma,'g.-', t,vma,'r.-');
legend('ATP','ADP')

I don't know how to modify the code, but I do totally understand how to use the Crank-Nicolson scheme.

POST adjusted by mentor using double # around latex expressions instead of single \$.

Last edited by a moderator: