% The Mathematics of Patterns 2013 % Turing Instabilities II: Gierer Meinhardt system in 1D % This MATLAB code will show you how to numerically simulate the 1D Gierer % Meinhardt system, and make the (simplified) version of the animations in % the Turing Instabilities II section of the website. clear all close all %%%%%%%%%%%%%%%%%%%%% %%% Set-up %%% %%%%%%%%%%%%%%%%%%%%% % Parameter values bc = 0.35; Du = 1; Dv = 30; % Diffusion constants % Grid and initial data: % w = 10; %no pattern w = 80; % pattern Nx = 500; % How many points we want to discretize our domain with x = linspace(0,w, Nx); dx = x(2) - x(1); dt = 1; % size of our time step t = 0:dt:400; Nt = length(t); % Number of time points % Set up for the surface [X, T] = meshgrid(x, t); U = 0*X; V = 0*X; % Easier to deal with column vectors x = x(:); t = t(:); %Initial conditions: small perturbation away from steady state u = 1/bc*ones(length(x),1) + 0.01*rand(Nx, 1); v = 1/bc^2*ones(length(x),1); % Save initial conditions U(1,:) = u; V(1,:) = v; %%%%%%%%%%%%%%%%%%%%%%%%% %%% Making the matrix %%% %%%%%%%%%%%%%%%%%%%%%%%%% % To begin, let us recall how to set up the matrices used in the explicit % and implicit finite difference methods. %%% Forward (explicit) method %%% % We want a tridiagonal matrix (see notes for details) a = (1-2*Du*dt/dx^2); % values along the diagonal b = Du*dt/dx^2; % values in the off-diagonal main = a*sparse(ones(Nx,1)); off = b*sparse(ones(Nx-1,1)); Bu = diag(main) + diag(off,1) + diag(off,-1); %Still a sparse matrix % Satisfying boundary conditions Bu(1, end-1) = b; Bu(end, 2) = b; % To have a more numerically stable code, we use the implicit method. %%% Backward (implicit) method %%% % For u a = (1+2*Du*dt/dx^2); % values along the diagonal b = Du*dt/dx^2; % values in the off-diagonal main = a*sparse(ones(Nx,1)); off = -b*sparse(ones(Nx-1,1)); Bu = diag(main) + diag(off,1) + diag(off,-1); %Still a sparse matrix % Satisfying boundary conditions Bu(1, end-1) = -b; Bu(end, 2) = -b; % Same thing for v a = (1+2*Dv*dt/dx^2); b = Dv*dt/dx^2; main = a*sparse(ones(Nx,1)); off = -b*sparse(ones(Nx-1,1)); Bv = diag(main) + diag(off,1) + diag(off,-1); Bv(1, end-1) = -b; Bv(end, 2) = -b; %%%%%%%%%%%%%%%%%%%%%%%%% %%% Plotting %%% %%%%%%%%%%%%%%%%%%%%%%%%% figure(1); %create new figure plot(x,u,'g.-', 'linewidth',1); hold on; plot(x,v,'r.-', 'linewidth',1); hold off; axis([-1 80 -.01 15.01]) % Fix axis limits for j = 1:Nt % f and g are the reaction terms in the G-M system f = u.^2./v-bc*u; g = u.^2 - v; % At each step we need to solve the system u = Bu\(u + dt*f); % backward Euler v = Bv\(v + dt*g); % Plot plot(x,u,'g.-', 'linewidth',1); hold on; plot(x,v,'r.-', 'linewidth',1); hold off; axis([-1 80 -.01 15.01]) title(['t = ', num2str(j*dt)],'fontsize',24) drawnow; % Save for surface U(j,:) = u; V(j,:) = v; end %%%% Plotting the surface %%%% figure(2); s = surf(x, t, U) set(s, 'EdgeColor', 'none', 'FaceColor', 'interp'); % Sets up the colors xlabel('x') ylabel('t') zlabel('u') %%%% contour plot %%% figure(3); p = pcolor(x, t, U); set(p, 'EdgeColor', 'none', 'FaceColor', 'interp');