%
% Ricardo Carretero 20/Oct/2003
% Integrates a reaction diffusion eq in 1D using finite differences
%
% u_t = f(u) + D u_xx
%
% Ex: PDE_integrator(N,dt,tf,L,R,disp,savesnaps)
%
% npts=N, dt:time step, tf:t_final, domain:[L,R]
%
%Ex:
%global dx dxx D oos allt allx allu;
%N=100;dt=0.0005;tf=20;L=0;R=10;snaps=100;disp=100;D=1;
%PDE_integrator(N,dt,tf,L,R,disp,snaps)
%wave_plotsurf(allt,allx,allu',R,1,N+1,'no')
%
function [] = PDE_integrator(N,dt,tf,L,R,exposh,saveframes)
global allt allx allu;
global dx dxx oos;
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);
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('tf \t\t= %f\n',tf);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% GENERATING ICs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%uini=0.01*exp(-(x-(R-L)/2).^2).*(0.1*rand(1,N)+1+sin(20*pi*(x-L)/(R-L)));
uini=0.01*(0.0*rand(1,N)+sin(2*pi*(x-L)/(R-L)));
save=0;
t=0;
u=uini;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot IC
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1); clf;
set(gca,'FontSize',[16]);
plot(x,u,'LineWidth',3);
title('u(x,0)');
pause
save=save+1;
allu(save,:) = u;
fprintf('\nt=%5.2f ',t);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Main loop
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for kk=1:maxstep
t=t+dt;
uN = u+dt*deriv_u(u,N);
%Fixed BCs
uN(1)=0;
uN(N)=0;
u = uN;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(fix(kk/stopdisp)==kk/stopdisp&kk>1)
figure(1); clf;
set(gca,'FontSize',[16]);
mytext=['T=(' num2str(kk*dt) ')/(' num2str(tf) ')'];
plot(x,u,'LineWidth',3);
text(R-0.2*(R-L),max(u)*0.9,mytext);
pause(0.0001)
end
end
fprintf('\n');
return;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%function rkZn = ode_rk(rkZ)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% global dx dxx oos;
%
% 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,N)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global dx dxx D oos;
%% Periodic BCs:
% up = [ u(1) u(1:N-1) ];
% um = [ u(2:N) u(N) ];
% Fixed BCs
up = [ 0 u(1:N-1) ];
um = [ u(2:N) 0 ];
du = f(u) + D*dxx*(up-2*u+um);
return;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function fu = f(u)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global dx dxx D oos;
fu = 1*u;