%
% Ricardo Carretero 20/Oct/2003
% Integrates a reaction diffusion eq in 1D using finite differences
%
% u_t = f(u) + D u_xx
%
% for the Spruce Budworm model:
%
% f(u)=r*u*(1-u/q) -u^2/(1+u^2)
%
% Ex: SpruceBudworm(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=10;L=0;R=4;snaps=100;disp=100;D=1;
%N=100;dt=0.0005;tf=10;L=0;R=6;snaps=100;disp=100;D=1;
%SpruceBudworm(N,dt,tf,L,R,disp,snaps)
%wave_plotsurf(allt,allx,allu',R,1,N+1,'no')
%
function [] = SpruceBudworm(N,dt,tf,L,R,exposh,saveframes)
global allt allx allu;
global dx dxx D oos r q allt allx allu;
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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
r = pi^2/25;
q = 10;
Lc = sqrt(D)*pi/sqrt(r)
fprintf('\nr \t\t= %f\n',r);
fprintf('q \t\t= %f\n',q);
fprintf('R-L \t\t= %f\n',R-L);
fprintf('Lc \t\t= %f\n',Lc);
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).^2).*(0.1*rand(1,N)+1+sin(20*pi*(x-L-0.1)/(R-L)));
%uini=0.01*(0.0*rand(1,N)+sin(2*pi*(x-L)/(R-L)));
fp=3.8853;
fp=0.45476301375605;
fp=5.66054422717771;
fp=0;
uini=1*uini+fp;
uini(1)=fp;
uini(N)=fp;
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);
u = uN;
% u(1)=fp;
% u(N)=fp;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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 r q;
%% Periodic BCs:
% up = [ u(1) u(1:N-1) ];
% um = [ u(2:N) u(N) ];
% Fixed BCs
up = [ 2*u(1)-u(2) u(1:N-1) ];
um = [ u(2:N) 2*u(N)-u(N-1) ];
du = f(u) + D*dxx*(up-2*u+um);
return;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function fu = f(u)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global dx dxx D oos r q;
fu = r*u.*(1-u/q) -u.^2./(1+u.^2);