Reaction Diffusion model leading to pattern formtaion

This model is usually referred to as the Scott-Grey model in the literature. Here I use a Crank-Nicolson scheme for the diffusion and RK4 along with Strang Splitting for the reaction. The RK4 is homemade and is not the built in Matlab version. For more information look for Reaction-Diffusion models.

Parameters

Grid and constants.

N=256
z=linspace(-1,1,N+1);
z=z(1:end-1);
L=2.5;
[x,y]=meshgrid(L*z,L*z);
t=0;
f=0.038;
kpar=0.06;
ru=2e-5;
rv=1e-5;
N =

256

RK4 functions.

f1=@(u,v) -u.*v.*v+f*(1-u);
f2=@(u,v) u.*v.*v-(f+kpar)*v;

Initial conditions.

r2=x.*x+y.*y;
u0=ones(size(x)); v0=zeros(size(x));
u0=u0-(0.5).*exp(-(r2.^8)/((0.25*L)^8))+0.02*(rand(size(r2))-0.5);
v0=v0+(0.254).*exp(-(r2.^8)/((0.25*L)^8))+0.02*(rand(size(r2))-0.5);
u=u0;
v=v0;

Spectral parameters.

dk=pi/L;
ksvec(1)=0;
for ii=2:(N/2)
ksvec(ii)=ii-1;
ksvec(N/2+ii)=-N/2 + ii -1;
end
ksvec=ksvec*dk;
[k2,l2]=meshgrid(ksvec.^2,ksvec.^2);
ksvec(N/2+1)=0;
[k,l]=meshgrid(ksvec,ksvec);
k2=k.*k; l2=l.*l;
lapu=-ru*(k2+l2);
lapv=-rv*(k2+l2);

Time and output parameters.

dt=0.2;
dto2=dt/2;
dt2=dt*dt;
numsteps=500;
numouts=5;

Plot ICs

figure(1)
clf
set(gcf,'DefaultLineLineWidth',3,'DefaultTextFontSize',12,...
'DefaultTextFontWeight','bold','DefaultAxesFontSize',12,...
'DefaultAxesFontWeight','bold');
subplot(2,1,1)
pcolor(x,y,u); shading interp
xlabel('x','fontweight','bold','fontsize',12);
ylabel('y','fontweight','bold','fontsize',12);
colorbar
mystr=['time = ' num2str(t)];
title(mystr);
subplot(2,1,2)
pcolor(x,y,v); shading interp
xlabel('x','fontweight','bold','fontsize',12);
ylabel('y','fontweight','bold','fontsize',12);
colorbar
drawnow

Set up Strang splitting.

myfactue=dt*lapu;
myfactve=dt*lapv;
myfactu=(1+dto2*lapu)./(1-dto2*lapu);
myfactv=(1+dto2*lapv)./(1-dto2*lapv);

Main Loop

for ii=1:numouts
for jj=1:numsteps
t=t+dt;

Strang split two steps at a time reversing order.

u=ifft2(myfactu.*fft2(u));
v=ifft2(myfactv.*fft2(v));

RK4.

rkk1=f1(u,v);
rkl1=f2(u,v);
rkk2=f1(u+0.5*dt*rkk1,v+0.5*dt*rkl1);
rkl2=f2(u+0.5*dt*rkk1,v+0.5*dt*rkl1);
rkk3=f1(u+0.5*dt*rkk2,v+0.5*dt*rkl2);
rkl3=f2(u+0.5*dt*rkk2,v+0.5*dt*rkl2);
rkk4=f1(u+dt*rkk3,v+0.5*dt*rkl3);
rkl4=f2(u+dt*rkk3,v+0.5*dt*rkl3);
kforrk=(dt/6)*(rkk1+2*rkk2+2*rkk3+rkk4);
lforrk=(dt/6)*(rkl1+2*rkl2+2*rkl3+rkl4);

First part.

u=u+kforrk;
v=v+lforrk;

Second step in reverse order now.

t=t+dt;
% RK4
rkk1=f1(u,v);
rkl1=f2(u,v);
rkk2=f1(u+0.5*dt*rkk1,v+0.5*dt*rkl1);
rkl2=f2(u+0.5*dt*rkk1,v+0.5*dt*rkl1);
rkk3=f1(u+0.5*dt*rkk2,v+0.5*dt*rkl2);
rkl3=f2(u+0.5*dt*rkk2,v+0.5*dt*rkl2);
rkk4=f1(u+dt*rkk3,v+0.5*dt*rkl3);
rkl4=f2(u+dt*rkk3,v+0.5*dt*rkl3);
kforrk=(dt/6)*(rkk1+2*rkk2+2*rkk3+rkk4);
lforrk=(dt/6)*(rkl1+2*rkl2+2*rkl3+rkl4);

Second part.

u=u+kforrk;
v=v+lforrk;

Diffusion.

u=ifft2(myfactu.*fft2(u));
v=ifft2(myfactv.*fft2(v));
end

Plot

figure(1)
clf
set(gcf,'DefaultLineLineWidth',3,'DefaultTextFontSize',12,...
'DefaultTextFontWeight','bold','DefaultAxesFontSize',12,...
'DefaultAxesFontWeight','bold');
subplot(2,1,1)
pcolor(x,y,u); shading interp
xlabel('x','fontweight','bold','fontsize',12);
ylabel('y','fontweight','bold','fontsize',12);
colorbar
mystr=['time = ' num2str(t)];
title(mystr);
subplot(2,1,2)
pcolor(x,y,v); shading interp
xlabel('x','fontweight','bold','fontsize',12);
ylabel('y','fontweight','bold','fontsize',12);
colorbar
drawnow
end