summerlife163
2009-04-13, 21:01
function [c,M,cf]=befordif(D,K,dt,dx,h,fo,bi,N,A,V,k1,k2,AA)
D=3.2e-9;K=6665;dt=1e-3;dx=1e-3;
h=3.21;
fo=D*dt/dx.^2;bi=h*dx/D;
N=94.9;
A=16.994;
V=A*2.5;
M(200001)=0;
k1=0.97;k2=0.19;
AA=62.404;
c(1:9,1)=3.19e7;
cb(200001)=0;
cf(200001)=0;
for i=2:200001
c(1,i)=2*fo*(c(2,i-1)+bi*cf(i-1))+(1-2*bi*fo/K-2*fo)*c(1,i-1);
cb(i-1)=c(1,i)/K;
for x=2:8
c(x,i)=fo*(c(x-1,i-1)+c(x+1,i-1))+(1-2*fo)*c(x,i-1);
end
c(9,i)=2*fo*c(8,i-1)+(1-2*fo)*c(9,i-1);
M(i)=dt*k1*cf(i-1)+(1-dt*k2)*M(i-1);
cf(i)=dt/V*(h*A*cb(i-1)+k2*AA*M(i-1))+(1-dt/V*(h*A+N*V+k1*AA))*cf(i-1);
end
D=3.2e-9;K=6665;dt=1e-3;dx=1e-3;
h=3.21;
fo=D*dt/dx.^2;bi=h*dx/D;
N=94.9;
A=16.994;
V=A*2.5;
M(200001)=0;
k1=0.97;k2=0.19;
AA=62.404;
c(1:9,1)=3.19e7;
cb(200001)=0;
cf(200001)=0;
for i=2:200001
c(1,i)=2*fo*(c(2,i-1)+bi*cf(i-1))+(1-2*bi*fo/K-2*fo)*c(1,i-1);
cb(i-1)=c(1,i)/K;
for x=2:8
c(x,i)=fo*(c(x-1,i-1)+c(x+1,i-1))+(1-2*fo)*c(x,i-1);
end
c(9,i)=2*fo*c(8,i-1)+(1-2*fo)*c(9,i-1);
M(i)=dt*k1*cf(i-1)+(1-dt*k2)*M(i-1);
cf(i)=dt/V*(h*A*cb(i-1)+k2*AA*M(i-1))+(1-dt/V*(h*A+N*V+k1*AA))*cf(i-1);
end