MATLAB: Diffusion equation with crank nickolson method

#equation #diffusion #crank #nicolson #pde #1d

I was solving a diffusion equation with crak nickolson method
the boundry conditons are :
I think i have a problem in my code because i know that ∆n(t) for a constant x should be a decreasing exponential curve but i found this
% Parameters
D=0.054;  alpha=41000; taux=300e-9; 
L=270e-9;Tmax=5e-7;
nx=6; nt=11;
dx=L/nx; dt=Tmax/nt;
%constant
a = (dt*D)/2*(dx*dx); 
b=1+2*a+(dt/(2*taux));
c=1-2*a-(dt/(2*taux)); 
%bluid matrix
M=zeros(nx-2,nx-2);
A=(D*dt)/(2*(dx*dx));
main_diag = (1+2*A+(dt/(2*taux)))*ones(nx-2,1);
aux_diag =(-A)*ones(nx-2,1);
M = full(spdiags([aux_diag ,main_diag,aux_diag],[-1,0,1], nx-2, nx-2));
for i=1:nt
    t(i)=i*dt;  
end 
% Boundary conditions and Initial Conditions 1
for i=1:nx
    x(i)=i*dx
   deltan0 (i)= 1e10*exp(-alpha*x(i));
end 
ligne=1
deltan(ligne,:)=deltan0; %store the first lign of final matrix with the boundry condition 1
% Loop over time
for k=2:nt
   for i=1:nx-2
         d(i)=a*deltan0(i)+c*deltan0(i+1)+a*deltan0(i+2);
   end  % note that d is row vector
 deltan_curr=M\d';
 deltan0=[1 deltan_curr' 0]; % to start over
 deltan(k,:)=deltan0;% Store results for future use
end
for i=1:nt
    c(i)=deltan(i,2);
end
loglog(t,c);
I used to work with same method with this pdf https://matlabgeeks.weebly.com/uploads/8/0/4/8/8048228/crank-example_with_matlab_code-v3__doc_.pdf
please help me

Best Answer

That doesn't look like C-N. Since C-N uses the values at t_i + dt/2 to calculate the gradients you end up with a system of equations for \Delta n with two matrices:
Mlhs * deltan_{t+1} = Mrhs*deltan_{t} + Q(t+1/2) % Loose notation...
You should be fine implementing your solution straight from: Crank-Nicholson at wikipedia, check that you correctly handle the boundary conditions, I couldnt read the code as typed in so, you should consider editing your question to make your code show up as code.
HTH
Related Question