# Crank-Nicholson in 2D with MATLAB

• MATLAB
• Aldo Leal
In summary, the conversation was about the Sel'kov reaction-diffusion model and the desire to modify or write a 2D Crank-Nicolson scheme to solve the equations. The code provided is a MATLAB simulation of the Sel'kov model in 1D, with parameters and initial conditions specified. The code uses an implicit scheme and solves for both the ATP and ADP concentrations over time. Several graphs and contour plots are generated to visualize the results. The conversation concludes with the mention of MATLAB tutorials for solving PDEs.

#### 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.

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:

## 1. What is Crank-Nicholson in 2D with MATLAB?

Crank-Nicholson in 2D with MATLAB is a numerical method used to solve partial differential equations (PDEs) in two dimensions. It combines the Crank-Nicholson method, which is a finite difference method, with the power of MATLAB programming language.

## 2. How does Crank-Nicholson in 2D with MATLAB work?

This method works by discretizing the PDE into a grid and approximating the derivatives using finite differences. It then iteratively solves the equations at each grid point using the Crank-Nicholson scheme, which is a combination of the backward Euler and forward Euler methods. This results in a more accurate and stable solution compared to using either method alone.

## 3. What are the advantages of using Crank-Nicholson in 2D with MATLAB?

One advantage is that it can handle stiff and non-stiff equations, making it suitable for a wide range of PDEs. It also has a higher accuracy and stability compared to other numerical methods, such as the explicit Euler or Runge-Kutta methods. Additionally, the use of MATLAB allows for easy implementation and customization of the code.

## 4. What types of PDEs can be solved using Crank-Nicholson in 2D with MATLAB?

This method can be applied to a variety of PDEs, including diffusion equations, convection-diffusion equations, and reaction-diffusion equations. It is also suitable for solving PDEs in various fields such as fluid dynamics, heat transfer, and population dynamics.

## 5. Are there any limitations to using Crank-Nicholson in 2D with MATLAB?

One limitation is that it requires the PDE to be linear and have constant coefficients. It may also be computationally expensive for large grids or when solving time-dependent problems. Additionally, the accuracy of the solution may be affected by the choice of grid size and time step.