%
% Ricardo Carretero 20/Oct/2003
% Integrates a SYSTEM of 2 reaction diffusion eqs in 1D using finite differences
%
% u_t = f(u,v) + D1 u_xx
% v_t = g(u,v) + D2 v_xx
%
% Ex: PDE_2pop(N,dt,tf,L,R,disp,savesnaps)
%
% npts=N, dt:time step, tf:t_final, domain:[L,R]
%
%Ex:
%global aa bb dx dxx D1 D2 oos allt allx allu allv;
%N=75;dt=0.002;tf=200;L=0;R=100;snaps=100;disp=100;
%aa=1;bb=3.65;D1=1;D2=335;
%aa=1;bb=3.65;D1=1;D2=345;
%PDE_2pop(N,dt,tf,L,R,disp,snaps)
%wave_plotsurf(allt,allx,allu',R,1,N+1,'no')
%wave_plotsurf(allt,allx,allv',R,1,N+1,'no')
%
function [] = PDE_2pop(N,dt,tf,L,R,exposh,saveframes)
global allt allx allu allv;
global dx dxx D1 D2 oos u0 v0;
global aa bb;
ran=1;
noise=0;
maxstep=tf/dt;
dispscreen= fix(maxstep/10);
if(exposh>0) stopdisp= fix(maxstep/exposh);
else stopdisp=9e+99;
end
if(saveframes>0)stopsave= fix(maxstep/saveframes);
else stopsave=9e+99;
end
% Spatial computational domain defined
dx=(R-L)/N;
x=L+(0:N-1)*dx;
dxx = 1/dx^2;
oos=1/6;
rand('state',0) %this fixes the seed for the random generation
strlength=0;
allt = 0:stopsave*dt:tf;
allu = zeros(saveframes+1,N);
allv = zeros(saveframes+1,N);
allx = x;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\nL \t\t= %f\n',L);
fprintf('R \t\t= %f\n',R);
fprintf('dt \t\t= %f\n',dt);
fprintf('dx \t\t= %f\n',dx);
fprintf('N \t\t= %f\n',N);
fprintf('D1 \t\t= %f\n',D1);
fprintf('D2 \t\t= %f\n',D2);
fprintf('aa \t\t= %f\n',aa);
fprintf('bb \t\t= %f\n',bb);
fprintf('tf \t\t= %f\n',tf);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% GENERATING ICs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
u0 = aa+bb;
v0 = bb/(aa+bb)^2;
mode=16;
bell = exp(-0.00001*(x-(R-L-2)/2).^2);
uini=u0+ 0.1*u0*bell.*(0.0*rand(1,N)+cos(mode*pi*(x-L)/(R-L-dx)));
vini=v0-0.001*u0*bell.*(0.0*rand(1,N)+cos(mode*pi*(x-L)/(R-L-dx)));
uini=u0+0.1*u0*(rand(1,N)-0.45);
vini=v0+0.01*v0*(rand(1,N)-0.5);
save=0;
t=0;
u=uini;
v=vini;
% %No flux BCs
% u(1)=u(2);
% u(N)=u(N-1);
% v(1)=v(2);
% v(N)=v(N-1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot IC
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1); clf;
subplot(1,2,1)
set(gca,'FontSize',[16]);
plot(x,u,'LineWidth',3);
title('u(x,0)');
subplot(1,2,2)
set(gca,'FontSize',[16]);
plot(x,v,'LineWidth',3);
title('v(x,0)');
pause
save=save+1;
allu(save,:) = u;
allv(save,:) = v;
fprintf('\nt=%5.2f ',t);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Main loop
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for kk=1:maxstep
t=t+dt;
uN = u+dt*deriv_u(u,v,N);
vN = v+dt*deriv_v(u,v,N);
% %Fixed BCs
% uN(1)=u0;
% uN(N)=u0;
% vN(1)=v0;
% vN(N)=v0;
% %No flux BCs
% uN(1)=uN(2);
% uN(N)=uN(N-1);
% vN(1)=vN(2);
% vN(N)=vN(N-1);
u = uN;
v = vN;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% include noise
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(noise>0)
u = u+noise*(rand(N)-0.5);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display time on screen
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if (fix(kk/dispscreen)==kk/dispscreen)
fprintf('%5.2f ',t);
strlength=strlength+1;
if strlength>=10
strlength=0;
fprintf('\n');
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% store solution
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if (fix(kk/stopsave)==kk/stopsave)
save=save+1;
allu(save,:) = u;
allv(save,:) = v;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(fix(kk/stopdisp)==kk/stopdisp&kk>1)
figure(1); clf;
mytext=['T=(' num2str(kk*dt) ')/(' num2str(tf) ')'];
subplot(1,2,1)
set(gca,'FontSize',[16]);
plot(x,u,'LineWidth',3);
text(R-0.2*(R-L),min(u)+(max(u)-min(u))*1.1,mytext);
subplot(1,2,2)
set(gca,'FontSize',[16]);
plot(x,v,'LineWidth',3);
text(R-0.2*(R-L),min(v)+(max(v)-min(v))*1.1,mytext);
pause(0.0001)
end
end
fprintf('\n');
return;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%function rkZn = ode_rk(rkZ)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% global dx dxx oos u0 v0;
% global aa bb;
%
% k1 = dt*deriv_u(rkZ );
% k2 = dt*deriv_u(rkZ+0.5*k1);
% k3 = dt*deriv_u(rkZ+0.5*k2);
% k4 = dt*deriv_u(rkZ+ k3);
%
% rkZn = rkZ + oos*(k1 + 2*k2 + 2*k3 + k4);
%
% return;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function du = deriv_u(u,v,N)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global dx dxx D1 D2 oos u0 v0;
global aa bb;
%% Periodic BCs:
% up = [ u(1) u(1:N-1) ];
% um = [ u(2:N) u(N) ];
%% Fixed BCs
% up = [ u0 u(1:N-1) ];
% um = [ u(2:N) u0 ];
% No flux BCs
up = [ u(1) u(1:N-1) ];
um = [ u(2:N) u(N) ];
du = f(u,v) + D1*dxx*(up-2*u+um);
return;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dv = deriv_v(u,v,N)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global dx dxx D1 D2 oos u0 v0;
global aa bb;
%% Periodic BCs:
% vp = [ v(1) v(1:N-1) ];
% vm = [ v(2:N) v(N) ];
%% Fixed BCs
% vp = [ v0 v(1:N-1) ];
% vm = [ v(2:N) v0 ];
% No flux BCs
vp = [ v(1) v(1:N-1) ];
vm = [ v(2:N) v(N) ];
dv = g(u,v) + D2*dxx*(vp-2*v+vm);
return;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function fuv = f(u,v)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global dx dxx D1 D2 oos u0 v0;
global aa bb;
fuv = aa-u+u.^2.*v;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function guv = g(u,v)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global dx dxx D1 D2 oos u0 v0;
global aa bb;
guv = bb-u.^2.*v;