Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

I Crank-Nicholson in 2D with MATLAB

  1. May 23, 2016 #1
    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.

    Code (Matlab M):
    % 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 $.
     
    Last edited by a moderator: May 23, 2016
  2. jcsd
  3. May 23, 2016 #2

    jedishrfu

    Staff: Mentor

Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?
Draft saved Draft deleted



Similar Discussions: Crank-Nicholson in 2D with MATLAB
  1. Crank Nicholson scheme (Replies: 0)

  2. Crank-Nicolson method (Replies: 14)

Loading...