diffusion equation with crank nickolson method
8 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
syrine hidri
am 13 Mai 2019
Kommentiert: Bjorn Gustavsson
am 14 Mai 2019
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
0 Kommentare
Akzeptierte Antwort
Bjorn Gustavsson
am 13 Mai 2019
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
2 Kommentare
Bjorn Gustavsson
am 14 Mai 2019
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.
Weitere Antworten (0)
Siehe auch
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!