Exploring the 2D GrayScott Model with MATLAB

[bolly]

Dear Community,

I was trying to reproduce the 2D GrayScott model given either here: http://blogs.mathworks.com/graphics/2015/03/16/how-the-tiger-got-its-stripes/

or here: http://www.joakimlinde.se/java/ReactionDiffusion/index.php?size=0

The reason was to create a nice screen-saver (so I am not really interested in the mathematics behind). So I tried to use MATLAB for a first trial however my interpretation of the given MATLAB code (first URL) neither gives results as shown in both the urls. I would be very happy if someone could check the following code – maybe I just forget something important:function [t, A, B] = initial_conditions(n)
t = 0;
% Initialize A to one
A = zeros(n); %ones(n);
% Initialize B to zero which a clump of ones
B = zeros(n);
B(50:80 ,50:80) =10;
% A(50:80 ,50:80) =6;
%B(61:80,71:80) = 20;
% --- Executes on button press in pushbutton1.
function pushbutton1_Callback(hObject, eventdata, handles)
figure;
width = 128;
[t, A, B] = initial_conditions(width);

% Setup image
ih=imagesc(B);
set(ih,'cdatamapping','direct')
%colormap(hot);
axis image off;
th=title('');
set(gcf,'position',[80 70 512 512],'color',[1 1 1],'menubar','none')
% Create 'Quit' pushbutton in figure window
uicontrol('units','normal','position',[.45 .02 .13 .07], ...
'callback','set(gcf,''userdata'',1)',...
'fontsize',10,'string','Quit');
f=.55;
k=.062;

da = 1;
db = 0.1;
dt = .001;
done = 0;
count = 0;
while ~done
% anew = A + (da*my_laplacian(A) - A.*B.^2 + f*(1-A))*dt;
% bnew = B + (db*my_laplacian(B) + A.*B.^2 - (k+f)*B)*dt;
%dAdt= A + (da*del2(A) - A.*B.^2 + f*(1-A))*dt;
% dbBdt= B + (db*del2(B) + A.*B.^2 - (k+f)*B)*dt;

dAdt = da*del2(A) - A.*B.*B + f*(1-A);
dBdt= db*del2(B) + A.*B.*B - (k+f)*B;%- k*B;

B = B + dBdt*dt;
A = A + dAdt*dt;

% avoid overflow
if A > 1000
done = 1;
A(60,60)
'over'
end
if A < 5.5000e-005
done = 1;
A(60,60)
'over'
end

if(mod(count,50) == 40)
set(ih,'cdata',A); %,B);
set(th,'string',sprintf('%d %0.2f %0.2f',count,A(60,60),B(60,60)));
drawnow;
end;
count = count + 1;
if ~isempty(get(gcf,'userdata')), done=1;
end % Quit if user clicks on 'Quit' button.
end

 

[Image Analyst]

Mike's blog gives the code. Why not just use his code?

 

[bolly]

Dear Community,

after try- and error-ing a while I finally found a suitable set of parameters which define the most fascinating manifestation of the GrayScott model: The self replicating spots- simulation. However since it takes a lot of calculation this approach is not well suited for my initial goal to make a screen saver out of it– however if someone else likes to experiment himself a little bit – here is a working code snipped (copied from here and there - majorly from Mikes Blog).

function [t, A, B] = initial_conditions(n)
t = 0;
% Initialize A to one
A = ones(n);
% Initialize B to zero which a clump of ones
B = zeros(n);
B(51:60 ,51:70) = 1;
B(61:80,71:80) = 1;

function out = my_laplacian(in)
out = -in ...
+ .20*(circshift(in,[ 1, 0]) + circshift(in,[-1, 0]) ...
+ circshift(in,[ 0, 1]) + circshift(in,[ 0,-1])) ...
+ .05*(circshift(in,[ 1, 1]) + circshift(in,[-1, 1]) ...
+ circshift(in,[-1,-1]) + circshift(in,[ 1,-1]));

% --- Executes on button press in pushbutton1.
function pushbutton1_Callback(hObject, eventdata, handles)f = 0.02;
k = 0.079;% Diffusion rates
da = 1;
db = 0.4;%0.5;%.5;
% Size of grid
width = 128;

dt = .25; %.15
[t, A, B] = initial_conditions(width);
figure;
% Setup image
ih=imagesc(B);

set(ih,'cdatamapping','direct')
%colormap(hot);
axis image off;
th=title('');
set(gcf,'position',[80 70 512 512],'color',[1 1 1],'menubar','none')
% Create 'Quit' pushbutton in figure window
uicontrol('units','normal','position',[.45 .02 .13 .07], ...
'callback','set(gcf,''userdata'',1)',...
'fontsize',10,'string','Quit');
% Add a text object to show the current time.
ht = text(3,width-3,'Time = 0');
ht.Color = [.95 .2 .8];
done = 0;

%count = 0;
while ~done
anew = A + (da*my_laplacian(A) - A.*B.^2 + f*(1-A))*dt;
bnew = B + (db*my_laplacian(B) + A.*B.^2 - k*B) *dt; %;- (k+f)*B)*dt;
A = anew;
B = bnew;
t = t+dt;
if(mod(count,100) == 40)
set(ih,'cdata',B.*30); %,B); %multiplying B with a factor was important for visualization
drawnow;
end;
% count = count + 1;
if ~isempty(get(gcf,'userdata')), done=1;
end % Quit if user clicks on 'Quit' button.
end

@ Image analyst: I certainly tried Mikes code but using my old Matlab 6.5 there were no waves visible (probably because ‘image’ has no autoscale feature). Furthermore I wanted ‘Selfreplicating spots’ shown on Joakim Lindes’ page ( – however without Mikes great Blog I wouldn’t be able to run a GrayScott model myself).