%
%
% Ricardo Carretero 20/Oct/2003
% Integrates a reaction diffusion eq in 2D using finite differences
%
% u_t = f(u) + D (u_xx+u_yy)
%
% Ex: PDE_2D(Nx,Ny,dt,tf,Lx,Rx,Ly,Ry,disp,savesnaps)
%
% npts=N, dt:time step, tf:t_final, domain:[L,R]
%
%Ex:
%global D allt allx ally allu ran;
%Nx=20;Ny=20;dt=0.0005;tf=5;Lx=0;Rx=pi+1;Ly=0;Ry=pi+1;snaps=100;disp=100;D=1;
%PDE_2D(Nx,Ny,dt,tf,Lx,Rx,Ly,Ry,disp,snaps)
%
function [] = PDE_2D(Nx,Ny,dt,tf,Lx,Rx,Ly,Ry,exposh,saveframes)
global allt allx ally allu;
allx=[]; ally=[]; allt=[]; allu=[];
ran=1;
noise=0;
%This next line is here so I do not forget to globalize
if(isempty(ran))
fprintf('\n');
error('Error: first do: global global dx dxx dy dyy D oos allt allx ally allu ran;')
end
if(Nx~=Ny)
fprintf('\n');
error('You need Nx=Ny, Nx not= Ny has not been implemented yet')
end
maxstep=tf/dt;
numdispscreen=100;
dispscreen= fix(maxstep/numdispscreen);
if(dispscreen==0) dispscreen=1;end;
if(exposh>0) stopdisp= fix(maxstep/exposh);
else stopdisp=9e+99;
end
if(saveframes>0)stopsave= fix(maxstep/saveframes);
else stopsave=9e+99;
end
if(stopsave==0) stopsave=9e+99;end;
% Spatial computational domain defined
x=linspace(-Lx,Rx,Nx);
y=linspace(-Ly,Ry,Ny);
X=kron(x,ones(Ny,1));
Y=kron(y',ones(1,Nx));
dx=x(2)-x(1);
dy=y(2)-y(1);
dxx=1/dx^2;
dyy=1/dy^2;
oos=1/6;
rand('state',0) %this fixes the seed for the random generation
strlength=0;
ifig=5;
allx = x;
ally = y;
allt = 0:stopsave*dt:tf;
allt = allt;
allu = zeros(saveframes+1,Ny,Nx);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% GENERATING ICs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% uini= 0.01*exp(-(X-(Rx-Lx-2)/2).^2).*(0.1*rand(Ny,Nx)+1+sin(20*pi*(X-Lx-0.1)/(Rx-Lx)));
% uini=uini.*0.01*exp(-(Y-(Ry-Ly-2)/2).^2).*(0.1*rand(Ny,Nx)+1+sin(20*pi*(Y-Ly-0.1)/(Ry-Ly)));
uini = 0.01*(0.1*rand(Ny,Nx)+sin(3*pi*(X-Lx)/(Rx-Lx)));
uini = uini .*sin(2*pi*(Y-Ly)/(Ry-Ly));
uini=uini+0*0.001;
uini(1,1:Nx)=0;
uini(Ny,1:Nx)=0;
uini(1:Ny,1)=0;
uini(1:Ny,Nx)=0;
save=0;
t=0;
u=uini;
allu(1,:,:) = u;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% setup
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
save=1;
kk=0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('Lx \t\t= %f\n',Lx);
fprintf('Rx \t\t= %f\n',Rx);
fprintf('Ly \t\t= %f\n',Ly);
fprintf('Ry \t\t= %f\n',Ry);
fprintf('dt \t\t= %f\n',dt);
fprintf('dx \t\t= %f\n',dx);
fprintf('dy \t\t= %f\n',dy);
fprintf('Nx \t\t= %f\n',Nx);
fprintf('Ny \t\t= %f\n',Ny);
fprintf('tf \t\t= %f\n',tf);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot initial data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
high = max(max(u));
low = min(min(u));
viewUtop=[0 90];
viewU=[-38 30];
figure(1+ifig); clf;
set(gca,'FontSize',[16]);
zl='u(x,y,t)';title='';
mytext=['T=(' num2str(kk*dt) ')/(' num2str(tf) ')'];
plotopU='mag'; axU=[Lx Rx Ly Ry low high];
plotU(x,y,u,axU,zl,title,mytext,plotopU);
view(viewU);
figure(2+ifig); clf;
set(gca,'FontSize',[16]);
zl='u(x,y,t)';title='';
mytext=['T=(' num2str(kk*dt) ')/(' num2str(tf) ')'];
plotopU='mag'; axU=[Lx Rx Ly Ry low high];
plotU(x,y,u,axU,zl,title,mytext,plotopU);
view(viewU);
figure(3+ifig); clf;
set(gca,'FontSize',[16]);
zltop='u(x,y,t)';title='';
plotopUtop='height'; axUtop=[Lx Rx Ly Ry low-5*(high-low) high+5*(high-low)];
axis(axUtop);
plotU(x,y,u,[],zltop,title,mytext,plotopUtop);
view(viewUtop);
colorbar
pause(0.001);
% fprintf('Press any key...');pause;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Main loop
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for kk=1:maxstep
t=kk*dt;
u = rk4_2d(u,Nx,Ny,dt,oos,dxx,dyy);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% include noise
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(noise>0)
u = u+noise*(rand(Ny,Nx)-0.5);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display time on screen
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if (fix(kk/dispscreen)==kk/dispscreen)
fprintf('%5.2f ',kk*dt);
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;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(fix(kk/stopdisp)==kk/stopdisp&kk>1)
figure(1+ifig); clf;
set(gca,'FontSize',[16]);
mytext=['T=(' num2str(kk*dt) ')/(' num2str(tf) ')'];
plotU(x,y,u,[],zl,title,mytext,plotopU);
view(viewU);
figure(2+ifig); clf;
set(gca,'FontSize',[16]);
mytext=['T=(' num2str(kk*dt) ')/(' num2str(tf) ')'];
plotU(x,y,u,axU,zl,title,mytext,plotopU);
view(viewU);
figure(3+ifig); clf;
set(gca,'FontSize',[16]);
plotU(x,y,u,[],zltop,title,mytext,plotopUtop);
view(viewUtop);
% colorbar
pause(0.001);
end
end
figure(3+ifig);colorbar
fprintf('\n');
return;