Crank-Nicholson in 2D with MATLAB

[Aldo Leal]

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.

% 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

% Información inicial y mallado
% 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);
% Condiciones iniciales salvadas.
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;%%%%%%%%%%%%%%%%%%%%%%%%%
%%%    Resultados     %%%
%%%%%%%%%%%%%%%%%%%%%%%%%

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;
     legend('ATP','ADP')
     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 $.

 

[jedishrfu]

MATLAB has several tutorials on solving PDE's. Here is one such reference I found:

http://www.mathworks.com/help/matlab/math/partial-differential-equations.html