% 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.
%%%%%%%
% This Matlab Script was sourced from an original script written by
% Yayoi Teramoto Kimura. The Matlab script can be found in his report
% titled "The Mathematics of Patterns The modeling and analysis of
% reaction-diffusion equations."
%%%%%%%
% Associated curriculum and script organization developed by Abdalkader Alsabawy under supervision of Mr. Michael Ralph
% working with the STEM Learning Center at the University of Kansas.
%%%%%%%
clear all
close all
% Set-up
% Parameter values
bc = 0.35;
Du = 1; Dv = 20; % Diffusion constants
% Grid and initial data:
w = 80; % Width of the space being simulated (in units)
Nx = 500; % Number of points discretized in our domain
x = linspace(0,w, Nx);
dx = x(2) - x(1);
dt = 1; % size of our time step
t = 0:dt:800;
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)
hi = image(v);
hi.CDataMapping = 'scaled';
hi.CData = v;
drawnow limitrate
subplot(2,3,2:3)
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)
% 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