% This code will create the plots that were used to make the animations shown in Turing II.
% 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.
% MATLAB code adapted by Michael Ralph for usability in Octave.
clear all
close all
% Set-up
% Parameter values
bc = 0.35;
Du = 1; Dv = 20; % Diffusion constants
% Grid and initial data:
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:300;
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.5*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
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
subplot(2,3,1:3)
plot(x,v);
axis([-1 w 5.01 15.01])
title(['t = ', num2str(j*dt)],'fontsize',24)
% Save for surface
U(j,:) = u;
V(j,:) = v;% space-time stripchart
subplot(2,3,4:6)
imagesc(x, 1:Nx, u');
drawnow()
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');
colorbar; % View and set current colormap