diffusion equation with crank nickolson method
8 ビュー (過去 30 日間)
古いコメントを表示
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 件のコメント
採用された回答
Bjorn Gustavsson
2019 年 5 月 13 日
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 件のコメント
Bjorn Gustavsson
2019 年 5 月 14 日
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.
その他の回答 (0 件)
参考
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!