x0 = 0.0;
xEnd = 1.0;
t0 = 0.0;
tEnd = 0.1;
dx = 0.001;
dt = 0.0001;
D = 0.00001;
c = D * dt/(dx^2);
xSpace = x0:dx:xEnd;
tSpace = t0:dt:tEnd;
volume = (xEnd - x0)^1;
Nx = numel(xSpace);
Nt= numel(tSpace);
initialMass = 100.0;
initialConc = initialMass/volume;
sourceDiscreteLocation = ceil(Nx/2);
sourceMap=zeros(1,Nx);
evenConcMap=zeros(1,Nx);
oddConcMap=zeros(1,Nx);
saveMap=zeros(1,Nx);
ImpulseResponseMap=zeros(1,Nx);
exponent = zeros(1,Nx);
saveInterval = 50;
sourceMap(sourceDiscreteLocation)=initialConc;
evenConcMap(1,:) = sourceMap;
for t = 1:Nt
if(mod(t,2)==0)
evenConcMap(2:Nx-1)=oddConcMap(2:Nx-1)+c*(-2*oddConcMap(2:Nx-1)+oddConcMap(3:Nx)+oddConcMap(1:Nx-2));
evenConcMap(Nx)=0;
evenConcMap(1)=0;
else
oddConcMap(2:Nx-1)=evenConcMap(2:Nx-1)+c*(-2*evenConcMap(2:Nx-1)+evenConcMap(3:Nx)+evenConcMap(1:Nx-2));
oddConcMap(Nx)=0;
oddConcMap(1)=0;
end
if(mod(t,saveInterval) == 0)&&(mod(saveInterval,2) == 0)
saveMap=evenConcMap;
end
if(mod(t,saveInterval) == 0)
saveMap=oddConcMap;
end
end
if(mod(Nt,2) == 0)
result = evenConcMap;
else
result = oddConcMap;
end
sourcePhysicalLocation = (Nx/2)*dx;
alpha = 4*D*tEnd;
exponent = exp(-((xSpace - sourcePhysicalLocation).^(2.0))./alpha);
ImpulseResponseMap=c*initialConc*(exponent)./((pi*alpha)^((1.0)/(2.0)));
ImpulseResponseMap(1)=0;
ImpulseResponseMap(Nx)=0;
figure;
hold on
semilogx(xSpace,result);
semilogx(xSpace,ImpulseResponseMap);
title('Diffusion vs Impulse response 1d when source at center');
legend('diffusion','impulse response')
xlabel('x - semilog scale');
ylabel('concentration(x,tEnd)');
hold off