FDM and impulse response for diffusion equation not yielding similiar results

1 visualizzazione (ultimi 30 giorni)
I am implementing forward euler and impulse response solutions for fick's diffusion equation as part of a larger project.
There's a constant difference between the results that seems more than quantization errors though i would be happy to be corrected.
I've maintained c<1/2 for forward euler stability and have tried playing around with x,dx and t,dt but the resulting difference is still the same shape. Thank you for your time.
x0 = 0.0;
xEnd = 1.0;
t0 = 0.0;
tEnd = 0.1;
dx = 0.001;
dt = 0.0001;
D = 0.00001;
% spatial step for 1d
c = D * dt/(dx^2); % stability requires c < 1/2 for 1d diffusion
%discretization:
xSpace = x0:dx:xEnd;
tSpace = t0:dt:tEnd;
volume = (xEnd - x0)^1;
Nx = numel(xSpace);
Nt= numel(tSpace);
%% initial condition:
initialMass = 100.0;
initialConc = initialMass/volume;
sourceDiscreteLocation = ceil(Nx/2);
% spatial step for 1d
% Initialize spatial profile
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

Risposta accettata

Or Bovan
Or Bovan il 22 Mag 2020
This is a question of both increasing dx while maintaining stability thus one should be mind full of dx to dt ratio as per euler stability. accuracy seems to improve on certain ratios but i've no scientific deduction beyond trial and error.
Answered myself.

Più risposte (0)

Prodotti

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by