diffusion equation with crank nickolson method

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

Answer 1

[Bjorn Gustavsson]

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

[syrine hidri]

I edit my question it's shown as code now ,you can check it

i have also worked with an another solution ,if you can verify both of them i

i need a help

%parameters
 D=0.054;
taux=300e-9;
Tmax=2e-7;
N=10;
L=270e-9;
dt=Tmax/N;
dx=L/N;

format short e

for i=1:N+1 
    t(i)=(i-1)*dt;
    x(i)=(i-1)*dx;
end 

%bluid the matrix M
cst=(dt*D)/2*(dx*dx);
M=zeros(N+1,N+1);
main_diag = (-4*cst/dt)-(1/taux)*ones(N+1,1);
aux_diag =D/(dx*dx)*ones(N+1, 1 )
M = full(spdiags([aux_diag ,main_diag,aux_diag],[-1,0,1], N+1, N+1)


deltan_out = zeros(N+1,N+1);

%generate iteration matrice
%A=I-0.5*dt*M 
%B*I+0.5*dt*M
A=full(speye(N+1))-0.5*dt*M;
B=full(speye(N+1))+0.5*dt*M;



curr_slice = 1;
deltan_curr=zeros(N+1,1);

%first boundary condition
for i=1:N+1
   deltan_curr (i)= exp(-alpha*x(i));
end 
deltan_out(:,curr_slice)= deltan_curr;
%store the first bondary condition in the final matrix

for i= 1 : N

      deltan_curr= inv(A)*B*deltan_curr;
      curr_slice = curr_slice +1;
      deltan_out(:,curr_slice)= deltan_curr ;   
end

%second boundary condition
for i= 1 : N+1
deltan_out(N+1,i)= 0 ; 
end
for i=1:N+1
c(i)=deltan_out(1,i);
end
 plot (t,c);

[Bjorn Gustavsson]

Make yourself a small (as in number of spatial grid-points) version of your problem, step yourself through the problem time-step by time-step. Have a look at how your matrix

inv(A)*B

looks, have a look at how your boundary conditions are respected by your solution. Typically for fixed boundary values you have to set the first an last rows of A and B to zeros except the diagonal elements that is one. Just a couple of hints, haven't got time for more detailed help.