% 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');